The signal processing method applied to the waist-worn accelerometer measurements to determine the ON/OFF state of PD patients relies on characterising motor symptoms, similarly to the methods presented in the related work section. In this sense, two specific algorithms, which analyse the presence of dyskinesia and bradykinetic gait, are used. Their output is then merged based on a hierarchical algorithm that eventually provides the motor state estimation. This section describes each one of the methods used to estimate the motor state of PD patients.
In a previous work of the authors [
53], a previous version of the ON/OFF detector was presented, which was only based on the analysis patients’ gait and it was validated in a different database than the one used. This new work presents a novel approach in which two algorithms are combined: a bradykinetic gait detector and a choreic-dyskinesia detector that, combined by means of a decision tree, perform the detection of the motor states. Notably, the algorithmic basis for the detection of bradykinesia in [
53] and this work, with some modifications, is based on a previous work by the authors [
52]. On the other hand the algorithmic basis of the dyskinesia block is based on another previous work [
33]. The way to combine these algorithms and the database where they are validated are completely new. In addition, the processing required for extending bradykinetic gait detection into 10-min periods and the method presented in
Section 4.3 to self-tune the corresponding thresholds are new.
Taking into account these observations, we aim to obtain a set of signal processing methods that identify the presence of bradykinesia and dyskinesia during daily life activities, without requiring patients to perform specific movements, as some methods presented in the related work do. Related to choreic dyskinesia, this movement disorder can be evaluated without requiring any exercise, since it is an involuntary repetitive movement that patients manifest in any body segment with specific frequencies that have been shown to be up to 4 Hz [
36,
55]. Bradykinesia requires a more complex solution since, in general, it is manifested as a slower than normal movement that may be presented in any body segment. In order to automatically evaluate bradykinesia, movements belonging to the activities of daily living are considered, since they can be automatically assessed through wearable sensors. In this sense, gait is one of the movements involved in many of these activities and, in addition, it is an automatic movement that is also performed slowly by PD patients due to the effects of the disease [
56]; hence, authors consider it the optimal way to analyse bradykinesia. This way, a signal processing method that analyses gait to determine the presence of bradykinesia is considered, which will enable the monitoring of low dopamine levels without requiring patients to do any specific movement. A general schematic representation of the algorithms is shown in
Figure 2.
4.1. Dyskinesia Detection
Dyskinesia detection is the first algorithm presented for the ON-OFF motor states detection. This algorithm was designed based on many previously works that relate dyskinesia to an increased power spectra of some specific frequency bands. The algorithm was developed by analysing the frequency spectra of inertial signals obtained from PD patients while performing different activities, either while presenting dyskinesias and without presenting them. As a result [
52], a specific frequency band in which power spectra increases with dyskinesia was identified; in addition, other activities were found to also increase the power spectra in some bands that are overlapped with the dyskinesia one; thus, these other activities could provoke false positive detections of the symptom. These overlapping frequency bands are also examined in the algorithm in order to avoid inaccurate dyskinesia detections.
Dyskinesia frequency band was identified within the range of 0 to 4 Hz [
57]. As mentioned previously, other activities with high power spectrum in the same band were found; for example, the natural frequency of gait and going upstairs and downstairs ranges from 0.5 to 6 Hz [
58,
59], being overlapped with the dyskinetic band. However, these activities have strong harmonics whose frequency reaches 20 Hz, which is not overlapped with the dyskinetic band. On the other hand, posture transitions span the band from 0 (not included) to 0.68 Hz [
60]. A symptom that could introduce harmonics at frequencies of interest is tremor; nonetheless, according to a consensus of the Movement Disorders Society [
55] the frequency of Parkinsonian tremor goes from 4 Hz to 9 Hz. This way, the upper limit of the dyskinesia band is set to the lowest frequencies of tremor, and the possible increment in the power spectra, caused by the tremor, would incorporate in the non-dyskinetic band. In consequence, dyskinesia algorithm relies on the calculus of three power spectra values: first, the one corresponding to dyskinetic band (
Pd), considered to be in the (0.68, 4] Hz range; second, non-dyskinetic band (
Pwalk), considered to cover [8, 20] Hz; and postural transition band (
PPT), which is (0, 0.68] according to [
57]. The power spectra in a given band is computed as the summation of the corresponding harmonic amplitudes among the three axis.
The sampling frequency has been determined by the maximum frequency of interest following the Nyquist-Shannon sampling theorem [
61]. Given a sample rate
f’s, the complete reconstruction of a continuous signal is guaranteed for a frequency band limit below
f’s/2. In consequence,
f’s is set to 40 Hz, since 20 Hz is our maximum harmonic of interest in the previous frequency bands and we want to minimise the resources used by the algorithms. On the other hand, the window length is set to 128 samples since it enables the evaluation of postural transitions (below 0.68 Hz) and dyskinesias. In consequence,
= 128/
f’
s = 3.2 s.
Dyskinesia features
Pd,
Pwalk and
PPT are obtained in each window of 128 samples, i.e., 3.2 s. They are used, as Equation (1) shows, to determine if a patient manifests dyskinesia, does not manifest it, or performed a movement that does not allow to evaluate its presence (i.e., the output is
Unknown). The latter case refers to a patient who walks or performs a postural transition, in which cases these movements do not enable the detection of dyskinesia since they have overlapped frequency bands. This way, the detection of dyskinesia in a certain window
is defined by
according to:
where
,
,
are the thresholds for dyskinetic and non-dyskinetic bands (posture transition band and walking band) respectively,
is the logical
OR operation and
is the logical
AND.
The values found for these thresholds are
= 1.75,
= 0.95, and
= 1. These thresholds have been set based on the previously described study [
52] (see
Section 3.3.2) and they are used in a generic way for any patient.
Dyskinesia is a movement alteration that commonly appears during several minutes. Nonetheless, short windows are being used (
= 3.2 s). In order to determine the presence of dyskinesia in a more appropriate time interval, it is proposed to collect the output of several windows under a period of
T = 60 s. Each window is overlapped with the previous one by 64 samples, i.e., a new window starts every half a window. The algorithm, thus, provides an output once per minute obtained from the information included in its
windows, being
the floor function, according to:
where
are the outputs represented in Equation (1) corresponding to the windows evaluated in minute
,
is the number of time windows in which the condition
was not held,
is the minimum rate of dyskinetic windows in the analysed period, and
tc is the threshold that represents the minimum rate of windows to analyse in order to validate the detection.
According to Equation (2), the algorithm output in a one-minute interval is determined to be dyskinetic () provided that most of the analysed window outputs are dyskinetic. This means that a minute period is considered dyskinetic if the rate of positive outputs of Equation (1) in this period is greater than . However, as these band’s power spectra might be increased by other activities, and not only by the appearance of dyskinesias, we must add a parameter of confidence in which we ensure that the patient is not performing activities that might cause false detections in this band. Then, a confidence index is defined as and represents the number of windows that were rejected due to the condition . A low confidence index indicates an unreliable detection assessment because only few windows could be analysed. For this reason, the confidence index is required to be greater than threshold tc.
Threshold values were found based on an optimisation procedure on the signals collected in the previously mentioned study with 20 patients [
52], in which several values for the thresholds were evaluated in 10 patients and the accuracy was measured onto the signals from other 10 patients. Values found were
= 0.4 and
= 0.3.
Figure 3 shows a schematic representation of the dyskinesia algorithm block.
4.2. Bradykinesia Detection
Bradykinesia detection analyses Parkinsonian gait, as previously presented. The signal processing algorithm exploits the fact that gait, as an automated movement, is slowed in Parkinson’s patients during low-dopamine level periods. In consequence, it is considered that the signal processing method has to, first, determine that patients are walking; second, identify gait cycles from the accelerometer signals; and, third, characterise gait cycles through a measurement that correlates to the presence of bradykinesia. The complete bradykinesia detection method consists of a five-step characterisation method, as described below.
The first step consists in detecting gait and it is based on an SVM classifier. SVM is chosen, first, given the bi-classification nature of the problem at hand (detecting if the patient walks or not), which matches the bi-classification problems that SVM solve; second, due to the high performance that SVM have reported when dealing with this kind of classification problems; and, finally, because SVM allows us to obtain a global optimal solution, as opposed to Artificial Neural Networks that may provide suboptimal ones. Thus, given the signal contained in a time window of
= 3.2 s, as in the dyskinesia case, it is needed to determine whether the patient is walking or not. The SVM is trained with a Radial Basis Function (RBF) kernel and vectors {
p1, …,
pn}, where
pi = [
] and
are the features which characterise the window of the accelerometer signal. Hyper-parameters γ and
C were tuned through Cross-Validation by testing the following values: {10
−3, 10
−2, …, 10
3}. In a previous study [
52], 800 frequency features were analysed using data from 10 patients. Two features were finally selected for walking detection as those that maximised inter-class distance and minimised intra-class distance according to Relief algorithm [
62]. These two features (
k = 2) selected for gait detection were the tri-axial power spectra between in the frequency bands [0.1, 3] Hz and [0.1, 10] Hz, which are noted as
h1 and
h2, respectively.
Each vector
pi has a label
yi = {1, −1} according to video observations used to label the video:
yi = 1 corresponds to those windows whose corresponding video labels were walking and
yi = −1 to the remaining windows. Dataset elements are denoted as {(
p1,
y1), …, (
pn,
yn)} and were employed to find the SVM classifier that allows the detection of gait by solving:
where
,
b is the hyperplane bias,
is the hyperplane that separates both classes and
are the slack variables. Parameters
C and γ are determined as the values that maximise the accuracy among the values 10
−2, 10
−1, …, 10
2 in a 10-fold Cross-Validation [
63], which were found to be 10 and 0.1, respectively.
The label of a new window
p is then obtained by:
where the set of α are the Lagrangian multipliers of the dual problem formulation of the SVM.
The second phase is focused on detecting patients’ strides and only those windows whose feature representation
pj = [
] satisfies
l(pj) = 1 (walking) are analysed. The principles used to detect strides are based on Zijlstra et al. work [
64]. Although segmentation techniques can be employed to detect strides [
65,
66,
67], we restrict to biomechanical properties of gait and the way they are observed in the acceleration signals to do so. More concretely, the beginning of the support phase of gait, that is when the heel touches the ground, can be detected by a local minimum in the front acceleration measured from the bottom of the trunk [
64]. This event in the gait cycle is known as ‘initial contact’ and it is regarded as the beginning of step. However, due to lateral particularities of PD [
1], our interest focuses on strides; i.e., the signal comprised between two consecutive steps of one feet. The discrimination between left and right steps can be performed by analysing the relative extrema of the lateral acceleration in the waist, which approximately describes a sinusoidal period during gait cycle [
64].
The third step consists in characterising the detected strides in the previous phase with the aim of analysing the presence of bradykinesia. The basis of this step is the previously mentioned study [
52] where several statistics were applied and evaluated in 20 patients. In this study, several features that characterised strides were analysed into their ability to linearly separate the presence of the interest symptom and, also, intuitively represent the fluidity of movement. From the conclusions of this work, the best feature to characterise the fluidity of movement was the power spectra in the bands of (0, 10] Hz of the stride. Given a stride detected on the accelerometer signal of a certain patient, the (0, 10] Hz power spectra of the stride is represented by
.
The considerations for the fourth step are related to the fact that bradykinesia is a symptom examined during gait, as it is an automatic movement. In consequence, with the aim of analysing gait during its highest degree of automation, the inertial parameters of gait are analysed after walking started and before the patient stopped. Therefore, we consider a sequence of
U consecutive walking windows, i.e.,
l (
px) = −1,
l (
px+1) = 1, …,
l (
px+U) = 1,
l (
px+U+1) = −1, called
walking stretch, during which
S strides are detected; first and last two strides are not considered because they are not done under a high automation control, i.e., only
4 strides are considered. Thus, the result of averaging the fluency characteristics obtained from the strides within a walking stretch
k is denoted by
and is defined as:
where
is the number of strides detected in the walking stretch
k.
Considering that, when bradykinesia appears, it remains present for long periods of time, we designed an aggregation strategy that provides an algorithm output per minute, in order to simplify the final evaluation of the presence or absence of symptoms. In this aggregation, algorithm’s output for a given minute h, noted as , is computed as the average of the bradykinesia values associated with the strides contained in the walking stretches within that minute. The standard deviation of these bradykinesia values is represented by .
Finally, the algorithm’s fifth step aims to give more robust outputs every minute by considering longer periods. In consequence, a weighted average of the last 10 min is proposed. Note that the output of the algorithm persists once per minute. Weighted aggregation is based on the importance each bradykinesia value has; for example, the value obtained in a minute from only 2 strides should not have the same weight as the value obtained in one minute with 20 strides detected. Furthermore, a very high standard deviation within a minute means a large scatter in the data and, in case of very high values of
, the average may not be significant. Furthermore, it has been observed that when the patient goes upstairs or downstairs the standard deviation grows above the usual values. An estimation of the maximum value of
has been empirically determined through the signals from the previous study [
52]. This value was determined by studying the values of standard deviation presented by patients who had made the dubious activities. From the presented premises, values per minute are filtered in the following way: those minutes
h in which
is greater than threshold 1.7, which is formalised through
coefficients that are set according to:
where
is the number of steps in the minute
h.The results presented in [
52] show that the algorithm works best when considering the average walking steps in stretches over 5 strides. From this result the value of the minimum number of steps in a minute to be able to estimate the presence of bradykinesia is 2. In order to take into account the number of strides, a weight function is included. The function selected is the sigmoidal, considering that the minimum weight (0) should be with 0 steps and the maximum weight is 1. Whereas the maximum number of strides in a minute will be around 30–40, it is considered that from 20 strides the confidence in this minute should be high and therefore the weight should be the maximum. According to these considerations, the weight of data
in a minute
h is represented by
. Finally, the fluidity value representing a period of 1 min, but weighted with the last 10 min, is calculated through the next function:
where
if
.
This value is finally used to determine the bradykinesia presence in the last-minute period under analysis. Diagnosis of bradykinesia is set differently for the first minute. The existence of bradykinesia evaluated for the first minute (
) is defined by:
where
is the patient-dependent threshold to determine the presence (1) or absence (−1) of bradykinesia. This threshold is unique for each patient and must be particularised, as it is described in next subsections.
From this first minute, in order to avoid constant changes of diagnosis in intermediate states, a minimum variation from the threshold must be considered. This minimum variation is consequently determined by the maximum allowable standard deviation. The presence or absence of bradykinesia for the next minutes (
j > 1) is set as follows:
Thus, the output of the bradykinesia algorithm in a given one-minute period
j, noted as
, is
whenever the patient did not walk in the corresponding minute,
in case of bradykinesia being detected, and
whenever not bradykinetic gait was present.
Figure 4 shows a schematic representation of the bradykinesia algorithm block.
4.3. Self-Adapting Bradykinesia Detection Algorithm
The threshold applied to the output of the bradykinesia algorithm allows determining the presence of the symptom by dividing the range of possible values into two zones, one for each motor state. However, the selection of the threshold is very critical. Given that bradykinesia values depend on the way of walking of each individual, a young person without any pathological movement would provide high values; nonetheless, with older patients and/or with the presence of diseases such as arthritis, lower values would be obtained. Similarly, ON and OFF motor states are very patient-dependant. In consequence, bradykinesia values from each patient must be analysed in order to establish an optimal separation threshold.
In our previous works [
15,
52], the adaptation of the threshold to each patient was performed by a customisation process which is, in practice, long and complex. This process requires that the patient visits the clinical setting without medication, a fact which is already difficult to accomplish in many cases, but also implies that the medical team should perform a double clinical examination and a double assessment as the patient performs a series of exercises (mainly walking) both in OFF and ON states. Arguably, it is a methodology completely inapplicable in the clinical practice. In this section, a new methodology to calculate the threshold is presented.
This new methodology is based on automatically analysing the distribution of the bradykinesia fluency values () obtained during few days. Ideally, their histogram could present two clearly separate distributions representing each motor state. In this case, the bradykinesia optimal threshold lies within the gap between the two distributions. This type of clear distributions is found only in some patients, but it does not occur in most of them, where the difference between states is not so obvious. In these cases some empirical rules, that allow optimally adjusting the threshold, are applied.
More specifically, bradykinesia weighted values are first collected during few days (from 1 to 3 days). A histogram is then obtained in order to analyse the data distribution. Histogram bins are arranged to cover bradykinesia values from 2 to 15 since, from our experience, fluidity values from PD patients in both motor states are contained within this range. Given this histogram, the special case in which two different distributions are found, i.e., one for each motor state, is determined by locating empty bins. In order to standardise it, it is considered that both distributions must be separated by at least 0.5 points and, in addition, both of them must contain at least 10% of the total data, in order to avoid identifying a double distribution from merely isolated data. For this case, the value of the threshold is set to that value in the middle between distributions.
However, the most common case consists of overlapping distributions. In this case, the premise for calculating the threshold for this group of patients consists in considering the lower values of the distribution corresponding to the OFF state and the highest ones to the ON state. Then, the threshold can be set based on the percentage of frequencies remaining on either side of the distribution. To implement this approach, the value of the histogram’s bin that has the largest absolute frequency, i.e., the mode, is obtained. Then, the bin that is located immediately below the mode and whose frequency is higher than 60% of the mode’s frequency, is selected as threshold .
4.4. ON/OFF Motor States Detection
PD patients manifest motor fluctuations as an alternation between ON and OFF states. As previously described, provided that specific symptoms and movement alterations appear in each motor state, a hierarchical algorithm is designed to estimate the motor state of PD patients by combining the output of the previously presented methods. This way, dyskinesia and bradykinesia algorithms’ outputs are merged.
In order to work in a time unit closer to the gold standard, which are the annotations given by patients that are commonly provided every 30 or 60 min, a similar time basis is proposed. This time unit is 10 min since it is considered long enough to give accurate estimations and short enough to avoid mixing different motor states. More concretely, the motor state classifier first computes the presence of bradykinesia and dyskinesia into the period of 10 min, noted as
and
, respectively, according to Equations (10) and (11):
where
counts the number of elements satisfying the within condition.
Bradykinesia-algorithm’s output is 1 when the symptom is present, −1 if it is absent, 0 if an intermediate state has been detected and Unknown (U) if there is not any gait period detected in the last 10 min. On the other hand, dyskinesia-algorithm’s output is 1 when the symptoms are present, 0 if absent and Unknown whenever the patient walked or performed a posture change in most of the 10 min, which is unlikely.
Once the 10-min output has been obtained, the motor state estimated by the algorithm at time
, which is noted as
, is defined by Equation (12):
OFF state is considered whenever a patient has a bradykinetic period. On the other hand, ON state is estimated if non-bradykinetic gait or dyskinesias have been detected. Finally, an intermediate state is also added as a consequence of the intermediate bradykinetic state. Once several consecutive outputs of the ON/OFF decision-tree defined in Equation (12) are obtained, a small filter is then applied. Considering three consecutive outputs of the decision-tree, if a blank period is found between two periods that are equal, the empty period is then set to the same state.
Figure 5 shows a schematic representation of the complete ON-OFF algorithmic block.