Assessing Periodicity of Periodic Leg Movements During Sleep

Background: Periodic leg movements (PLM) during sleep consist of involuntary periodic movements of the lower extremities. The debated functional relevance of PLM during sleep is based on correlation of clinical parameters with the PLM index (PLMI). However, periodicity in movements may not be reflected best by the PLMI. Here, an approach novel to the field of sleep research is used to reveal intrinsic periodicity in inter movement intervals (IMI) in patients with PLM. Methods: Three patient groups of 10 patients showing PLM with OSA (group 1), PLM without OSA or RLS (group 2) and PLM with RLS (group 3) are considered. Applying the “unfolding” procedure, a method developed in statistical physics, enhances or even reveals intrinsic periodicity of PLM. The degree of periodicity of PLM is assessed by fitting one-parameter distributions to the unfolded IMI distributions. Finally, it is investigated whether the shape of the IMI distributions allows to separate patients into different groups. Results: Despite applying the unfolding procedure, periodicity is neither homogeneous within nor considerably different between the three clinically defined groups. Data-driven clustering reveals more homogeneous and better separated clusters. However, they consist of patients with heterogeneous demographic data and comorbidities, including RLS and OSA. Conclusions: The unfolding procedure may be necessary to enhance or reveal periodicity. Thus this method is proposed as a pre-processing step before analyzing PLM statistically. Data-driven clustering yields much more reasonable results when applied to the unfolded IMI distributions than to the original data. Despite this effort no correlation between the degree of periodicity and demographic data or comorbidities is found. However, there are indications that the nature of the periodicity might be determined by long-range interactions between LM of patients with PLM and OSA.

constitutes periodicity may have a significant impact on the final result. According to the new standards PLM occur in a series of at least four LM and are separated by 5-90 s (Zucconi et al., 2006). An increased number of PLM during sleep with a PLM index (PLMI, i.e., the number of PLM per hour of total sleep time) of five or more is considered abnormal (American Academy of Sleep Medicine, 2005). Though the PLMI is taken to reflect periodicity in PLM in most studies, periodicity is not always obvious and PLMI may fail to report it. The definition of a period implicit in the PLMI encompasses inter movement intervals (IMI) from 5 to 90 s. This range of IMI may not be sufficiently restrictive to identify periodicity in LM and could allow randomly occurring LM to be considered periodic. On the other hand, IMI may change across the night (Coleman et al., 1980;Pollmächer and Schulz, 1993;Nicolas et al., 1998;Ferri et al., 2009) and PLM will be missed if period criteria are too narrow. Varying IMIs during night may lead to broad IMI distributions or to a distribution showing several local maxima. This makes a wide range of periods appear similarly prevalent, and thus hides intrinsic periodicity. Another way intrinsic periodicity may become difficult to detect is when pooling data from patients with different average periods of PLM. The definition of periodicity of LM is problematic and may not be sensitive enough.

IntroductIon
Periodic leg movements (PLM) during sleep are a sleep related phenomenon with periodic episodes of repetitive stereotypical movements of the lower extremities (American Sleep Disorder Association, 1993). PLM are more frequent in several sleep disorders, including restless legs syndrome (RLS, Montplaisir et al., 1997), narcolepsy (Moscovitch et al., 1993), rapid eye movement sleep behavior disorder (Fantini et al., 2002) as well as in disorders not primarily affecting sleep such as end stage renal disease (Benz et al., 2000). PLM during sleep are also documented in normal subjects (Pennestri et al., 2006) giving rise to the debate about functional relevance.
To determine the functional significance of PLM, the criteria used to score PLM during sleep are pivotal. Thus the criteria of detecting single leg movements (LM) as well as determining what 2. Does periodicity differentiate between the clinical conditions of RLS, of OSA and of PLM without RLS or OSA? 3. Does data-driven clustering separate patients into different groups and are such groups clinically homogeneous?
Answering these questions could potentially provide further information about the clinical relevance of PLM during sleep.

PatIents and settIng
Thirty patients were retrospectively evaluated at the Inselspital Bern Sleep Laboratory, all having a PLMI >15/h. Ten patients were diagnosed with PLM and OSA (group 1) with 10 or more obstructive apneas or hypopneas per hour of total sleep time. Ten patients were diagnosed with PLM without RLS or OSA (group 2). Five of these patients fulfilled the diagnostic criteria for PLMD (American Academy of Sleep Medicine, 2005). The last 10 patients had in addition RLS (group 3), fulfilling the minimal criteria accepted for the diagnosis of RLS . The presence of the following medical conditions was recorded: chronic cardiac or pulmonary disease, snoring, systemic hypertension, depression, obesity, Parkinsonism, headaches, fatigue, epilepsy, non-rapid eye movement parasomnia, idiopathic hypersomnia, and insomnia. The case history analysis revealed that 17 of 30 patients were taking at least one medication on a regular basis at the time of evaluation. In group 1 (PLM with OSA) 2 patients out of 10 were not taking any medication, 2 patients were treated with opioids, 1 patient was taking betablocker and 1 patient paroxetin. The remaining 4 patients were taking medication not known to influence PLM. In group 2 (PLM without OSA or RLS) 6 patients out of 10 were not taking any medication, 1 patient was taking dopamine, 1 patient venlafaxine, trimipramin and betablocker and 1 patient mirtazapine and olanzapine. In group 3 (PLM with RLS) 5 patients out of 10 were not taking any medication, 1 patient was taking gabapentin and 1 patient was treated with mirtazapine and trazodone. The remaining 3 patients were taking medication not known to influence PLM. Tricyclic antidepressants, selective serotonin reuptake inhibitors and betablockers are considered to exacerbate PLM, whereas dopamine agonists, opioids, gabapentin, apomorphine are considered to reduce PLM .
Characteristics of patient groups are shown in Table 1 and Figure 1. The studied cohort shows a broad distribution of demographic data and sleep parameters. To assess statistical differences between the groups first a Kruskal-Wallis H-test for different central tendencies in the three groups (Siegel, 1956) was performed. If the differences were significant pairwise Mann-Whitney-Wilcoxon U-tests (Siegel, 1956) were subsequently performed in order to reveal where medians are different. Both tests were performed on significance level α = 0.01. For the Epworth Sleepiness Scale score (ESS) the H-test was not significant and for the PI it was only marginally significant (p = 0.049). With a few exceptions (see Figure 1) most pairwise U-tests were not significant. The reason is that the within group variability of the small groups (N = 10) was rather large and overshadowed the between group differences in most cases. However, these characteristics are representative for a typical patient cohort in a sleep clinic. Note that the highly significant differences with respect to AHI are due to the clinical a priori definition of the groups. This lack of sensitivity could be one of the reasons for the ongoing debate about their functional relevance (Ancoli-Israel et al., 1991;Mendelson, 1996;Hornyak et al., 2004;Carrier et al., 2005;Högl, 2007;Mahowald, 2007).
A new statistical approach to analyze PLM was developed by Ferri et al. (2006b) who recommend to study the entire distributions of IMI rather than parametric measures as mean and standard deviation. They examined periodicity with a metric, termed the Periodicity Index (PI) measuring the percentage of LM appearing in series of periodic events. Instead of focusing on IMI frequency distributions the PI accounts for certain dependencies between subsequent IMIs. The PI was investigated for homogenous patient groups with RLS, narcolepsy/cataplexy and rapid eye movement sleep behavior disorder and without further potentially sleep relevant comorbidities or medication (Ferri et al., 2006a,b;Manconi et al., 2007). Based on different levels of the PI patients were divided into subgroups. For patients with narcolepsy/cataplexy and RLS a significantly different median of the PI was found (Ferri et al., 2006a). Very recently increased PI was also detected in a subgroup of patients with insomnia, for which no apparent cause could be found and enabled diagnosis of periodic leg movements disorder (PLMD, Ferri et al., 2009). Quantifiable aspects of PLM such as periodicity may serve to differentiate between clinical conditions and may eventually help capturing the nature of the underlying biological phenomenon (Rye, 2006). However, as afore mentioned, quantitative assessment of the periodicity of LM/PLM as measured by IMI is difficult because the empirical probability distribution may obscure underlying periodicity.
To date no study uses methods to enhance and reveal intrinsic periodicity before analyzing the statistics of IMI. We here employ the so-called unfolding procedure, which potentially uncovers intrinsic periodicity of LM data that is confounded by external influences. The unfolding procedure was originally developed in a branch of statistical physics called Random Matrix Theory (RMT) to study and compare the correlation between energy levels of heavy atomic nuclei and other complex quantum systems (see Brody et al., 1981;Guhr et al., 1998;Mehta, 2004, for comprehensive reviews). In the context of neurophysiology the concept was recently applied to the eigenvalue spectrum of the cross-correlation matrix calculated from multichannel EEG (Seba, 2003;Müller et al., 2006;Baier et al., 2007).
In the present study we use the unfolding concept to assess the intrinsic periodicity of IMI of PLM examining a typical clinical patient cohort with PLM during sleep and different comorbidities from our database. Instead of studying the absolute IMI sizes, IMI are locally normalized to unit mean event spacing. Under these conditions the distribution of the spacings of randomly distributed events decays exponentially (Poisson distribution). In contrast, for strictly periodic systems the events are equally spaced. Fitting intermediate distributions, which interpolate between these extreme cases, to unfolded IMI distributions we quantify PLM periodicity.
Specifically, the goal of our study is to answer the following questions: 1. Does applying the unfolding procedure enhance and reveal periodicity that is not obvious from the measured IMI distribution? Rummel et al. Assessing periodicity of PLMS

PolysoMnograPhy
Each patient underwent a diagnostic full night polysomnography (PSG). Patients had given their written informed consent for using their data for scientific purposes. Retrospective use of clinical data was approved by the local ethic committee. Parameters recorded were electroencephalogram (EEG, minimal leads C4-M1, C3-M2, O2-M1, O1-M2, and additional F4-M1, F3-M2) respiration (nasal and oral airflow with thermistors, nasal pressure with a cannula, thoracic and abdominal respiratory movements with strain gauges, and oxygen saturation with finger oximetry), left and right electrooculogram, and three electromyogram (EMG) channels (submentalis muscle and right and left tibialis anterior muscles). EMG signals were sampled at 200 or 500 Hz, high pass filtered at 10 Hz and low-pass filtered at 45 or 100 Hz, respectively. Before the beginning of a recording a sleep technician checked that the amplitude of the EMG signal from the pendulum remains a periodic system. Solely, its intrinsic periodicity is obscured by an outer influence [its varying length l(T)] such that it can hardly be detected by measuring the spacing distribution P 1 (S) in the most straight forward way.

Formal development
Instead one has to make use of a method designed for revealing and quantifying intrinsic fluctuations regardless of systematic trends of possibly unknown nature and origin. Originally it was introduced in the context of nuclear physics to put energy spectra of different atomic nuclei on the same footing before making comparison of certain fluctuation properties. Nowadays this method is extensively used in the field of quantum chaos. The crucial point consists in performing the so-called "unfolding procedure" (Brody et al., 1981;Haake, 1992;Guhr et al., 1998;Mehta, 2004) where the event times are transformed T i → t i such that the average distance (or time interval) between adjacent events becomes unity. Whereas a simple normalization can attain this goal on the global level, unfolding is designed to assure unit average spacing also locally. Different to the spacings S i of the directly measured event times T i the spacings s i = t i+1 − t i of the unfolded t i can be used to reveal intrinsic periodicities.
To develop the unfolding procedure formally, one starts from the density of n discrete events -in our case LM times: Here δ(·) denotes the delta-distribution, i.e., a function that is very narrowly peaked at zero with unit integral. In the second identity we have split the event density into a smoothly T-dependent part that will subserve the unfolding distribution and a fluctuating residual. It is obvious that the probability P 1 (S) of finding a neighboring event at a distance S from a selected one crucially depends on the event density r(T). Neighboring events will in general be closer if the event density is high than for low event density. To eliminate this influence, which may affect the inference of the system's periodicity, one has to correct all spacings for the local smooth event density r smooth (T). Estimates of the event density r(T) from a finite amount n of data typically involve parameters like bin width and bin number. To avoid such influences it is more convenient to work with the parameter free cumulated event density: that counts the number of events in the interval [ −∞,T], see Figure 2A for an illustration. At every event time T i the so-called "staircase function" makes one step up. The unfolded events t i are obtained from the smooth part of the cumulated event density by the transformation: By definition the mean spacing 〈s〉 of the t i is locally equal to one along the whole unfolded event series. The crucial issue of unfolding an event series is to find an appropriate estimate of N smooth (T). In some physical cases an analytical formula for r smooth (T) or N smooth (T) is known and can two tibialis anterior muscles was below 2 μV at rest and impedance was kept less than 5 kΩ. Sleep stages were scored by an experienced scorer (H. Gast) visually in 30 s epochs according to Rechtschaffen and Kales (1968).

leg MoveMent detectIon and coMPutatIon of IndIces
Leg movements were first detected by the software Somnologica, Version 5.0.1. (Embla N7000 Recording Systems, Embla, Broomfield, CO, USA). The detections proposed by the automatic analysis were then edited by an experienced scorer (H. Gast) before computing the various parameters. All LM during sleep were included irrespective of association with arousals or apneas or hypopneas. LM while awake were not considered and thus a PLM sequence was terminated by intervening wakefulness. Otherwise the WASM standards for recording and scoring PLM during sleep were applied (Zucconi et al., 2006). A LM duration was at least 0.5 s and no longer than 10 s and IMI were calculated as the time elapsed from onset to onset of subsequent movements. LMs in the right and left leg were considered as simultaneous LMs and counted as only one LM if the separation between termination of the earlier and onset of the latter was equal or less than 0.5 s. LM were included in periodicity analysis if they occurred in a series of four or more events and if they were separated by IMI of more than 5 s and equal or less than 90 s.

the unfoldIng transforMatIon
To characterize and quantify the intrinsic periodicity of PLM regardless of confounding systematic influences the "unfolding procedure," which originates from RMT, was introduced to the field of PLM analysis for the first time. In the following a simple example is given, illustrating why unfolding may become necessary to reveal intrinsic periodicity. Thereafter the unfolding transformation is developed formally. Finally its usefulness is illustrated by applying it to model data.

Motivation
The ideal pendulum often serves as the paradigm of a periodic system. Event times T i (i = 1,…,n) are defined as whenever the pendulum crosses the equilibrium from the left to the right. In the absence of any disturbing influences the events are equally spaced in time: T i+1 = T i + S. Here S = ∆T is the fundamental period of the pendulum which is related to the frequency f and the length l of the pendulum by ∆T f l = ∼ 1/ . According to this equation a shorter pendulum l′ < l has a higher frequency f′ and a smaller period ∆T′. If the pendulum is disturbed by environmental influences of any kind the spacings S i = T i+1 − T i are no longer strictly equal but follow a distribution P 1 (S), which we assume centered at ∆T. The characteristic of a (quasi-)periodic system is the narrowness of this distribution P 1 (S). In contrast, broad distributions P 1 (S) are identified with aperiodic systems where a characteristic frequency is absent.
To extend the example, suppose that the length of the pendulum becomes time dependent: l = l(T). If the length of the pendulum suffers a sudden change from l to l′ < l the distribution P 1 (S) becomes bimodal with two centers at ∆T and ∆T′. If, however, a continuous but slow variation of l takes place over a large time period a broad distribution may result for P 1 (S). Without doubt, the Here we proceed as follows, see Figure 2: according to the definition of PLM during sleep (Zucconi et al., 2006) we select LM runs of n ≥ 4 events occurring with IMI between 5 and 90 s from the whole series of LM times. An odd low order polynomial is fitted independently to each of these pieces. In order to prevent overfitting the staircase function with a polynomial of too high degree m, this parameter is given a logarithmic dependence on n: Here the floor function b·c denotes the largest integer smaller than the argument. Eq. 4 guarantees that 1 ≤ m < n. It must be checked be used in Eq. 3. Otherwise N smooth (T) must be approximated from the experimentally measured data. To this end different strategies can be used, see e.g., discussion by Bruus and Anglés d'Auriac (1997). A common procedure consists of the so-called Gaussian broadening, where the delta-distributions in Eq. 1 are replaced by Gaussians. The constant or event-dependent widths are used as fit parameters. It has been found that Gaussian broadening tends to artificially introduce periodicity into event series at larger scales (Guhr et al., 1999;Gómez et al., 2002). To avoid this problem one can use a low order polynomial fit to the experimentally or numerically obtained cumulated event density N(T). The polynomial that minimizes square deviation from the staircase function of Eq. 2 is then used to perform the unfolding transformation Eq. 3, see e.g., Flores et al. (2000). ing situation can be seen as a smooth version of the abrupt one with equal steepness at both ends of the curve. Comparison with Figure 2A shows that similar variations of the cumulated event density indeed occur in LM data. In Figures 2D,G the measured IMI distributions P 1 (S) are shown. Here and in subsequent histograms the optimal bin width (Scott, 1979) is used, which balances bias and variance by minimizing the mean squared error between measured data and histogram bins. It scales like n −1/3 where n is the number of available data. The uncertainty of the number of events n(s) falling into a histogram bin centered at s is estimated by n s ( ) and symbolized as gray bars. Finally histograms (and errors) are downscaled in order to have unit sum ∑ = s n s ( ) . 1 In the abrupt situation ( Figure 2D) two distinct, narrow peaks are visible. For the linearly changing case ( Figure 2G) both have merged into a common, broader distribution. Different to the directly measured IMI distributions the unfolding procedure is able to reveal the intrinsic periodicity of the underlying system in both cases, see Figures 2E,H. Here clear peaks at s = 1 become visible after normalizing to the locally defined event density. Note that a simple global normalization would not succeed in either case. Ferri et al. (2006b) have found that log-normal distributions fit the pooled measured IMI distributions best. Here it was investigated how the unfolded IMI distributions of individual patients can be fitted by standard distributions. To be independent of any kind of binning the least-squares fit was performed for the cumulated carefully that the dependence of the results on the particular choice of m is very weak. This can for example be done by lowering/rising m of Eq. 4 by ∆m ∈ [−4, 4] (but respecting m ≥ 1), see Figures 3B-D.

QuantIfyIng the PerIodIcIty of the unfolded IMI dIstrIbutIons
In addition to the m-dependence of the nearest neighbor spacing distribution P 1 (s), which is extensively discussed in the present paper, also the next nearest neighbor distribution P 2 (s) and the number variance Σ 2 (l) (Guhr et al., 1998;Mehta, 2004) should be considered. It measures long-range correlations between events by quantifying the variance of the number of events in pieces of (unfolded) length l and is known to be much more sensitive than the spacing distributions. Figure 2B shows the unfolded cumulated LM density after application of Eq. 3. Systematic trends are absent here. The data fluctuate around the function f(t) = t and the average event distance is 〈s〉 = 1 throughout the whole data set (reflecting local normalization).

Application to model data
The performance of the unfolding is illustrated at the example of the pendulum with variable length l and period ∆T discussed already above: n = 400 events were sampled from a pendulum subject to strong environmental noise, which was modeled by a Gaussian distribution of width 7.5 s. In the first case the center of the distribution is abruptly changed from 15 to 60 s after 200 events, whereas in the second case the change is continuously modeled by linear interpolation between the limiting values. The cumulated event densities are different, see ences between the groups. To quantify the similarity of two IMI distributions P μ and P ν the KS test was employed again, now in its version for two samples (Siegel, 1956;Press et al., 1992). As before, it results in a distance d µν KS ≥ 0 and a significance 0 Here p µν KS depends on d µν KS as well as on the sizes n μ and n ν of the distributions. In the present paper p µν KS was used to quantify the similarity of two distributions with a single number in [0, 1] (0 for non-overlapping distributions and 1 for identical ones) in a size independent manner.
For automated similarity clustering on basis of the matrix with entries p µν KS the contrast turned out too weak. To deal with this problem the transformation: was performed, where 〈 〉 p µν KS and σ p denote the average and standard deviation of the off-diagonal p µν KS , respectively. This monotonic transformation sharpened the contrast by enlarging elements larger than average and decreasing elements smaller than average.
On the basis of these x μν the recently developed cluster participation vector (CPV) algorithm (Rummel et al., 2007Rummel, 2008) was employed. It operates on the eigenvectors of the symmetric similarity matrix with elements 0 ≤ x μν ≤ 1 (with ones on the diagonal) and performs two tasks in a fully data-driven manner. First, the number K of clusters present in the data can be estimated from features of all eigenvectors rather than pre-specifying it on the basis of previous knowledge. For this step another Kolmogorov-Smirnov test must be carried out, for which a significance level of α = 0.01 was chosen. Second, the data sets are attributed to the clusters on basis of the components of suitable linear combinations of the eigenvectors belonging to the K largest eigenvalues of the similarity matrix. The CPV algorithm has been tested extensively and was shown to outperform traditional approaches especially in situations with clusters that show sizable inter-cluster similarities.

results unfoldIng In sIngle PatIent data
Applying the unfolding procedure Eq. 3 revealed periodicity and allowed pooling of data. Figure 3A shows the example of a measured IMI distribution of a 75-year-old male with PLM and OSA (BMI = 27.3, AHI = 63.4, LMI = 77.1, PLMI = 70.7, PI = 0.62). The distribution of measured IMI is broad and almost flat between 5 and 60 s with only a small peak near IMI = 30 s. This distribution suggests a non-periodic system with all periods between 5 and 60 s appearing similarly often. Note that a small sample of IMI drawn with a random number generator (the prototype of a non-periodic system) from a uniform distribution between 5 and 60 s (and zero outside) could lead to a similar distribution. By definition such a LM series has PI = 1 although it is purely random (Ferri et al., 2006b).
Unfolding the data showed a clear contrast with a distinct peak at 1, revealing periodicity of LM that was previously hidden (Figures 3C,D). The intrinsically periodic character of the data was obscured in Figure 3A by the different slopes of the staircase function (see this patients' event times in Figure 2A), which represent dif-event densities. Like the unfolded IMI distributions P 1 (s) all used fit distributions had mean 〈s〉 = 1. All of them depended on a single parameter, which controls the width of the distributions in monotonous manner. Consequently, these parameters could be used to quantify PLM periodicity. If an unfolded IMI distribution P s a 1 ( ) ( ) was more narrow than a second distribution P s b 1 ( ) ( ), we concluded that the PLM of P s a 1 ( ) ( ) were intrinsically more periodic. For the log-normal distribution the fit parameter σ ln determines the width of the distribution directly, see Appendix for details and formulae. The log-normal distribution neither contains the exponentially decaying Poissonian case as a limit nor is it used in RMT for description of spacing distributions. Here it was used in order to tie in with previous work by Ferri et al. (2006b). Excluding one possible parametrization of the IMI distribution in favor of another might help to learn about the nature of the interaction (e.g., long-range vs. short-range) between subsequent LM in PLM series. Therefore we probed in addition two RMT spacing distributions that can be derived rigorously for classical one-dimensional particle gases on rings (Dyson gases), see Appendix for details and formulae. The Scharf-Izrailev distribution (Izrailev, 1988;Scharf and Izrailev, 1990) describes the spacings for Dyson gases with long-range interactions between all particles at temperature 1/β SI . Similarly, starting from Dyson gases with short-range interactions between neighboring particles at temperature 1/β sP one arrives at the semi-Poisson distribution (Bogomolny et al., 1999). For both distributions a decreasing temperature (increasing fit parameters β SI and β sP in the present context) translates into decreasing width. Both distributions interpolate between the exponentially decaying Poissonian (β SI = β sP = 0) and peaked distributions (β SI , β sP > 0).
To investigate which distribution fits the empirical distributions of unfolded IMI best, the non-parametric Kolmogorov-Smirnov (KS) test (Press et al., 1992) was employed. It gives a distance d KS ≥ 0 between the cumulated density of the sample distribution and the one of the functions given in explicit form in the Appendix (Eq. 7 for log-normal, the numerical integral of Eq. 8 for Scharf-Izrailev or Eq. 11 for semi-Poisson) and a significance 0 ≤ p KS ≤ 1. The latter depends on d KS as well as on the size n of the distribution and gives the probability that a distance d KS or larger is found although the null hypothesis of equal distributions is correct. We considered a fit as a good one if p KS ≥ 0.1 and discarded it otherwise.
In order to rank the performance of the fit distributions for clinical patient groups and automatically defined clusters (see Automated Similarity Clustering) we followed two approaches. The first one used the average p-value 〈p KS 〉 of a fit distribution over the individuals of the selected group. For the second one we gave credits to the fit distributions in the following way: zero credits if for a patient's unfolded IMI distribution p KS < 0.1, one credit if p KS ≥ 0.1, three credits for the best fit and two credits for the second best fit. In order to make the total number of credits comparable across groups or clusters it was divided by the number of group or cluster members.

autoMated sIMIlarIty clusterIng
Besides interpreting the distributions of measured and unfolded IMI for patients and clinically defined patient groups it was in addition investigated if it is possible to partition the total of the 30 data sets in a way that maximizes similarity within and differ-The shape of all three distributions was clearly different. In contrast, the unfolded IMI distributions in groups 2 (PLM without OAS or RLS) and 3 showed much more similar shape.
The majority of the unfolded IMI distributions of individual patients belonging to the three patient groups could be fitted well by the functions given in the Section "Methods." Patient group 3 (PLM with RLS) could be fitted best by the log-normal distribution according to both ranking criteria ( Table 2). This result confirmed a finding by Ferri et al. (2006b) who successfully fitted the pooled measured IMI distributions of RLS patients by log-normal distributions. For patient group 1 (PLM with OSA) the Scharf-Izrailev distribution turned out as a much better fit than the semi-Poisson or the log-normal distribution, implying that long-range interactions between LM could play a role in this patient group. For patient group 2 (PLM without OSA or RLS) there was disagreement between the maximum number of credits and the maximal average p-value 〈p KS 〉 thus implying fundamental differences in the shape of the unfolded IMI distributions within these patient groups. These differences were accounted for best by describing the data of different patients with different fit distributions. Overall, unfolded IMI distributions were fitted best by the Scharf-Izrailev distribution closely followed by the semi-Poisson distribution, see Table 2. With the exception of the RLS group the log-normal distribution led to the worst fits.
A Kruskal-Wallis H-test was carried out for group-wise different central tendency of the accepted (p KS ≥ 0.1) fit parameters σ ln , β SI , and β sP without finding significant differences (p > 0.1). Finally, it was found that the fit parameter σ ln of the N = 17 good fits of the log-normal distribution anti-correlated significantly (non-parametric Spearman correlation, p < 0.01) with the patients' ferent kinds of external influences on IMI duration such as, e.g., sleep stage or time elapsed after falling asleep. Both, the number variance ∑ 2 ( ) l ( Figure 3B) and the unfolded IMI distribution were rather insensitive to the particular choice of the polynomial degree of Eq. 4. Increasing m for every piece by ∆m = 4 led to very similar results. For both versions of the unfolded IMI distribution the p-value of a  in Appendix] was largest as compared to the remaining fit distributions and exceeds 0.1. Therefore the Scharf-Izrailev fit was accepted for further analysis, see below.
Pooled IMI from different patients may become broad and therefore suggest high randomness, i.e., non-periodic distribution of IMI. The measured IMI distribution of two patients showing PLM without OSA or RLS are given in Figures 3E,F (Figures 3G,H) the IMI distribution of both patients became similar and the estimated values of the lognormal fit parameters σ ln were very close. Pooling them would not alter the character of the distribution. The log-normal distribution (see Eq. 6 in Appendix) turned out as the best fit distribution in both cases and is indicated as fully drawn lines in Figures 3G,H. Although the log-normal distribution fitted the data better than the alternatives, both p-values remained below 0.1. Therefore these examples were excluded from our further analysis of the fit parameter σ ln .

unfoldIng In PatIent grouPs 1 to 3
The pooled measured IMI of patient group 1 (PLM with OSA) showed the most peaked distribution whereas patient group 3 (PLM with RLS) showed the broadest distribution, see To address question 1 the CPV clustering algorithm on basis of the similarities between the measured IMI distributions was employed (see Methods). Figure 5A shows the similarity matrix of the original pairwise Kolmogorov-Smirnov significances p µν KS . Data was ordered according to the clinically defined groups: PLM with OSA (a to j), PLM without OSA or RLS (k to t), PLM with RLS (u to δ). Apart from the diagonal and some exceptions in the group PLM with RLS, similarities were small. Application of Eq. 5 sharpens the contrast (Figure 5B). However, using these x μν for CPV clustering did not reveal a distinct cluster pattern. Determining the number of clusters in data-driven manner led to K = 10 clusters, the strongest one comprising the data sets o, w, y, and δ (data not shown). The second strongest cluster comprised data sets t, u, and x. Three clusters contained only two data sets, like, e.g., m and α. Figure 5B shows large elements of the similarity matrix x μν connecting these data sets. However, most importantly, none of the clusters was homogeneous consisting of patients belonging to the same patient group or sharing any other clinical or demographical characteristic. In summary, clustering on the basis of the similarities of measured IMI distributions did not seem to lead to clinically conclusive results.
Cluster participation vector clustering was repeated on the basis of the similarity of the unfolded IMI distributions. The elements p µν KS and x μν of the similarity matrices shown in Figures 5C,D were in general much larger than for the measured data. For better visual representation the data sets were rearranged according to K = 5 clusters that were found by the CPV algorithm. The clusters comprised data sets a to δ (cluster A, seven members), e to α (cluster B, six members), c to t (cluster C, eight members), y to x (cluster D, four members) and d and f (cluster E, two members). Data sets q, w, and b remained unclustered. Data set β (female, 21 years, RLS) had only slightly stronger similarities with cluster A than with B and therefore connected both clusters. Similarly, clusters A and E as well as C and D were strongly interconnected. Figure 6 shows the pooled unfolded IMI distributions for the clusters A to E of Figures 5C,D. Clusters B and E had the narrowest distribution notwithstanding a weak connection. Cluster E was more strongly connected with cluster A via data set d (male, 49 years, OSA), see Figure 5D. The strongly connected clusters C and D had similar distributions.
Analyzing medical conditions or medication that exacerbate or reduce PLM within these clusters did not reveal any common cluster specific features. Only the Kruskal-Wallis H-test for clusterwise different central tendency of the PI (Ferri et al., 2006b) was significant on significance level α = 0.01. For the cluster medians of the PI significant pairwise differences between clusters A and D as well as between clusters B and D (U-test) were found. Note that these cluster pairs were especially weakly interconnected in Figure 5D. In addition PI medians were marginally different for the following combinations: C vs. D (p = 0.016) and C vs. E (p = 0.044). The medians of the demographic data were only marginally significantly different. BMI and PLMI: cluster A vs. D (p = 0.042); AHI: cluster B vs. D (p = 0.038). It is interesting to note that Ferri's PI was the only quantity that showed significant differences between the automatically found similarity clusters. Mind that for the clinically defined groups no significant difference of PI was found (Figure 1). BMI, suggesting more overweight patients show a more narrow distribution, see Table 3. Investigation on a larger patient cohort is necessary to decide whether age, subjective sleepiness and PI (which turn out marginally significant, p < 0.05) are also candidates for significant anti-correlation. Very similar tendencies toward non-parametric correlation are present for the N = 19 accepted fit parameters β SI of the Scharf-Izrailev fit (now correlation instead of anti-correlation). In contrast no significant correlation with demographic data was found for the N = 18 accepted semi-Poisson fits.

data-drIven clusterIng
Visually comparing the measured and unfolded IMI distributions within clinically defined patient groups showed large dissimilarities occasionally, see Figure 3E and F for an example of measured data. On the other hand very similar IMI distributions were found for patients from different clinically defined groups. Therefore the following questions arose: 1. Is it possible to define new groups according to the shape of the IMI distributions using cluster algorithms? 2. Is there a common/are there shared clinical or demographic features that may cause the cluster formation?   For clusters B and D the log-normal distribution delivers clearly the best fit. However, this does not imply the fit parameter σ ln was the same in both clusters (see Figure 6): σ ln = 0.25 ± 0.03 for cluster B and σ ln = 0.55 ± 0.06 for cluster D (marginally significant difference of medians, U-test: p = 0.029). Note that for the clinically defined patient groups differences in σ ln were not significant (p > 0.1). Note also that the apparently similar unfolded IMI distributions of clusters B and E were best fitted by different Subsequently the analysis of the goodness-of-fit of the oneparameter fit distributions given in the Section "Methods" for the automatically defined similarity clusters was repeated, see Table 4. With the exception of cluster C there was agreement between both ranking methods for all clusters. In addition now the cluster-wise average p-value of the best fit distribution was larger, implying that clustering automatically groups patients whose IMI can be fitted best by the same distribution. x μν (B) of measured IMI distributions. Data sets were ordered according to patient groups and group boundaries are indicated by fully drawn lines: PLM with OSA (a to j), PLM without OSA or RLS (k to t), PLM with RLS (u to δ). The similarity in the group PLM with RLS was higher than for the remaining groups. Data-driven clustering predominantly depicted small clusters (see Methods). Bottom: Similarities p µν KS (C) and x μν (D) of the unfolded IMI distributions were much larger. Data-driven clustering found several large clusters. For illustration purposes, the channels were ordered according to five automatically detected clusters comprising data sets a to δ (cluster A, seven members), e to α (cluster B, six members), c to t (cluster C, eight members), y to x (cluster D, four members) as well as d and f (cluster E, two members). Cluster boundaries are again indicated by fully drawn lines.

Log-normal
their characteristics (Figures 3G,H). Technically, it is important to note that our choice of the degree m of the fit polynomial according to Eq. 4 prevents overfitting the data. Allowing for variations ∆m ∈ [−4, 4] we have shown in addition that the results are robust against variation of the particular choice of m (Figures 3B-D).
Here we investigated PLM periodicity by focusing on the shape of IMI distributions through the unfolding procedure. In patients with PLM and RLS Ferri et al. (2006b) found a bimodal distribution of measured IMI with the lower peak considerably below 10 s. This differs from our finding of a single-peaked distribution of measured IMI in patient group 3 ( Figure 4C). The reason for this difference might be that in contrast to our analysis these authors based their calculation of LM distribution on all IMI intervals >0.5 s in contrast to the WASM definition of PLM during sleep applied to our data (Zucconi et al., 2006).
PLM associated with OSA and RLS show distinct IMI characteristics (Culpepper et al., 1992;Briellmann et al., 1997). Therefore our second objective was to test whether periodicity could differentiate between the clinical groups containing patients with PLM and OSA (group 1), PLM and RLS (group 3) and PLM without RLS or OSA (group 2). PLM periodicity was assessed by fitting oneparameter distributions to the unfolded IMI distributions of individual patients. The pooled IMI distributions of the three groups appeared very similar (Figures 4D-F) and none of the patients' fit parameters σ ln , β SI , and β sP was significantly different between the groups. These results suggest that the degree of periodicity as quantified by these parameters is not related to the comorbidity (OSA or RLS). However, patients within groups were heterogeneous in terms of medication and comorbidities, thus representing a typical patient cohort in a sleep clinic. This heterogeneity contrasts with other studies that investigated patients with RLS and/or PLM excluding patients with medication influencing motor activity during sleep or significant sleep disorder or major comorbidities Garcia-Borreguero et al., 2004;Hornyak et al., 2004;Ferri et al., 2006aFerri et al., ,b, 2009Manconi et al., 2007).
Non-parametric correlation between the fit parameters and demographic patient data showed a significantly more narrow IMI distribution, i.e., a more pronounced periodicity, for overweight patients. BMI of patient group 1 (PLM with OSA) was significantly higher compared to groups 2 (PLM with RLS) and 3 (PLM without RLS or OSA), see Figure 1.
In contrast to the degree of periodicity (as defined by the values of the fit parameters σ ln , β SI , and β sP ), its nature (as defined by the best fit distribution) was different between patient groups. We fit distributions. This illustrates how automatic clustering reveals features of the data that cannot be obtained by visual inspection. For the unfolded IMI distributions of patients belonging to clusters A and E the Scharf-Izrailev distribution had the best fit performance. Here the fit parameters had different tendencies in both clusters (A: β SI = 5.21 ± 1.26, E: β SI = 9.76 ± 0.34), however, due to the small cluster sizes the difference did not reach significance (U-test: p = 0.095).

dIscussIon
In the present study, we applied for the first time the unfolding procedure of RMT to analyze PLM data. Our first objective was to test whether unfolding would enhance and reveal intrinsic periodicity that is not obvious from the measured IMI distribution. Indeed, we found that applying this procedure in single patients improved detection of periodicity, which was hidden within a broad distribution of measured IMI (Figures 3A,C,D). Thus, this result underlines the usefulness of unfolding for the detection of periodicity in neurophysiological time series. The broad range of measured IMI may be due to their dependence on the time of night and sleep stage (Coleman et al., 1980;Culpepper et al., 1992;Pollmächer and Schulz, 1993). In addition, the IMI probably may also depend on the time of intake of medication and its pharmacokinetics influencing motor activity in sleep.
Because unfolding may be necessary to enhance and reveal periodicity, we propose this method as an important pre-processing step before analyzing PLM. Note that this recommendation holds independently from the method used for recording the LM, be it EMG (as in our study) or other methods like the recently applied Emfit sensor (Rauhala et al., 2009). Pooling of measured IMI may also broaden a group's IMI distribution, whereas applying the unfolding procedure before pooling the patients' data retains   Table 2. The cluster-wise best fit is displayed in boldface.

Log-normal Scharf-Izrailev Semi-Poisson
Cluster long-range and short-range interactions, respectively. As a consequence, excluding one distribution in favor of the other one, could help to understand mechanisms underlying LM series.

Log-normal distribution
Similar to Ferri et al. (2006b) the log-normal distribution is used (see e.g., Limpert et al., 2001) with mean 〈s〉 = 1 and a fit parameter σ ln that controls the width:

Scharf-Izrailev distribution
A heuristic formula for interpolating between the exponentially decaying Poisson and the peaked Wigner distribution results from the spacing statistics of Dyson gases with long-range interactions at temperature 1/β SI (Izrailev, 1988;Scharf and Izrailev, 1990): Here increasing parameter β SI (decreasing "temperature") implies decreasing width of the distribution. The constants A SI and B SI must be calculated numerically for given β SI by demanding normalization ∫ = ∞ 0 1 1 ds P s ( ) and mean unit spacing 〈 〉 = ∫ = ∞ s ds sP s 0 1 1 ( ) . Also the cumulated density is integrated numerically. For integer values of β SI the Scharf-Izrailev distribution coincides with well known distributions of RMT. For β SI = 0 the Poisson distribution and for β SI = 1 the Wigner distribution are reproduced.

Semi-Poisson distribution
The semi-Poisson distribution describes the spacing statistics of Dyson gases with short-range interaction at temperature 1/β sP (Bogomolny et al., 1999) as well as incomplete spectra of complex systems (Hernández-Saldaña et al., 1999). In addition it has been found empirically that distributions of (unfolded) car distances on highways can be fitted by the semi-Poisson distribution (Krbalek et al., 2001): Also here increasing parameter β sP (decreasing "temperature") decreases the distribution's width. The normalization is given by