Effects of Ageing and Sex on Complexity in the Human Sleep EEG : A Comparison of Three Symbolic Dynamic Analysis Methods

Symbolic dynamic analysis (SDA)methods have been applied to biomedical signals and have been proven efficient in characterising differences in the electroencephalogram (EEG) in various conditions (e.g., epilepsy, Alzheimer’s, and Parkinson’s diseases). In this study, we investigated the use of SDA on EEGs recorded during sleep. Lempel-Ziv complexity (LZC), permutation entropy (PE), and permutation Lempel-Ziv complexity (PLZC), as well as power spectral analysis based on the fast Fourier transform (FFT), were applied to 8-h sleep EEG recordings in healthy men (n=31) and women (n=29), aged 20-74 years. The results of the SDA methods and FFT analysis were compared and the effects of age and sex were investigated. Surrogate data were used to determine whether the findings with SDA methods truly reflected changes in nonlinear dynamics of the EEG and not merely changes in the power spectrum. The surrogate data analysis showed that LZC merely reflected spectral changes in EEG activity, whereas PE and PLZC reflected genuine changes in the nonlinear dynamics of the EEG. All three SDA techniques distinguished the vigilance states (i.e., wakefulness, REM sleep, NREM sleep, and its sub-stages: stage 1, stage 2, and slow wave sleep). Complexity of the sleep EEG increased with ageing. Sex on the other hand did not affect the complexity values assessed with any of these three SDA methods, even though FFT detected sex differences. This study shows that SDA provides additional insights into the dynamics of sleep EEG and how it is affected by ageing.


Introduction
Neuronal interactions in the brain are highly nonlinear, making nonlinear time series analysis methods particularly appropriate for the characterisation of electroencephalogram (EEG) recordings [1].In fact, the introduction of nonlinear methods for the analysis of EEG signals [2,3] made the study of complex neural networks in the brain possible in ways that were not feasible with linear methods, such as the Fourier transform (FT).However, many of the classic methods used to this aim can provide spurious results when computed from noisy time series, such as EEG recordings [4].Therefore, research into nonlinear methods that are well-suited to the analysis of noisy complex time series, such as the EEG, is required.
One possible approach consists of using symbolic dynamic analysis (SDA) methods.SDA measures the complexity of a signal by quantifying the emergence of patterns within a symbol series of the original signal [5].In SDA a signal is represented using a series of symbols; this symbolic sequence is then characterised in ways that are dependent on the algorithm used.The symbolisation provides two main advantages over other techniques [4,6].First, the computational efficiency is better than that of other nonlinear 2 Complexity methods.Second, the robustness against noise is improved, even when applied to short time series [6].Furthermore, SDA requires no prior knowledge about how the signal of interest is generated.Once a good correspondence between the symbol series and the original time series is reached, SDA techniques provide information about the mechanisms underlying the generation of these signals [6,7].
In the current study, the SDA methods Lempel-Ziv complexity (LZC), permutation entropy (PE), and permutation Lempel-Ziv complexity (PLZC) were used.LZC was chosen because this technique was successful in characterising EEG signals in different conditions (e.g., epilepsy, Alzheimer's, and Parkinson's disease) [8][9][10][11][12][13][14][15].LZC was introduced by Lempel and Ziv [16] and characterises complexity in Kolmogorov's sense [17] (i.e., the smallest binary program capable of reproducing an information containing sequence).A description of the algorithm is given in [18], with signals being converted into finite symbol series using symbolic sequence decomposition.Different techniques can be used to determine a threshold value for the symbolisation of the original time series, including k-means value, mean value, median value, or mid-point [19,20].The median has proven to be robust to outliers within the signal and is commonly used in LZC applications [21][22][23].Following the partitioning of the original time series, symbol series are formed and scanned from the start to identify the appearance of new patterns or words in the symbolic sequence [18].The total number of patterns, normalised to take into account the sequence length, is a measure of complexity in the time series [18].In addition to being a complexity indicator, LZC also describes the dynamics of a signal (e.g., periodic signals result in lower LZC measures whereas complex signals result in higher LZC).
The second SDA technique used in the current study, PE, was selected because of its symbolisation algorithm, which provides information about the temporal order of the signal.PE is an entropy measure first introduced by Bandt and Pompe [24].It computes Shannon's entropy of the occurrence of patterns of a certain length in the time series.The symbolisation algorithm used in PE takes into account the order of the data points.Patterns are also influenced by the nature of the signal (e.g., a less regular signal contains a wider range of patterns compared to a regular one).PE values depend on input parameters such as the vector length (i.e., embedding dimension or embedded vector length ) and the lag (i.e., time delay ).PE has been applied to the analysis of various biological signals including cardiac [25,26], brain (EEG) signals in epilepsy [27][28][29] or under anaesthesia [30,31].Many of these studies focused on optimal parameter selection to characterise changes associated with a disorder.There also are numerous studies on how to select the values of the input parameters for PE computation (i.e.,  and ) [32][33][34][35][36].
The third SDA technique used in the current study is PLZC, a recently introduced method that brings together the symbolisation used in PE and the complexity algorithm used in LZC.PLZC was first proposed as Lempel-Ziv permutation complexity (LZPC) [37] but later standardised and renamed [38].Differences between the LZC, PE, and PZLC were assessed by using both synthetic and real EEG data [37].In these studies, the logistic map was used to characterise the methods' response to a well-known synthetic signal to create asymptotic bifurcations.It was found that LZC did not detect bifurcations in high chaotic regions but these bifurcations were detected with PE and LZPC.Even for low embedding dimensions (i.e., d = 3 which is denoted as  throughout this paper), PE and PLZC were able to detect the inherent oscillatory behaviour [37].The three methods resulted in low measures of complexity for nonchaotic regions and high for chaotic regimes.When these methods were applied to real EEG data comprising four epileptic seizures, LZC was not able to detect any changes in the data, whilst PE (m = 4,  = 1) and LPZC (m = 4,  = 1) detected seizures.In addition, PLZC enables detection of changes in a signal associated with the permutation process which reflects the relations of the signal points (i.e., the samples) [39,40].
In spite of the aforementioned evidence on the usefulness of nonlinear analysis techniques, including SDA methods, for EEG studies, it is necessary to validate results to ascertain if any changes observed arise from nonlinear properties of the signal and do not merely correspond to changes in the spectral content of the EEG [4,41].To that aim, numerous algorithms of surrogate data analysis have been introduced for the characterisation and modelling of nonlinearity in nonstationary and high dimensional time series [42][43][44][45][46].They allow investigating whether the time series of interest truly contains nonlinear structures or could equally well be described by linear models [42][43][44][45][46]. Thus, the results from these surrogate data analysis methods can determine whether the SDA results are not merely reflecting changes to the power spectrum of the signal but indeed reflect changes in the (nonlinear) dynamics of physiological signals.
EEG studies are particularly relevant in the context of sleep.Sleep is a fundamental physiological process that is subdivided in two states, nonrapid-eye-movement (NREM) and REM sleep which can be identified using polysomnography (PSG) which includes EEG recordings [47,48].Sleep and related EEG hallmarks show changes across the lifespan, as well as sex differences, which have been documented using power spectral density (PSD) measures [49][50][51][52][53].In sleep studies, the standard quantitative EEG analysis approaches have been based on FT algorithms [49][50][51][52][53].These studies quantified properties such as predominant brain activity in specific sleep stages (e.g., delta activity in slow wave sleep [54]), significant reduction in slow wave sleep with age and reduction in EEG power in the delta and theta frequencies particularly in NREM sleep [55], or sex related changes in the slow frequency range (e.g., elevated slow wave activity in women [56]).However, these changes in the EEG are not uniform and, thus, these recordings during sleep should be treated as nonlinear since transitions between vigilance states (VS) can also be random on top of the nonuniform physiological changes due to ageing and gender differences [48].
In the current study we applied and compared three SDA methods, i.e., LZC, PE, and PLZC, to characterise changes in brain activity during sleep based on PSG recordings.It was hypothesised that SDA methods would highlight changes in the nonlinear dynamics of the EEG during sleep and that they would be sensitive to ageing and sex differences.PSD and surrogate data analyses were used to validate whether the SDA methods reflected nonlinear dynamic changes in brain activity.1).Baseline recordings were obtained after the participants had spent at least one 'adaptation' night in the laboratory.PSG recordings were performed according to the International 10-20 system of electrode placement.We used EEG recordings from two central (C3 and C4), two occipital (O1 and O2), and two reference electrodes (A1 and A2) with C3-A2, C4-A1, O1-A2, and O2-A1 derivations.The EEG recordings started at lights-off at 23:00 until lights-on at 07:00 h.A band-pass filter (0.3-70 Hz) and a notch filter at 50 Hz were applied during the recording of the signals.In addition to these filters, at the time of signal processing, a band-pass filter of 0.5-40 Hz was also applied.The sampling rate was set to 256 Hz [56][57][58].The data were stored in European Data Format (.edf) and vigilance state (VS) classification was performed by trained sleep technicians using a standardised software application (Profusion PSG within Nexus Control, Compumedics Limited, Melbourne, Australia).Sleep staging was performed according to the Rechtschaffen and Kales [59] criteria on each 30-s epoch.NREM sleep stages 3 and 4 were combined to form slow wave sleep (SWS).Thus, the VS were denoted as wakefulness, NREM sleep stage 1, stage 2, SWS, and REM sleep.
(3)  and  (substrings of the symbol series) are allocated with the first and the second symbols of the symbol series, respectively, and complexity counter () is set to 1.
(6) c(n) is the complexity of the symbol sequence {()} which denotes the number of distinct words found in the sequence.The total number of sub-sequences present in {s()} had an upper bound [8] denoted as C(n) and would be quantified using the following equations.
2.3.Permutation Entropy.PE is an entropy measure which considers the order of the data points which reflects the temporal organisation of the time series [24].An example of the computation of PE is shown in Figure 1 where input parameters embedded vector length m = 3 and time delay  = 1 were used on an arbitrary signal.Embedded vectors are patterns within the signal where the order of the time points corresponds to the so-called motifs.Relative frequency of each motif is equal to the appearance of a motif within the time series.PE evaluation was based on the algorithm described in [26]: (2) Then, by using a time index defined as , elements of the embedded vectors (i.e., pattern) were rearranged into increasing order.Furthermore, these time indexes were used to create a symbol series with  the possible motifs that can be obtained from the embedded vectors.Additionally, to safeguard the computation of the PE in equal values of the time series, the same symbols were assigned to these, i.e., [+( 1 −1) = [+( 2 −1)], and then the patterns were put in a row vector () as follows.
(3) PE was calculated as in (6), where relative frequencies of the motifs (Figure 1) within symbol series () were computed using Shannon's entropy measure and represented as   .
(4) Last, permutation entropy was normalised as in (7) where a PE value independent of the embedded vector dimension was obtained.
(2) Complexity is then computed using steps (4)-( 6) of the LZC algorithm.The upper bound denoted as () is a combination of the complexity counter () and the number of distinct patterns found in the symbol series can be estimated as follows [8].
(3) PLZC output is then normalised using (9) where  denotes the total length of the symbol sequence.This step ensures calculated measures to be independent of time series and motif lengths.When  is very large, is used to compute the complexity measure.
2.5.Surrogate Data Analysis.Surrogate data (i.e., surrogate time series) were created according to the algorithm introduced by Palus and Hoyer [46], which can be summarised as follows: (1) Fourier coefficients of the signal were obtained by fast Fourier transform.(2) Phases of these coefficients were randomised but the amplitude coefficients were not altered to ensure spectral consistency.(3) Next the inverse Fourier transform was used to create a surrogate time series with the same spectral characteristics as the original one, but different dynamics.
LZC, PE, and PLZC were then applied to these surrogate time series (STS) and the real time series (RTS).Then, we tested whether the complexity measures obtained both from the RTS and STS were significantly different from each other.

Results and Discussion
3.1.Parameter Selection for SDA Techniques.LZC, PE, and PLZC were computed on each sleep epoch (30-s) for each VS.These complexity measures were then averaged over thirds of the sleep period and grouped according to age and sex.For LZC, the median was used as threshold in the symbolisation of the time series due to its robustness to outliers.In the case of PE, there are no definite guidelines on how to select the input parameters (vector length and time delay).The time delay  was chosen as 1 to ensure no downsampling of the time series occurring in the symbolisation.Regarding m, the use of larger vector lengths improves regularity estimation [35,36].
On the other hand, the length of the time series must be larger than m! [23].PE was evaluated with different combinations of input parameters.Table 2 summarises the percentage of epochs showing significant differences (Related Samples Wilcoxon Ranked Sign Test) between PE values for RTS and STS.Based on these results, m = 6 and  = 1 were chosen to compute PE in this study.In this way, the percentages of epochs showing significant differences between RTS and STS is maximised and, at the same time, the use of a large vector length would improve regularity estimation with PE.
As PLZC and PE use the same symbolisation algorithm, the input parameter values chosen for PLZC analysis were also m = 6 and  = 1 for consistency with PE.

Effects of Vigilance States, Ageing, and Sex on Complexity
Measures.Complexity values for the C3-A2 derivation were used in the figures illustrated and described in this section as a preliminary analysis showed no significant interaction for factors 'electrode position' x 'third of the night' ( (6,240) = 1.39, p = 0.218).
In Figure 2, LZC, PE, and PLZC (left panels) values in different VS were illustrated together with PSD (right panels) measures of the same VS.The three SDA techniques were able to identify significant differences between VS across the whole night sleep (p values reported in Table 3).LZC, PE, and PLZC values in wakefulness (a) and REM (i) were higher compared to NREM sleep stages ((c), (e), and (g)).A gradual decrease in complexity from stage 1 to 2 and to SWS was also found and illustrated in Figures 2(c), 2(e), and 2(g) and also in Figure 3 across thirds of sleep.Additionally, EEG power spectra for all VS (Figure 2, right panels) and similar trends in PSD across VS were observed in all age and sex groups.PSD amplitude was lower in older males compared to all other groups.In older males, complexity was high in all VS (Figure 2, left panels).This inverse result between PSD and complexity measures was also illustrated in Figure 3 across thirds of the night sleep.SWA decreases during the course of the night; the decrease in LZC is consistent with this finding.By contrast, the increased complexity found with PE and PLZC likely reflects the increased time spent in REM sleep states and reduced SWS (i.e., increased activity in the brain) during the course of the night.
Power was measured in the low delta (0.5-4.5 Hz) range to characterise changes in SWA (Figures 3(j)-3(l)).As NREM sleep stage deepens, i.e., when sleep progresses from stage 1 to 2 and to SWS, SWA increases, in accordance with previous reports [53][54][55][56][57]. Also in accordance with previous reports, SWA decreases from the beginning to the end of sleep.In NREM sleep (Figure 3), LZC complexity values in first, second, and third thirds of the night were markedly lower in the young than in other age groups (p = 0.010 and p = 0.001; p = 0.022 and p = 0.002; p = 0.023 and p = 0.003 middle and old age groups, respectively).PE yielded significantly lower complexity values in the young age group compared to middle and older age (p = 0.002, p = 0.045).In the first third of sleep, PLZC also identified significant differences (p = 0.044) between young and middle age groups (Figure 3).Lower LZC, PE, and PLZC values were obtained when an increase in SWA was found.This is due to the fact that SWA reflects in the EEG in which the rate of change in the dynamics is not as rapid as in activated states, which lowers the complexity.This creates an inverse relation between the PSD and the complexity measures.However, particularly during SWS, this is not observed in older males during the last third of the night which might be due to the reduced number of SWS episodes or increased number of transition episodes between SWS and REM sleep, both of which need further investigation.
Furthermore, LZC values in females were higher than males in all VS (see Figure 2, left panels, and     but significantly different in sleep states (i.e., NREM and REM sleep) indicating higher variations in the EEG, high amplitude brain activity.This variation in the EEG is not consistent with previous findings showing increased SWA in women using FT [50][51][52] which would be expected to be highlighted as lower complexity values.By contrast to LZC analysis, sex differences were only found in REM sleep for PE and PLZC analyses.The differences in PE and PLZC values in females in NREM sleep stages were not significant (Figure 3).In SWS, on the other hand, complexity computed with all three methods showed the lowest values suggesting that ageing has a greater effect on the complexity of brain activity than gender.
Surrogate data analysis of the current data set showed that the complexity estimated by LZC merely reflected spectral changes in the signal.LZC captured the changes caused by ageing and sex differences similar to traditional spectral sleep EEG analysis, which may be primarily a consequence of changes in the amplitude of EEG activity (i.e., SWA).However, comparison of PE and PLZC results obtained from the real and surrogate data analyses showed that these methods highlight the nonlinear dynamic changes in brain (EEG) activity.Moreover, the LZC results were consistent with previous sleep studies based on the Fourier analysis [50][51][52].Our findings therefore corroborate that LZC measures changes in the brain activity which are prone to be reflecting changes in the EEG power spectra.This includes the characterisation of sex differences.
Overall, all three methods captured differences in brain activity in all VS, although differences between the methods were identified (e.g., LZC showed decrease in complexity whereas PE and PLZC showed increase in the course of the night, Figure 3).

Comparison between LZC, PE, and PLZC.
In the current study, PE and PLZC yielded similar results.Both methods could distinguish between brain activities in different VS and detect physiological changes (e.g., ageing) that affect the structural mechanisms underlying the brain activity, in agreement with previous studies [52].These results show the potential usefulness of these methods in sleep research, where nonlinear dynamic changes in the brain activity could be characterised by relative increases and decreases in complexity.Furthermore, we have corroborated our findings with surrogate data analysis.In contrast to LZC, PE and PLZC reflect changes in the brain dynamics and not mere changes in the spectral content.
The simplicity of SDA methods, their robustness to noise, and their computational efficiency have been reported in the literature [5].Nevertheless, the main issue with SDA techniques is the discretisation of the signal: when a limited number of symbols are drawn from the finite alphabet, the symbolised time series may not capture the dynamics of the real time signal accurately [61].This is an issue in particular with LZC, which makes the method amplitude dependent due to the use of limited number of symbols (i.e., mean and median technique, symbols "0", "1") [62,63].Thus, PE is a superior SDA technique compared to LZC as it considers the temporal organisation of the samples within the signal.This is also an advantage of PLZC as the symbolisation algorithms of these methods are similar.All of the methods are robust against noise and no knowledge about the signal dynamics is required [5,6].However, PE and PLZC require two input parameters for their computation and the values chosen have an impact on the PE and PLZC values.Even though there is no consensus on the optimal input parameters, it has been suggested that the embedded vector length m must be chosen taking into account the length of the signal being analysed [24,[33][34][35][36]. Therefore, we selected a large value of m, to improve regularity estimation [35] whilst at the same time meeting the epoch length requirements [24].The effects of time delay have also been discussed in the literature and it has been suggested to test different combinations starting from  = 1.However, increasing the time delay is equivalent to downsampling the time series, which could lead to a loss of information by not following Nyquist criterion.Furthermore, the combination of embedded vector length and time delay also affects the maximum frequency range coverage (10) that can be achieved.The maximum frequency that can be measured in relation to the sampling frequency f s and input parameters m and  is given by [64]   =   ( * ) .
In this study, EEGs were sampled at 256 Hz.Therefore, with the chosen combination of m = 6 and  = 1 frequencies up to 42.67 Hz would be analysed.This is an important point to consider when deciding the input parameter values, especially if changes happening at different frequencies need to be captured when computing PE and PLZC.
Another difference between LZC, PE, and PLZC is the computational cost of these techniques.In our study, MATLAB5 2013a was used for all mathematical operations including signal extraction, nonlinear signal processing, and epoch-by-epoch statistical analysis between RTS and STS.As a by-product of the study, the computational performances between these SDA techniques could be compared (device features; Vostro Desktop 470MT, i7 microprocessor, memory 8GB).LZC was the fastest analysis method (e.g., 1 sleep period, 4 electrodes took ∼20 min to compute).PLZC (under the same conditions it took ∼2h to compute) was the second most efficient method, whilst PE (∼6h to compute under the same conditions) was the slowest of all three SDA methods.Taking into account the symbolisation and complexity computation steps of each method, LZC and PLZC are faster compared to PE due to the complexity evaluation step.Given that PLZC revealed similar characteristics of the signal to PE and that LZC differences were merely highlighting changes in the spectral content of the EEG, it could be argued that PLZC is the most efficient method of the three to characterise nonlinear changes in sleep EEG.
Our study has some limitations that should be mentioned.Firstly, the sample size for some age/sex subgroups was relatively small with different number of males and females (e.g., 5 females and 14 males in the young group; 12 females and 6 males in the middle age group).As a result, our findings, especially those highlighting differences related to gender, are preliminary and our study should be extended on a larger sample size with balanced groups.However, the total number of n=29 for women and n=31 for men is sufficient in terms of statistical power when the age groups are aggregated.Secondly, the selection of input parameters for the computation of PLZC and PE was based on an empirical approach that depended on the characteristics of the particular database used in this study.Therefore, further work with different EEG databases must be carried out to check the validity of our choice of input parameter values in different scenarios.Last, but not least, there are different implementations of PE.Thus, it would be worthwhile evaluating whether other PE algorithms show similar changes in EEG dynamics to those reported in this study.
In summary, LZC, PE, and PLZC proved to be useful in the analysis of EEG data for the characterisation of different aspects of sleep.LZC was found to reproduce the findings of traditional Fourier-based sleep research.On the other hand, PE and PLZC revealed nonlinear changes in the dynamics of the EEG and showed significant changes in complexity with regard to VS and ageing.

Conclusions
In this study, three nonlinear SDA techniques (LZC, PE, and PLZC) were used to characterise brain activity from EEG recordings during sleep, as well as to investigate effects of ageing and sex differences.Our findings suggest that LZC changes merely reflect changes in the PSD of EEG signals, but that PE and PLZC provide additional information about the nonlinear dynamics of the signal that cannot be captured with conventional linear analysis methods.This work could be extended with the application of SDA to characterise changes in the nonlinear dynamics of the EEG in different datasets (pharmacological manipulations of sleep, sleep in patients suffering from dementia, etc.).This could potentially reveal information that is contained in the nonlinear dynamics of the EEG providing insights into brain functioning under different conditions in ways that are not possible with linear methods.

Figure 1 :
Figure 1: Illustration of the calculation of permutation entropy (PE) for embedded dimension m = 3, and time delay  = 1.Possible motifs m! = 6 for m = 3 are given at the top panel (a) of the figure with translations of those from the arbitrary signal on the left panel (b).At the lower panel (c), a histogram of the occurrence of patterns can be seen which are used in PE computations.

Table 1 :
Number of men and women in three age groups from whom EEG recordings were obtained.

Table 2 :
Percentage of epochs showing significant differences (Related Samples Wilcoxon Ranked Sign Test) between RTS and STS when PE was computed with different input parameter values.Hz.PSD were averaged across to obtain an average PSD for each 30-s epoch.PSD measures of the whole EEG recordings were classified into VS in each age and sex group.Slow wave activity (SWA) was defined as power spectral density in the 0.5-4.5 Hz range.2.7.Statistics.Statistical analyses were performed using the IBM SPSS Statistics v.23 software.Probabilities of p < 0.05 were considered as significant.LZC, PE, and PLZC values were computed from both RTS and STS to detect dynamic changes in brain activity.A nonparametric Wilcoxon Signed Rank Test for related samples was performed on the LZC, PE, and PLZC values obtained from both RTS and STS.Effects of VS, thirds of sleep, age, and sex were determined using a mixed model repeated measures analysis of variance, ANOVA.Subsequently, post hoc Bonferroni and Dunnett's tests were applied to characterise main and interaction effects.

Table 3 :
Effects of thirds of the night, vigilance state, age, and sex on complexity measures assessed by repeated measures ANOVA.