Bispectral Analysis of Parkinsonian Rest Tremor: New Characterization and Classification Insights Pre-/Post-DBS and Medication Treatment

Rest tremor is a most common symptom of Parkinson’s Disease (PD), with diagnosis and severity estimation often being hindered by subjectivity and limitations of existing methods. Additionally, treatment strategies for PD face challenges in monitoring symptom fluctuations with clinical methods, like the MDS-UPDRS and motor diaries, suffering from subjectivity, limited sensitivity, and infrequent sampling, potentially impacting treatment effectiveness. Hence, methods that can yield explainable biomarkers that accurately describe properties of PD rest tremor, while accounting for the presence of ongoing treatments, such as Deep Brain Stimulation (DBS) and medication, are important for integration in accurate AI-powered wearable systems. To that end, a Higher Order Spectrum (HOS)-based analysis to extract features from index finger velocity recordings of 16 PD patients is proposed. Two different scenarios are implemented for characterizing and classifying treatment (medication/DBS) effectiveness and tremor severity (Low-/ High-Amplitude (LAT/HAT)), by means of statistical tests and a leave-one-subject-out cross-validation classification scheme, respectively. The proposed analysis resulted in area under the Receiver Operating Characteristics curve (AUC) score of 0.94 for the Medication treatment classification, 0.71 for the DBS treatment and 0.83 for LAT/HAT prediction. Our results demonstrate that the proposed approach can effectively assess the influence of medication and DBS and recognize rest tremor severity. Finally, our HOS-based methodology enables the establishment of new rest tremor classes, based on its nonlinearity and allows for new insights about the dynamic nature of the resting tremor production system.


I. INTRODUCTION
Parkinson's Disease (PD) is a common progressive neurodegenerative disorder that causes death of dopaminergic neurons in the substantia nigra [1].PD is mainly characterized by motor symptoms, such as rest tremor, bradykinesia, The associate editor coordinating the review of this manuscript and approving it for publication was Rongbo Zhu .muscle rigidity and postular instability, as well as non-motor symptoms, such as cognitive impairments (e.g., dementia), psychiatric disturbances and sleep dysfunction [1].Individuals with PD, typically develop mild symptoms at the beginning of the disease onset [2].One of the most common symptoms of PD is rest tremor and occurs when the muscles are relaxed and still [3].Another, less common type of PD tremor is action tremor, which occurs with the voluntary movement of a muscle [4].There are several sub-classifications of action tremor, such as postural, kinetic and isometric [5].
For alleviation of PD symptoms (e.g., rest tremor), medication treatments, such as L-Dopa (levodopa), are commonly administered to patients [6].In addition, individuals can take advantage of advanced treatment methods such as Deep Brain Stimulation (DBS) [7], to counterbalance the symptoms that they experience when the effect of a medication dose disappears [1].DBS involves surgical placement of unilateral or bilateral leads (wires) transcranially in the subthalamic nucleus (STN) or the globus pallidus interna (GPi) [1], [7].
The treatment strategy design varies and depends on symptom timelines and presentations.However, there are challenges in monitoring symptom fluctuations over time.On clinical settings, PD rest tremor is often classified based on amplitude and frequency characteristics [8], as these are depicted through the commonly used Movement Disorder Society-Unified Parkinson's Disease Rating Scale (MDS-UPDRS) [9], without the use of instruments.However, such screening methods are often limited by subjectivity of the expert or the subject's answers [10].Currently, motor diaries are the standard for tracking changes outside the clinic, but they have compliance issues, subjective errors, and limited sensitivity [11].Most importantly, such assessments are a mere snapshot of the patient's disease staging, in fact sampled at a very low frequency [12].As a result, the treatment plan might not be effective enough, resulting in insufficient symptom relief, or it could be overly aggressive, leading to side effects.
Thus, technological interventions that provide more flexible and accurate ways of evaluating rest tremor severity, are of great importance in devising more efficient treatment strategies and early diagnosis and prediction methods.Wearable ubiquitous sensing alleviates shortcomings of traditional clinical evaluations as it fosters continuous and fine-grained unobtrusive monitoring, hence enabling evaluation of the patient's real-time condition in the full spectrum of human activity [13].Consequently, technological advancements in the field of wireless sensing have brought forth a booming interest in artificial intelligence (AI) algorithms that can yield insights based on the recorded data [14].Although these approaches maintain an efficient outlook in successfully and ubiquitously monitoring PD patients, they are based on AI methods, that lack explainability and a deeper understanding of the inherent nature of the disease [14].Reliable disease biomarkers that can train these AI methods more efficiently can solve one of the most critical open problems in PD monitoring, that lays in characterization of PD across disease staging, treatment conditions and different daily life and activity contexts [14].
In this vein, our research aims to develop explainable biomarkers that can effectively assess the nonlinear properties of the tremor producing system, while being sensitive to alternating treatment states.For that cause, our work introduces the use of Higher-Order Statistics/Spectra (HOS), a signal processing framework that is related to and can be expressed in terms of the moments of a random process [15].One of the most commonly used members of the Higher-Order Spectra family is the Third-Order Spectrum, known as the Bispectrum (BS), which is defined as the two-dimensional Fourier Transform (2D-FT) of the third-order statistics [15].HOS have found applications in many scientific fields, especially in biological systems, where nonlinearity is an inherent component of biological mechanisms.To our knowledge, no other studies have used HOS-related features for characterization and identification of PD rest tremor.In this way, our work moves one step further than the conventional approaches of amplitude and frequency representations, delivering a new perspective of PD rest tremor in terms of nonlinearity.
The analysis of our results indicates that our proposed methodology holds potential in describing various important aspects of PD rest tremor.The main contributions of our work can be summarized as follows: • The proposed methodology achieves high performance, while maintaining a robust behavior across different conditions, in recognizing the effect of medication and DBS and classifying subjects to the HAT/LAT classes.
• In agreement with previous studies [8], [16], [17], we also show that medication performs better than DBS in suppressing tremor and that a tremor rebound occurs 15 minutes after DBS is switched off [8].We also observe a secondary rebound that was not reported in [8], 60 minutes after DBS is switched off.
• We show that system-obtained information from HOS, such as nonlinearity, provide a robust performance, regardless of tremor severity.
• A precise description of the dynamical nature of PD rest tremor using HOS is introduced, effectively modeling the changes in tremor properties, caused by the transitions into and out of the various DBS and medication conditions.
• Finally, we propose a more efficient characterization of PD resting tremor based on nonlinearity, rather than amplitude-and frequency-based representations, through a re-clustering of the HAT/LAT classes based on HOS features.
The structure of the paper is organized as follows.In Section II, related work on treatment effect and PD tremor characterization and identification are presented.The theoretical background and the proposed methodology are presented in Section III.In Section IV, the employed dataset is described along with the different analysis scenarios used in performance evaluation.In Section V, the results of the proposed analysis are described.In Section VI, the results are discussed, along with implications and extensions of our results, including a HOS-based clustering analysis to further unveil a new perspective about PD rest tremor characterization; Section VI finalizes with the limitations of our work.Finally, Section VII concludes the paper.

II. RELATED WORK A. EFFECT OF DEEP BRAIN STIMULATION AND MEDICATION ON PARKINSONIAN TREMOR
The significance of PD tremor-related biomarkers that can be easily extracted using readily available commercial devices, such as inertial measurement units (IMUs), has recently been emphasized.These biomarkers have the potential to facilitate the design of advanced and personalized treatment strategies and the study of treatment effects on PD tremor management.For instance, in [18], a closed-loop DBS control system, capable of adjusting DBS settings by analyzing simple tremor signals recorded from the patient's index finger was shown to outperform conventional non-adjustable systems.Furthermore, in an endeavor to enhance levodopa treatment for advanced PD, Teymourian et al. explored a novel autonomous wearable system to medication delivery and wearable symptom tracking [19].Such a system could enable real-time monitoring of levodopa levels, allowing for personalized dosing adjustments in response to symptom fluctuations.
To enhance DBS settings optimization, Gülke et al. introduced a semi-automatic system, worn on patients' fingers, combining IMU data and clinician assessments [20].Time and frequency-domain features were employed to assess tremor severity.The results from ten PD patients demonstrated significant symptom improvement, promising advancements in closed-loop DBS optimization systems.In another study by Heldman et al., a novel DBS programming optimization method was introduced, utilizing a motion device on the more affected hand's index finger [21].Similar to the previous approach, this device integrated IMU sensors and conducted time and frequency-domain analyses.The software systematically adjusted stimulation settings to successfully minimize tremor and bradykinesia severity, side effects, and battery usage.Moreover, Castaño et al. designed a fine-motor-skills task to extract time-domain features sensitive to DBS-induced changes in motor function when switching between ''On'' and ''Off'' states [22].The reported results indicated sensitivity to motor function changes caused by DBS alterations.
Smits et al. investigated the use of graphical tasks to assess changes in upper limb function in response to dopaminergic medication in PD patients [23].Tremor fluctuations were quantified using power spectral density features of IMU data.Results demonstrated the validity of the proposed approach and provided insights into medication response and tremor reduction.In a study by Pulliam et al. [12], the challenge of monitoring temporal patterns of symptoms in PD patients in response to levodopa was addressed through wearable IMUs placed on each wrist and each ankle.Spectral power in frequencies typical of tremor, along with multiple regression models, were used to quantify the dose response of tremor, bradykinesia, and dyskinesia during various activities.Results from 13 PD patients demonstrated the potential of the proposed system for monitoring motor fluctuations in daily life.In another study by Ricci et al. [24], the objective assessment of levodopa treatment in PD patients was enhanced by devising multiple motor tasks, each related to different PD symptoms.IMUs were used to measure motor tasks in 36 patients at therapy initiation, six months, and twelve months.The power and amplitude of tremor within the PD tremor frequency band were used to assess treatment effects, yielding insights about the levodopa therapy.
Furthermore, certain studies have explored the combined effects of DBS and medication.In a study by Beuter et al. [8], the authors investigated the impact of DBS and medication on PD rest tremor by analyzing IMU signals.While DBS did not significantly affect patients under medication, it had a notable impact on both the amplitude and frequency characteristics of rest tremor when medication was not administered.Additionally, Meka et al. recruited 27 PD patients and evaluated their tremor levels across four different states, which were combinations of the presence and absence of DBS and medication [25].IMU recordings from the patients' middle finger were used for tremor quantification, and the area under the 3 to 7 Hz band was calculated and compared within these four states.Results demonstrated that while both DBS and medication were effective when used individually, the most significant tremor reduction occurred when both DBS and medication were combined.

B. PARKINSONIAN TREMOR CHARACTERIZATION AND IDENTIFICATION FROM IMU DATA
Characterization and identification of PD tremor constitute established research topics and hence, several studies have explored PD tremor dynamics in order to create workflows that can accurately characterize PD tremor instances.In this vein, di Biase et al. effectively discriminated PD from essential tremor in IMU data with robustness in transitions across postures, by designing a ''tremor stability index'' through a logistic regression analysis that encompassed frequency and amplitude characteristics [26].Furthermore, Hssayeni et al. used data from wrist-and ankle-mounted IMUs from PD patients during daily life tasks in lab settings [27].Medication effect was incorporated in the data collection protocol, which overall aimed to assess the severity of PD tremor.By utilizing the gradient tree boosting and long-short term memory (LSTM) models, they successfully predicted the expert-assigned UPDRS-III severity annotations [27].
The advent of wearable devices and the development of AI, enabled the integration of wearable technology with algorithms for on-the-edge computing.Pioneering approaches emerged such as the one proposed by Patel et al., who evaluated the severity of tremor and other motor symptoms by means of IMUs mounted on the body [28].They extracted amplitude and frequency features and utilized support vector machines to distinguish between different motor symptoms and quantify their severity [28].Furthermore, Chén et al.
addressed identification and severity assessment of PD through a multi-modal smartphone behavioral data approach that included the modality of tremor [29].Through a custom-made ML model, they analyzed week-long data obtained over a 6-month interval and identified disease specific feature-profiles [29].
Other works have used deep learning models in an attempt to address PD identification.Oktay and Kocer leveraged an LSTM model to differentiate between PD and essential tremor from resting and postural IMU data [30].Similarly, Sun et al. designed TremorSense, a framework for discerning rest, postural and action PD tremors, through a convolutional neural network (CNN) [31].Papadopoulos et al. employed expert-derived UPDRS tremor annotations from IMU data of PD patients and healthy controls collected in-the-wild during phonecalls, to design a hybrid multiple instance learning -CNN model [32].Later, the same group designed a novel deep learning architecture to capitalize the information of IMUs and touchscreen typing data towards identification of PD symptoms including tremor in a PD and healthy control cohort, collected during users' regular interaction with their smartphones [13].Finally, in Lamprou et al., a previous work of ours, we have utilized the same in-thewild phonecall data as Papadopoulos et al. [32] to introduce DeepBispecI, a comprehensive tremor detection framework that outperformed other methods by means of HOS images that were fed to a CNN [33].
Ultimately, our research is based on the PD rest tremor dataset introduced by Beuter et al., that focuses in exploring the effects of DBS and L-dopa medication through IMU data from a velocity-transducing laser, aimed at the index finger of each subject [8].Works that utilized this dataset so far used deep learning to predict different attack stimulations in DBS [34], employed discrete non-Markovian stochastic processes to study the Parkinsonian pathological tremor [16], applied nonlinear analysis to characterize tremor properties as chaotic or stochastic [35], or considered Long-Term Correlations (LTC) and Multifractal (MF) properties to tackle three different classification problems [17].The Fast Fourier Transform (FFT) [36] and the Power Spectrum (PS) [8] have also been utilized to analyze tremor amplitude and frequency characteristics of resting tremor signals.
To summarize, although the related work in the effects of DBS and dopaminergic treatment, as well as in PD tremor characterization from IMU data has progressed significantly, the following gaps can be identified: (i) to the best of our knowledge, biomarkers that describe tremor properties and are correlated to disease severity, while accounting for simultaneous presence of DBS and medication, do not exist in the literature; (ii) amplitude, frequency and harmonicity expressions of tremor hold clinical value, yet recent research has proven that these properties are not sufficient for describing the nature of PD tremor.More advanced approaches such as the one by Sarbaz and Pourakbari [35], contributed significant characterization insights but did not yield biomarkers that can track dynamic fluctuations related to treatment regimes.
Our proposed methodology is designed to address the above significant gaps that exist in PD tremor identification and characterization, by designing explainable nonlinear biomarkers that are able to describe the tremor producing system and track dynamic alterations that accompany treatment regimes transitions.

III. METHODOLOGY A. THEORETICAL BACKGROUND
HOS have some appealing properties and are appropriate for the examination of complex nonlinear systems [15].HOS do not only reveal amplitude information about a random process, but also reveal its phase information.This is of great importance, since, as is well known, second-order statistics, are phase-blind.By taking the multi-dimensional Fourier Transform (FT) of cumulants, the space of higher-order spectra is defined [15].The most commonly used member of the higher-order spectra family is BS [15].When applied to a random process that either is non-Gaussian or originates from a nonlinear system, the BS is a powerful tool, because unlike the PS, it can detect non-linear interactions, such as quadratic phase coupling (QPC) [15].
The BS(f 1 , f 2 ) of a random process {x(k)} can be defined as: where E[•] is the expectation value, X (f ) is the FT of the random process {x(t)} at frequency f and X * (f ) is its complex conjugate [15].The normalized Bispectrum or Bicoherence is defined as: where BS(f 1 , f 2 ) is the bispectrum of the random process {x(k)} at frequencies (f 1 , f 2 ) and P(f ) is its power spectrum at frequency f [15].The magnitude of the bicoherence, |, is called the Bicoherence Index (BI) and takes values between 0 and 1. BI quantifies the presence of QPC between any two frequency components in a random process due to their non-linear interactions.In particular, two periodic components at frequency and phase pairs of (f 1 , φ 1 ) and (f 2 , φ 2 ), respectively, exhibit QPC when there exists a third periodic component at frequency and phase pair of (f 3 , φ 3 ), with its frequency and phase satisfying close to 1 translates into an almost perfect phase coupling, whereas close to 0 denotes the absence of QPC [15].

B. PROPOSED METHODOLOGY
The proposed methodology, presented in Fig. 1, aims to extract HOS-based features that are able to accurately discriminate between the HAT/LAT classes and capture the effects of DBS and L-dopa medication.Furthermore, a classification scheme is designed to test the efficiency of the extracted features.Following, the feature extraction procedure and the adopted classification scheme, are presented.• BS-based Features 1) Bispectral Entropy(BspecEnt1): where and L denotes the BS region of interest and m,n = 1,2,. . .,M , where M is the number of points in L. 2) Bispectral Squared-Entropy (BspecEnt2) at [2.5-8]Hz: q mn ln(q mn ), where and L denotes the BS region of interest and m,n = 1,2,. . .,M , where M is the number of points in L. 3) Area of Dominant Peak (MaxArea): A measure of the area around the dominant peak of the BS domain.Area is measured in units of bi-frequency bins, where one bi-frequency bin denotes the unitary area.4) Total Area (TotalArea): A measure of the total area of the BS domain.
ESR ranges from 0 to 1, with the value of 1 indicating a perfect elliptical shape of the BS content.The motivation behind the area-and shape-related features is that HAT subjects tend to have a sharper and more concentrated BS dominant peak (Fig. 2-right panel).Thus, we expect that, as the tremor enters the pathogenic region, the main peak will converge to a Dirac delta function, maintaining an elliptic shape as it shrinks.On the other hand, a more dispersed BS content, indicates that tremor 'escapes' from a pathogenic behavior and tends to become physiological.In this way, more frequency components are present in the signal, whereas the shape of the contour becomes more and more irregular (Fig. 2-left panel).These claims also remain valid for the NoPeaks, BspecEnt1 and BspecEnt2 features.In fact, a high entropy value indicates a BS distribution that approximates a uniform distribution, where randomness is maximized, whereas a low entropy value indicates that the BS distribution has a very concentrated form and approximates a Dirac delta function where randomness is zero.
• BIC-based Features 1) BIC Entropy (BicEnt1): where and L denotes the BIC region of interest and m,n = 1,2,. . .,M , where M is the number of points in L. 2) BIC Squared-Entropy (BicEnt2): q mn ln(q mn ), (10) where and L denotes the BIC region of interest and m,n = 1,2,. . .,M , where M is the number of points in L. 3) Total BIC (TotalBic): and L denotes the BIC region of interest.
In the examined problem, the selection of L = [2.5 − 8]Hz region is mainly considered, as HAT/LAT tremor activity that involves QPC interactions between frequency pairs, is located at this area (Fig. 2-third row).Moreover, in the case of HAT, the BS peak is even limited within the [4][5][6]Hz region and has an elliptical shape (Fig. 2-third row, right panel).On the contrary the BS content of the LAT subject is dispersed and has an irregular shape (Fig. 2-third row, left panel).In this vein, the S = [4 − 6]Hz region within the L region is adopted, defining more focused S region-related features.
• S Region BS-related Features 1) Relative Power (Power46(%)) of BS at S and L regions, i.e.: 2) Belongs at S (belongs46): This feature describes whether the majority of the BS density of a process belongs at the S interval.Takes only binary values, that is, 1 if the dominant peak of the BS domain is located at S band, and 0 otherwise.

2) CLASSIFICATION SCHEME
To assess the predictive ability of the introduced features in a supervised classification setting, a Leave-One-Subject-Out (LOSO) scheme was adopted, accounting for model hyper-parameter optimization and overfitting avoidance.Inside each LOSO loop, a standardization scaling is applied, followed by the ANOVA feature selection method, where the k best features are kept in each loop.For the proposed classification methodology, four types of classifiers are utilized, i.e., Gaussian Kernel Support Vector Machine (RBF-SVM), Linear Support-Vector Machine [37], k-Nearest Neighbors [38] and Logistic Regression [39].This procedure is repeated for 100 different random states.The Receiver Operating Characteristics (ROC) Curve [40] and the area under the ROC curve (AUC), as well as the Accuracy metric were used to evaluate all subjects against tremor level, and Sensitivity/Specificity values were calculated for all cases.receiving DBS to alleviate PD symptoms.During the tremor recordings, all patients had to be in the ''off'' condition, i.e., without any effect from medication (dopaminergic therapy, L-Dopa) for at least 12 h before the recordings.A more detailed description of the patients' demographics and their DBS/PD-related scores is given in Table 1.It should be noted that the MDS-UPDRS-III rest tremor subitem scores in Table 1 refer only to the ''off'' condition [8].We adopt the notation introduced in the original paper of Beuter et al. [8], that desribes the following tremor conditions: 1) 'ron' condition (DBS ''off'', medication ''on''): After refraining from medication at least for 12 hours, participants took 150% of their morning dose of L-Dopa.2) 'ref' condition (DBS ''on'', medication ''off''): Participants were recorded after refraining from medication at least for 12 hours, with the stimulation being in the ''on'' state.3) 'off' condition (DBS ''off'', medication ''off''): Participants were recorded without the effect of medication for at least 12 hours, after the DBS had been ceased for at least 1 hour.During this hour, rest tremor measurements were taken every 15 minutes for 60s, starting immediately after DBS cessation, resulting in 5 different ''off'' conditions, labeled as 'rof','r15of','r30of','r45of','r60of'.4) 'ren' condition (DBS ''on'', medication ''on''): The recordings of this condition were taken after participants had taken their medication.Stimulation had been in the ''on'' state for over an hour.It should be noted that the time sequence according to which the conditions were recorded is as listed above.As some patients were experiencing discomfort, not all patients were recorded in all conditions [8].

B. PREPROCESSING AND IMPLEMENTATION SETTINGS
All signals were recorded with a sampling frequency of 100 Hz, and had approximately 60s duration.As some signals had slightly longer duration, the first 60s from each recording were used for the analysis.Pathological PD rest tremor is thought to be located in the [4 -6] Hz region [41].Nonetheless, there are studies where the rest tremor band is taken over a broader region, e.g.[3.5 -7.5 Hz] [42], [4 -8] Hz [43], whereas physiological tremor is identified in the [8 -12] Hz band [42].In this vein, a FIR bandpass filter ([2.5 -12]Hz) was first applied to the recordings, to remove any noise originating from wrist or arm movements (< 2.5Hz) or any high frequency noise (> 12Hz).Since we have to deal with relatively low frequencies, a subsampling to 50 Hz was applied to the signals to reduce the computational cost.All signals were finally standardized to mitigate any effect due to amplitude variation.
In the process of BIC estimation (see ( 2)), values of P(f ) close to zero could create many spurious peaks in the resulting BIC.To account for this, usually a small constant ϵ is added to the denominator of BIC [44].In our case, the value of ϵ = 10 −7 was empirically selected, effectively eliminating many spurious peaks, without causing any alterations in the significant peaks of BIC.Moreover, to ensure that the estimated BS peaks have statistical significance, surrogate data sets that are explicitly Gaussian, linear, and share the same second-order statistics with the original data were produced.Estimation of the BS from the surrogate data defines the threshold that can be used to evaluate the significance of the BS peaks of the original data.Here we followed the computational implementation of this procedure based on [45]; the related hyperparameters were empirically set as P = 100 and M = 100.For the calculation of BS and BIC, the Higher-Order Spectral Analysis (HOSA) MATLAB Toolbox was used [46].

C. PERFORMANCE EVALUATION SCENARIOS
In order to evaluate the performance of the proposed approach, two specific evaluation scenarios were constructed, including features evaluation and classification pathways, as described below: In this scenario, we aim to assess the various treatment conditions of DBS and/or Medication and categorization into LAT/HAT classes under treatment.To that end, we use the 'ren', 'ref', 'ron' and 'rof' conditions (48 recordings from 12 subjects in total), where all possible combinations of treatment are involved (i.e., both DBS and medication, DBS only, medication only, neither DBS nor medication).To examine the effect of DBS, medication and tremor amplitude, we apply a statistical analysis for feature characterization and design a classification pathway, for each of these factors.

Features evaluation:
We use these 48 signals to test for statistically significant differences between the proposed features, utilizing as grouping variables, the presence of DBS, the presence of medication and categorization into LAT/HAT classes.Furthermore we test for statistical significance using as grouping variables the presence of DBS and medication, inside each amplitude tremor group.For feature characterization in this scenario, we employ two different statistical tests.The t-test is used as a parametric test, whereas to account for the small sample size and deviations from Gaussianity, the non-parametric Wilcoxon rank sum test is used.The threshold for the p-values to be statistically significant considered is set to 5%.In a comparative process, we compare the proposed features with those produced in [17], where Long-Term Correlation (LTC)-based features, i.e., the Hurst Exponent (H) and the MultiFractal Spectral Width (MFSW), along with Welch PS, combined with Principal Component Analysis (PCA) and kernel-PCA (kPCA)-based features are used.
Classification pathways: Three different classification pathways are implemented, considering the presence of DBS, the presence of medication and categorization into LAT/HAT classes.For each classification pathway, we use only features that were successfully evaluated at least in one statistical test in the features evaluation step, as well as the binary NoPeaks and belongs46 features.In the case where a condition does not present any significant feature, all features are used for feature selection.In order to enable a fair comparison of our results with the ones in [17], the same procedure as in [17] is followed in this scenario.In particular, only the RBF-SVM classifier is used, in a simple non-nested LOSO scheme, and the ROC curve AUC score (with the 95% CI), along with the sensitivity/specificity metrics are estimated.

2) SCENARIO B
In this scenario, we aim to assess the categorization into the LAT/HAT classes under treatment-free conditions.In this vein, we use only recordings that belong to the ''off'' conditions (i.e., 'rof','r15of','r30of','r45of','r60of'), where there is absence of treatment, and the tremor can be studied in its purest form.To examine the tremor evolution in time, as well as the time independent tremor characteristics, we separately study each ''off'' condition and the grand average of the five ''off'' conditions.In the latter case, data from only 11 subjects that participated in all 5 ''off'' conditions are used.To that end, we apply a statistical analysis for feature characterization and design a classification pathway, for each of the six conditions ('rof','r15of','r30of','r45of','r60of', grand average).
Features evaluation: In each of these six conditions, we test for statistically significant differences between the features, utilizing as grouping variable the categorization into the LAT/HAT classes.The non-parametric Wilxocon rank sum test is only employed here, due to the reduced the sample size.The features evaluation process is further expanded with a feature selection procedure, which takes place in the classification scheme.
Classification pathways: Six different classification pathways are implemented, considering the categorization into the LAT/HAT classes, as in the features evaluation procedure.For each classification pathway, we evaluate subjects against tremor amplitude level, using a LOSO scheme with nested 4-fold cross validation.We also insert an additional oversampling step before the LOSO loop, where the initial subjects are oversampled for the minority class, using the SMOTE algorithm [47].It should be noted that this oversampling procedure does not take place in the 'rof','r30of' and 'r60of' conditions, where the HAT subjects to LAT subjects ratio is 6/8, 4/8, 4/8 respectively.However in the 'r15of' and 'r45of' conditions, this ratio drops to 3/8.As we employ a LOSO scheme, there are certain folds where there is a big imbalance between the classes, which is combined with the fact that in such a case, the minority class contains only two subjects; hence, this clearly could negatively affect the classification results.
To maintain the integrity of the classification and, simultaneously, counterbalance the aforementioned issue, we oversampled with the SMOTE technique to a ratio of 4/8 for the 'r15of' and 'r45of' conditions, by generating an artificial subject.All classifiers and all evaluation metrics of the proposed methodology are applied in this scenario.For computational implementation of the classification schemes, the Scikit-Learn Python package is used [48].

1) FEATURE CHARACTERIZATION
Table 2 tabulates the results from the features evaluation in terms of the estimated p values (bold values denote p < 0.05).From this table it is evident that all analyzed features do not present any statistical power in the case of the existence or not of DBS.As is also stated in [17], this result contradicts the fact that DBS is commonly used for alleviation of PD symptoms in clinical settings.However, some of them yield statistically significant differences when the existence or not of medication (Med On/Off) or tremor amplitude (HAT/LAT) is used as grouping variable.This result applies to both our HOS features and LTC-PS features.Taking into account the use or not of medication, BS features, as well as the Power46 feature, seem to perform effectively under both t-test and rank-sum test, while BIC features provide statistically significant differences only when considering the t-test.LTC-and PS-based features also capture the effect of medication.Considering the characterization of tremor amplitude, we can conclude that HOS provide a higher number of features that are significant at least in one test, compared to the ones from LTC-PS.In addition, BicEnt1, BspecEnt1, and BspecEnt2 perform effectively at both tests, while for LTC-PS features only kPC1 appears significant in both tests.
Figure 3 depicts the distribution of the statistically significant features based on Table 2.In particular, Figs.3a and 3b illustrate the distribution of the {BicEnt1, BspecEnt1, TotalArea, Power46} and {BicEnt1, BspecEnt1, TotalArea} features for Med On/Off and HAT/LAT settings, respectively.These features satisfy the condition of successfully meeting at least one statistical test (Table 2).In Fig. 3a it is noted that the use of medication causes a dispersion of the bispectral content.This effect is reflected both in BspecEnt1 and TotalArea features.Moreover, the S Region BS-based Power46 feature confirms the expected effect of medication which is the decrease of power in the pathological tremor band.Finally, it can be seen that the ''Med On'' condition of the BicEnt1 feature translates into higher values of entropy in the BIC domain, indicating that more frequency components interact nonlinearly, producing a more dispersed spectrum.In Fig. 3b, the BicEnt1, BspecEnt2 and TotalArea features exhibit higher values in LAT case.As mentioned before, this translates into a more widespread content in the bispectral domain for the LAT subjects, compared to HAT ones.
Table 3 tabulates the results from the features evaluation in terms of the estimated p values (bold values denote p < 0.05), by examining the differences of each amplitude tremor group in medication and DBS settings.LAT subjects do not present any differences, both in the DBS and Med groupings, whereas in the latter, two statistically significant features are provided by LTC-PS.This poor performance of the features in the LAT group originates from the fact that LAT subjects are considered to have low tremor severity, both in terms of amplitude and MDS-UPDRS-III rest tremor subitem score (Table 1).On the contrary, statistically significant differences in features are seen for HAT subjects, both for DBS and medication settings.HAT subjects are thought to have high tremor severity, and thus DBS and Medication have a greater effect in the HAT group.The most important result of Table 3 is that a meaningful characterization about the use of DBS arises, contrasting the results of Table 2.More precisely, HAT subjects manifest statistically significant differences when  the state of DBS is considered.This result is quite important, since it is delivered only through HOS features, whereas LTC and PS features cannot capture this difference.This result is visualised in Figure 3c, where the distribution of the statistically significant features, TotalBic and BicEnt2, is highlighted.We can conclude that both the quantity of the nonlinear interactions and the coupling strength of those interactions, are increased when the DBS is in the ''On'' state.Taken altogether, statistically significant differences that arise in Table 2 regarding DBS and medication settings, are probably due to differences in the HAT class (Table 3).

2) CLASSIFICATION
Table 4 tabulates the best classification results (AUC (95% CI), Sensitivity/Specificity) generated by the proposed HOS-based approach and compares them directly to those from [17].As it can be seen from Table 4, our methodology provides significantly better results than the LTC-PS approach in all problems considered.Regarding the Med On/Off case, we observe that HOS-based features achieve an almost perfect discrimination, which is expected since medication strongly affects nonlinear properties, as was indicated by the feature characterization (Table 2, Table 3).In the DBS On/Off case, an expected decline in performance occurs compared to Med On/Off, confirming that medication is more effective in suppressing rest tremor.Finally, our methodology delivers a satisfactory performance in predicting amplitude level.This result is quite important, since it shows that through HOS-based features, tremor severity can be predicted, even when factors, like DBS medication, are on effect.
Figure 4 shows the selection percentage of the k best selected features yielded by the feature selection procedure, where k represents the feature dimension used in each classification setting as presented in Table 4.In the Med On/Off setting, the BS-based features noPeaks, BspecEnt1 and the S Region BS-based Power46 feature have a selection percentage of 100%, while BspecEnt2 and maxArea are selected as frequently as 70% and 20% of the cases, respectively.Furthermore in the HAT/LAT classification, BS-based features BspecEnt1, BspecEnt2 and BIC-based feature BicEnt1 are always selected, while the BS-based feature noPeaks has a selection percentage of 85%.Finally, in the DBS On/Off setting, BIC-based feature BicEnt2 appears in 100% of the cases, while the BS-based ESR and BIC-based BicEnt1 are selected almost 100% of the cases.

1) FEATURE CHARACTERIZATION
In Table 5, the statistical significance of the proposed HOS-based features, as well as their selection percentage, when the k = 8 best features are considered, are tabulated.This percentage is taken over 100 iterations of the LOSO cross-validation process.It is worth noting that some features provide information only for certain ''off'' conditions.This observation implies that the phenomenon that occurs after cessation of DBS and in the absence of medication, is dynamical in nature, and each ''off'' condition has its own properties.The time evolution of this phenomenon delivered through the recordings under ''off'' conditions is further analyzed in the next section.Nevertheless, the fact that features like BS and BIC entropy provide predictive power in all ''off'' settings, highlights even more the value of HOSbased analysis.4. TABLE 5. Scenario B: selection percentage and statistical significance of HOS features.

2) CLASSIFICATION
Table 6 presents the classification results under the Scenario B. More specifically, Table 6 tabulates the Accuracy, AUC and Sensitivity/Specificity scores, along with the 95% CI where applicable, for the three classifiers, at each ''off'' condition, along with the grand average case.It is evident that HOS-based features provide a robust classification with steady performance in the ''off'' conditions, as well as in the grand average problem.Observing the accuracy scores, it can be seen that in all conditions and all three classifiers, an average of one subject is misclassified, with the exception of the kNN classifier in the grand average scheme, where in some random states two subjects are misclassified.After inspection of the results, we noticed that the subject that is misclassified in all six problems, is the s16 individual.To further investigate this phenomenon, an unsupervised clustering analysis was applied to discover hidden patterns in the data, that might explain this misclassification.This analysis is presented and discussed in the succeeding section.

VI. DISCUSSION
A novel HOS-based approach towards the understanding of the PD rest-tremor characteristics has been attempted, taking into consideration two commonly used methods for PD rest tremor alleviation, i.e., the DBS and the L-Dopa medication.From the presented results it is derived that HOS-based features are capable of efficiently discriminating between the LAT and HAT classes.In the examined Scenario A, we utilized statistical tests and revealed how medication and DBS affect the tremor producing system.This effect is much more prominent in the HAT group.In addition, an increase in the number and strength of nonlinear interactions was observed through BIC features, in the HAT group during the DBS ''On'' setting.HOS-related features also yielded significant discrimination when tremor amplitude (LAT/HAT) was considered.
Comparing the findings of the proposed HOS-based analysis with the ones based on LTC-PS introduced in [17], it is evident that the HOS-based features allow for more consistent discrimination between the examined classes, providing further new insights about the effect of DBS on the tremor production system.It should be noted though, that HOS analysis failed to detect the effect of Medication in the LAT cohort.However, in Scenario A, our HOS-based features outperformed LTC-PS representations in all categories, and indicated the influence of medication on nonlinear properties, and the superiority of medication over DBS.A significant performance was observed in tremor amplitude level prediction through HOS-based features, even when medication and DBS were on effect.
In Scenario B, we investigated further the predictive ability of our HOS-based analysis on tremor amplitude, and implemented a nested LOSO classification scheme focused on the ''off'' conditions.The performance of the classifiers VOLUME 11, 2023 used, reached an accuracy score of 0.93 in the 'rof' condition and a ROC AUC score of 1 in the 'r45of' and Grand Average cases.Moreover, the classification performance remained high, irrespectively of the time passed after the termination of DBS process.

A. UNSUPERVISED NONLINEARITY-BASED CHARACTERIZATION
As mentioned earlier, in all examined classification problems, a certain individual was misclassified.In order to analyze this behavior and examine PD rest tremor properties from a new perspective, we implemented an unsupervised clustering analysis.The absence of the class information that characterizes this analysis, could contribute to a better understanding of rest tremor, as projected in the HOS domain.Since the signals were standardized during the feature extraction procedure, it is expected that the amplitude factor of a signal will not have significant contribution to the BS and BIC characteristics.Instead, it is interesting to see if new, nonlinearity-based classes can arise from the bi-frequency domain.To that end, we performed a clustering using the K-Means algorithm with a multiplicative penalty scheme on the mean value of the Silhouette metric, to determine the optimal number of clusters ranging from two to six clusters [49].The best clustering was found by searching over 100 random states for the initialization of the centroids of the K-Means algorithm.
In Table 7, the clustering results are presented.The resulting new class of each subject is noted for each ''off'' condition, and a comparison with the HAT/LAT classes is also made through the Simple Matching Coefficient (SMC), while the performance of the clustering is evaluated through the mean Silhouette score.The SMC ranges from 0 to 1, corresponding to zero similarity and perfect similarity, respectively.The Silhouette score measures the similarity of an observation with its own cluster (cohesion), compared to other clusters (separation).The Silhouette ranges from −1 to +1, where a high value translates into a sample that matches well its own cluster and at the same time poorly matches other clusters.As it can be seen from the results tabulated in Table 7, the estimated SMC takes high values and we can conclude that the tremor amplitude characteristics are related to the bi-frequency derived characteristics.Moreover, it is evident that the mean Silhouette score indicates a good clustering performance in all conditions.
Moreover, Table 7 indicates that the optimal number of clusters is two, in all cases.Another important outcome is that the s16 subject that previously was an outlier for the LAT class, here joins the HAT group (class ''1'').It should be noted that the initial class of the s16 subject (LAT), is well defined by the clinicians that performed the experiment [8].However, the adoption of the classification metrics from the nonlinear space, could shed more light upon the rethinking of the classification process, instead of solely relying on the commonly used approaches of amplitude and frequency-based representations, of the tremor signal.Perhaps, consideration of the tremor signal representation at the HOS domain, could provide insights that ameliorate any propensity for misclassification.
Table 8 provides an illustrative and explanatory description of the new classes, by presenting values of some representative features in the Grand Average case of the ''off'' conditions.Subjects that belong in class ''1'', are characterized by low values of Bispectral and Bicoherence Entropy, which, in fact indicate, a smaller quantity of nonlinear interactions, compared to subjects of the class ''0''.This conclusion is reinforced by the Total Bicoherence values, where class ''1'' subjects present lower BIC values, and, hence, reduced quadratic coupling strength of nonlinear interactions.Power46, which constitutes a feature that represents the physiology of tremor and has higher practical interpretability, gets very high values in the case of class ''1''.The behavior of the Power46 feature regarding the new classes, matches quite well to the other three features that describe nonlinear interactions.This is an important finding, since it points out that nonlinearity is a fundamental property of PD rest tremor.Based on the aforementioned observations, we move forward and define the classes ''1'' and ''0'', as Low nonLinearity System (LnLS) and High nonLinearity System (HnLS), respectively.These new classes are more connected with the nonlinear behavior of the tremor production system, rather than the amplitude/frequency characteristics of the tremor signal itself.Such a characterization is quite important, as it reveals a new perspective in approaching PD rest tremor, which cannot be captured with conventional approaches that solely examine amplitude and frequency characteristics.
Previous literature is indicative of tremor signatures that are referred to in general terms as ''nonlinear'', regarding untreated tremor and tremor during medication and DBS effect.Timmer et al. used nonlinear dynamics and claimed that tremor signals are indeed nonlinear oscillations and not strictly periodic, in fact possessing nonlinear second order stochastic dynamics [50].Razlighi et al. examined the effect of different DBS settings on PD tremor properties through nonlinear dynamics and found that tremor exhibited chaotic behaviors, and that these behaviors changed with alteration of DBS voltage settings [51].Finally, Sarbaz et al. provided very interesting results on the same dataset that we use here [8], regarding the nature of PD rest tremor through nonlinear analysis; they characterized untreated tremor as chaotic, whereas they noted that effective treatment changes this behavior towards a more stochastic manifestation, thereby weakening chaos [35].
As our insights stem from the conditions where treatment is in the ''off'' states, we can attempt a matching with the insights from Sarbaz and Pourakbari, that characterized these states as chaotic.Moreover severity of the disease increases chaotic behavior [35].In our case, severity of the disease is measured by the MDS-UPDRS-III tremor subitem, which is generally higher in the HAT subjects (see Table1), that mostly make up the LnLS class.Therefore, severity of the disease in nonlinear terms (LnLS class) could have some relation to the chaotic behavior described by Sarbaz and Pourakbari, whereas the HnLS class could tend to a more stochastic behavior.These assumptions are supported by the nonlinear measures Higuchi Fractal Dimension and Approximate Entropy, which had lower and higher values, therefore lower/higher nonlinearity, in the chaotic and stochastic states, respectively [35], as in the LnLS/HnLS classes.Although we cannot make conclusions on the relation of our nonlinearity insights to the chaotic and stochastic insights of Sarbaz et al., it is worth commenting on the similarity of our theories, which should be addressed in future research efforts.

B. DYNAMIC EVALUATION OF TREMOR PROPERTIES
So far, the various conditions have been independently examined from each other.In a step further, it is interesting to examine if the proposed HOS-based features can accurately describe the dynamic change of tremor properties over the time course of the two-day experiment, where medication and DBS conditions alternate between the ''On'' and ''Off'' states.Results from this dynamic analysis are shown in Figure 5, where the most representative HOS-based features are displayed under the available experimental conditions, as they change over time, in the LnLS subject cohort.To begin with, the superiority of medication over DBS is depicted in the Power46 and BS and BIC Entropy features.Lower power and higher entropy values are observed in the Med ''On'' -DBS ''Off'' condition compared to the Med ''Off'' -DBS ''On'' condition.After the DBS is turned off, a gradual deterioration is observed, reaching its peak 15 minutes later, as was also reported in [8].30 minutes after DBS was stopped, a rebound of the DBS effect occurs, followed again by a gradual aggravation until 60 minutes.This behavior after DBS cessation has been previously examined by other studies, and seems to be varying depending on subject condition and study protocol.In [52], where the essential tremor is studied, a rebound of tremor is observed two minutes after the DBS closure, declining until 30 minutes, when it then disappears.Hence, it could be possible that the condition that is described in 15 minutes, was more intense a few minutes earlier.In [53], the essential tremor rebound is reported at approximately eight minutes after DBS was turned off.Therefore, it is not safe to draw general conclusions about the exact behavior of tremor after DBS is switched off.Nevertheless, it is certain that ''off'' conditions have distinct differences with the other conditions.In contrast to the other features, Total Bicoherence indicates that the effect of DBS is stronger than the effect of Medication, since the DBS ''On'' -Med ''Off'' condition has the greater separation, when compared to the ''Off'' conditions.This might seem strange; however, it is confirmed in Table 3 from the p-values in the DBS ''On/Off'' test for the HAT group, where a significant p-value occurred when the Total Bicoherence was considered.This, combined to the fact that Bicoherence Entropy remains lower in the DBS ''On'' -Med ''Off'' condition, is another indication that DBS increases the coupling strength of nonlinear interactions.It should also be noted that we chose to demonstrate these results for the LnLS class only, where the tremor severity is at its highest.Such a decision is justified by the fact that in this approach we handle only HOS-based features, so studying the LnLS class seems more suitable than the HAT class.
From this dynamic perspective of the PD rest tremor, it can be concluded that HOS-based features deliver a compact description of the sequence of conditions recorded in this dataset, and provide variables that can be used as biomarkers for real time assessment of PD severity and behavior.Moreover the results of the scenarios examined in this study, indicate that HOS possess high sensitivity in discriminating the effect of medication and DBS.Thus, HOS could be appropriate for incorporation in an adaptive treatment system, where medication dosage and DBS parameters are dynamically adjusted, according to PD tremor condition.It is also important to note that the LnLS and HnLS classes introduced here, should not be considered as a replacement, but as a complement to the HAT/LAT classes.

C. COMPARISON WITH PREVIOUS WORKS
A comparison of HOS-based approach to previous studies conducted on the same dataset, highlights the superiority and the comprehensiveness of our analysis.Specifically, when comparing the classification results for scenario A to those from [17], we observe improvements of 9%, 19%, and 22% for the Med On/Off, DBS On/Off, and HAT/LAT classification tasks, respectively.Furthermore, in scenario B, our classification results can be directly compared to those in [36], where the goal was to distinguish between HAT and LAT patients using recordings from the 'rof' condition.While their classification results match ours (misclassification of one subject in both cases), the absence of classification results during the rest 'off' states ('r15of', 'r30of', 'r45of', 'r60of') and in states where medication and/or DBS are still active, hinders a thorough comparison.Finally, in [35], it was demonstrated that effective treatments (DBS and medication) lead to more noticeable stochastic behavior.While direct comparisons between our approaches are challenging due to methodological differences, this enhancement of stochastic behavior in the presence of treatment may be related to the increased nonlinearity observed in our study.

D. LIMITATIONS
Despite the promising results of the proposed approach, some limitations could be identified.In particular, the number of cases examined in both Scenarios can be considered small.However, in order to reduce the effect of the small sample size to the generalization power of the findings, the feature characterization was carried out with two different statistical tests that account for such limitation, whereas, in all classification settings, a LOSO scheme was implemented.Moreover, in Scenario A, only the statistically significant features enter the classification pipeline.We also believe that the potential of the new tremor classes that are introduced here could be further reinforced in a larger subject cohort.Moreover, the ''off'' conditions available in this dataset are recorded with a time distance of 15 minutes between each other.However, when tremor level changes are considered, this time interval is quite large.As previously mentioned, other studies [52], [53], have examined tremor behavior after DBS cessation with much smaller time step than the one available here.Apparently, examination of the HOS-based features under various time steps between various conditions could lead to more conclusive results regarding their ability to clearly capture the dynamic behavior of PD tremor before and after DBS and medication treatment.Research efforts to further address such limitations have already been undertaken.

VII. CONCLUSION
Parkinsonian Rest Tremor characterization and identification under different treatment conditions has been presented here, by involving, for the first time, features from the Higher-Order Spectra domain.When this approach was applied to experimental data with different PD tremor amplitude/frequency characteristics and under various treatment conditions, new insights and discrimination potentialities were identified.Specifically, HOS-based features responded with high classification performance that outperformed existing methods applied to the same dataset, in two different analysis scenarios; in Scenario A, three classification problems to examine the effect of DBS and dopaminergic treatment were addressed, whereas in Scenario B six different classification problems related to tremor in untreated conditions were solved.Through these results, our proposed approach manages to connect the clinically valid amplitude/frequency classes of LAT/HAT to identification insights regarding differentiation in treated and untreated tremor states.Most importantly, the proposed HOS-based approach introduced an additional perspective in PD resting tremor analysis via the revelation of nonlinear behavior as quantitatively expressed in the HOS-based features.This can allow rethinking of the classification metrics and can reveal transitions between various degrees of nonlinearity, when alternating between different treatment conditions; hence, a more efficient evaluation of the DBS and/or Medication effect on PD tremor characteristics can be accomplished.Future applications can examine the possibility of incorporating HOS-based nonlinear biomarkers as inputs to smart AI-powered wearable technologies to capitalize on their explainability and sensitivity of PD tremor properties.

FIGURE 1 .
FIGURE 1. Flowchart of the proposed HOS-based methodology.Initially, the accelerometer data undergo preprocessing.Following, the bispectrum and the bicoherence are calculated and subjected to a feature extraction pipeline, designed to capture tremor related information.Consequently, two different scenarios are implemented for characterizing and classifying treatment (medication/DBS) effectiveness and tremor severity (LAT/HAT).
In order to classify rest tremor signals into the HAT/LAT classes, different features that are drawn from the BS/BIC domains are proposed.The motivation behind is to accurately describe the morphology and the unique characteristics of the tremor affected frequency bands, as represented in the bispectral domain.For this purpose, different features are designed, to account for the form, the amplitude, and the distribution of bispectral content in the tremor affected region.In Fig.2, an example of the representation of the acquired signals from subjects with LAT (Fig.2-left panel) and HAT (Fig.2-right panel), at the PS, BS, and BIC domains is presented.The noticeable morphological changes between the two classes seen at these domains justify the selection of the features described next.

FIGURE 2 .
FIGURE 2. Standard procedure starting from the preprocessed signal (first row), to the PS (second row), to BS (third row) and BIC (fourth row) for the s8 HAT subject in the 'r15of' condition (right panels) and the g13 LAT subject in the 'r15of' condition (left panels).It should be noted that in the top panels the signals have been normalised.

6 )
Ellipse Similarity Ratio (ESR): After the main body of the BS content is located, an ellipse (E) is fitted to the contour (C) of the BS content.Then, the ratio of the intersection to the union of E and C is estimated as follows:

FIGURE 3 .
FIGURE 3. Visualization of HOS features under Med and DBS settings and amplitude tremor grouping.In subfigure 2a, HOS features accurately describe differences between Med states.In subfigure 2b, Bicoherence and Bispectrum Entropies successfully separate the high and low amplitude tremor groups.In subfigure 2c, Bicoherence features deliver a significant difference between DBS states in the HAT subjects cohort.

FIGURE 4 .
FIGURE 4. Scenario A Feature Selection Results: Percentage of times each feature is selected in each classification setting, for the corresponding feature vector dimension, as presented in Table4.

FIGURE 5 .
FIGURE 5. Evolution over time of (a) BspecEnt1 and BicEnt1 and (b) TotalBic and Power46, for all LnLS subjects.The blue line represents the mean value across LnLS subjects feature values.The red dashed line represents the 95% Confidence Interval of the mean, calculated over 100 bootstrap iterations.

TABLE 4 .
Scenario a classification results.

TABLE 6 .
Scenario B classification results.The 95% CI was less than 0.02 in all cases.

TABLE 8 .
Grand average of feature values for each subject in the ''Off'' conditions.