A novel plantar pressure analysis method to signify gait dynamics in Parkinson’s disease

: Plantar pressure can signify the gait performance of patients with Parkinson’s disease (PD). This study proposed a plantar pressure analysis method with the dynamics feature of the sub-regions plantar pressure signals. Specifically, each side’s plantar pressure signals were divided into five sub-regions. Moreover, a dynamics feature extractor (DFE) was designed to extract features of the sub-regions signals. The radial basis function neural network (RBFNN) was used to learn and store gait dynamics. And a classification mechanism based on the output error in RBFNN was proposed. The classification accuracy of the proposed method achieved 100.00% in PD diagnosis and 95.89% in severity assessment on the online dataset, and 96.00% in severity assessment on our dataset. The experimental results suggested that the proposed method had the capability to signify the gait dynamics of PD patients.


Introduction
Parkinson's disease (PD) is a common neurodegenerative disease [1][2][3], which is responsible for a considerable amount of disability and deaths, leading to extremely high demand for health resources [4]. Globally, PD patients are probably much higher than currently treated [5].
The clinical symptoms of PD include tremor, rigidity, and bradykinesia [6]. Meanwhile, gait disturbance is a common and debilitating symptom of PD [7,8]. Although classic motor symptoms occur early and are the pillars of diagnostic criteria, the development of gait disturbance drives the progression of the disease [9,10]. Gait analysis is essential for PD diagnosis and severity assessment.
There are many gait analysis studies in PD based on motion capture, inertial measurement units, force plate and force insoles [11][12][13][14]. Force insoles have advantages of convenient installation. Furthermore, the plantar pressure is crucial for gait analysis since the only external force during walking is the ground reaction force. The typical pathological manifestations of PD affect gait performance, which can be represented in plantar pressure.
In recent studies, much attention was paid to plantar pressure. Some studies focused on the analysis of classical kinematic features. Time and frequency domain features were extracted from plantar pressure signals to distinguish PD patients and healthy controls [15,16]. PD patients' gait kinematics and dynamics features had distinct properties [17,18]. However, the dynamics features had not been well explored. The gait dynamics had essential features which were used to quantify disease progression [15,16,19].
This work proposed a plantar pressure analysis method to signify PD gait dynamics. A sub-regions division method was designed to ensure that the plantar pressure collected by different experimental equipment can be uniformly processed. A dynamics feature extractor (DFE) was established based on three-dimensional phase space reconstruction (3D-PSR) and an improved complete ensemble empirical mode decomposition with adaptive noise (ICEEMDAN). The intrinsic mode functions (IMFs) were extracted and composed of the gait dynamics feature. The radial basis function neural network (RBFNN) was trained to learn and store gait dynamics. Moreover, a minimum error mechanism for classification was proposed. An online dataset was used to evaluate the effectiveness of the proposed method. In the PD diagnosis experiment, the accuracy achieved 100.00%. The classification of PD's severity was designed according to H&Y scores. Experimental results showed that the proposed method achieved 95.89% and 96.00% classification accuracy on online and our datasets, respectively.

Preprocessing and sub-regions division
Pathological analysis and clinical evidence of PD suggest that plantar pressure can indicate gait disturbance. During walking, plantar pressure transitions from heel to toe, and changes in plantar pressure on the left and right sides have a specific regularity unique to each individual. Analyzing sub-regions plantar pressure can provide valuable information on the subject's gait dynamics.
We proposed a method of dividing the plantar pressure signals into sub-regions. The plantar pressure signals on the left and right sides were divided into five sub-regions from heel to toe, as illustrated in Figure 1. These sub-regions were defined as follows: • Toes (Sub-region 1, R1): The region of phalanges in front of the phalangeal metatarsal joint.
This sub-regions division method was based on the anatomical structure of the foot and the specific characteristics of gait analysis in PD. It allowed for a more detailed analysis of the plantar pressure during gait.  During the preprocessing step, each test's first and last two gait cycles' plantar pressure signals were excluded to eliminate any possible start-up and stopping effects. Additionally, the plantar pressure signals were normalized with respect to the subject's body weight to account for individual differences in weight distribution during gait.
The plantar pressure signals of a PD patient and a healthy subject are shown in Figure 2.

Plantar pressure analysis method
We proposed a plantar pressure analysis method to signify the gait dynamics of PD patients. We divided the plantar pressure into five sub-regions and applied DFE to extract the gait dynamics features. The extracted features were set as the input signals, and the derivative of the features was set as the output signals. We used a constant RBFNN to learn the features. Based on the approximation error, we established a classification method. The flowchart depicting our proposed method is shown in Figure 3.    Figure 3. The flowchart of the proposed gait dynamics analysis method.

Dynamics feature extractor
A DFE based on 3D-PSR and ICEEMDAN was designed to extract the gait dynamics features. The 3D-PSR is an effective method for searching gait dynamics features [20,21]. The Euclidean Distance (ED) is a practical feature extraction method. It is also an effective method to calculate the distance between feature points and origin in phase space [22]. The ED from the points in phase space to the origin was calculated to obtain the transformed signal.
In this study, the time series were written as Q = {q 1 , q 2 , q 3 , . . . , q K }, the phase space can be reconstructed as follows: where τ is the time delay, i = 1, 2, . . . , K − (d − 1)τ, K is the length of the series and d is the embedding dimension.
The distance vector D = {D 1 , D 2 , . . . , D i , } between the point R i and the origin point can be defined as Eq (2.2).
In gait analysis, 6-10 gait cycles were appropriate [23]. The plantar pressure signals from 10 gait cycles were used in this study, and K is the length of the plantar pressure signals. Previous studies founded d = 3 can best represent the attractor for human movement [21,24,25]. In this study, d = 3 was used. It has been shown that when the signal has a strong periodicity, the time of the first zero of the autocorrelation function can provide a good embedding. This lag is approximately the same as the one-quarter length of the average cycle [26]. Cause the natural gait frequency was 1Hz and the sampling time of the force insole was 0.01s, the time delay τ = 25 was set.
The example of a PD patient's 3D-PSR of five sub-regions is shown in Figure 4. The distance vector D was extracted as the transformed signal. The curves in each sub-region are shown in Figure 5. The ICEEMDAN can extract the main components as IMFs. The IMF component of the most critical characteristics of gait cycles will have most of the energy. In this study, the periodogram was used to estimate the power spectrum of the plantar pressure signals. The Fourier transform transforms the plantar pressure signals from the time domain to the frequency domain. It estimates the signal's power spectrum by computing the transformed result's squared modulus. The rectangle window was used, and the number of points used in the Discrete Fourier Transform was set as the next power of 2, greater than the signal length K. This fully utilizes the Fast Fourier Transform algorithm to achieve more efficient computation.
The power spectrum of these IMFs was calculated, and the IMF with most of the energy was extracted as the gait dynamics feature. Compared with EMD, ICEEMDAN can effectively solve the problem of the existence of residual noise and spurious modes.   10 -3 PD sub-region 5 Figure 6. The power spectrum of IMFs in five sub-regions of a PD patient and a healthy subject.
After decomposing the vector D with ICEEMDAN, which was extracted from sub-regions plantar pressure, the power spectrum of the IMFs was calculated as shown in Figure 6. The frequency of the extracted IMF, which has the most energy, is close to the frequency of the gait cycle. The IMF with the most energy im f p for each sub-region forms the gait features as Eq (2.3) for the following task.
It is essential to consider the gait of bilateral lower limbs in gait analysis because of the asymmetry of PD. This study considered both sides' plantar pressure signals after the preprocessing and sub-region division. The dynamics features were extracted from the left and right plantar pressure signals. In the gait feature vector x, row 1-5 is the left side's feature, and row 6-10 is the right side's feature.

Classification mechanism
The classification mechanism we propose is that if two subjects have similar gait dynamics, then one subject's RBFNN can fit with the other subject's gait dynamics with a small error. We define the label of the test set by using the label of the RBFNN that can provide a minimum error in the training set.
Consider the general gait dynamics in Eq (2.4). The gait dynamics system is a continuous function, and the RBFNN is superior in approximating continuous functions.
where x is the gait feature as shown in Eq (2.3), and p is a constant parameter. The RBFNN is as follows: where Z is the gait feature, W is the weight, N is the number of nodes, S is the Gaussian function as shown in Eqs (2.6) and (2.7).
x n ] is the state vector, A = diag[a 1 , . . . , a n ] is a diagonal matrix, a i > 0 is a constant. Moreover, neural weight is updated by the Eq (2.9). 10) The overall system is shown in Eqs (2.11) and (2.12).
Based on the study [27], the constant RBFNN can approximate gait dynamics as shown in Eq (2.13).
The RBFNN consists of three layers: the input, hidden, and output layers. The input layer consists of 10 units corresponding to the input feature x. The weights between the neurons from the input and hidden layers are all 1. The neurons in the hidden layer use radial basis functions as activation functions. The weights between the hidden and output layers were changed by training the specific net for each subject. Meanwhile, the output layer also has 10 units corresponding to the structure of the derivative of x. In the hidden layer, the maximum number of neurons was 100. The radial basis function's expansion rate is 1.
The training set consists of M gait features. As demonstrated above, the general gait system dynamics can be preserved in constant RBFNNWS (x). The M estimators of training gait features are constructed by using the learning knowledge acquired during the training stage, as shown below: x j =W j T S j (x j ) (2.14) where j = 1, . . . , M, x j = [x j1 , . . . , x jn ] T are the features of gait dynamics and n = 10. For the data in the test set, the gait dynamics features x t were extracted, and the state estimation error was defined as follows: where e j = [e j1 , . . . , e jn ] T is the test error vector, n = 10 is the number of channels. Furthermore, we compute each channel's L1 norm of the test error vector e j and get the sum value e t j as the classification indicator. The label for the subject in the test set label t is defined with the minimum error's label label m . For each subject in the training and test set, the input was the gait feature vector x, and the outpuṫ x was the derivative of the vector x. Each subject's gait dynamics were saved in an RBFNN in the training set, and the corresponding label was saved. The feature x was fed in the training set's RBFNN during testing, and the estimated outputx was calculated. We assume that the label corresponding to the network which can obtain the minimum estimate error is the label of the subject in the test set.

Dataset
The online dataset Physionet [28] is used. This dataset included 73 healthy subjects and 93 idiopathic PD patients. Their plantar pressure signals in the walking test are provided. The database includes subjects' vertical ground reaction force records as they walked at their usual, self-selected pace for approximately 2 minutes on level ground, which consists of about 100 steps, including start-up, ending, and turning [28,29]. Underneath each foot were 8 sensors (Ultraflex Computer Dyno Graphy, Infotronic Inc.) that measure the force (in Newtons) as a function of time. The output of each of these 16 sensors has been digitized and recorded at 100 samples per second. The plantar pressure signals is normalized by body weight in preprocessing. However, two healthy subjects' and two PD patients' weight information needed to be included, and their data were excluded from the experiment.
We performed a walking test and recorded plantar pressure signals in 63 patients with wireless force insoles (Podomed, ZIGUN AVOIN Company, Finland). The data acquisition rate was 100 Hz, and the measurement range was 100 N with 0.1 N/cm 2 resolution for each sensor. Each insole has 100 sensors, and the surface area of all the sensors is less than 6 × 8 mm.
We used the standard 10m walking clinical protocol. Each patient in our dataset walked with the protocol. When the subjects heard the beep, they started walking at their self-selected pace and stopped at the end. Patients did the walking test in the same room to reduce the impact of environmental factors, and patients were always under the protection of doctors through the walking test. The subjects provided their written consent to participate in the study. All 63 patients received the walking test and in-depth professional evaluation and treatment, and 50 patients were diagnosed with primary Parkinson's disease. The other 13 patients had some comorbidities and were excluded from the analysis.
Plantar pressure signals were collected at the Department of Neurorehabilitation, Tianjin Huanhu Hospital, Tianjin, China. This study was approved by the ethical committee of Tianjin Huanhu Hospital and registered in the Chinese Clinical Trial Registry (ChiCTR1900025372).
Gait analysis typically requires analyzing 6-10 gait cycles [23]. The start-up, ending, and turning gait were excluded in the preprocessing. Moreover, 10 gait cycles' plantar pressure signals were used both in the online and our datasets based on the requirements of gait analysis.
In this study, PD patients were divided according to H&Y scores. The division is based on the median, a useful way to divide the data. PD patients with a score of 2 or less are considered mild stage, and scores from 2 to 3 are considered moderate stage. Demographics of the datasets are shown in Table 1.

Classification results
The gait system dynamics were preserved in constant RBFNN, and the classification experiment relies on our designed classification mechanism. In RBFNN, the input is x, and the output isẋ. The size of x andẋ is the same: a matrix with 10 rows and K columns, where K is the length of the input signal. The classification process does not require retraining or weight updates of the network.
In classification experiments, 10-fold cross-validation was used. For the evaluation, five parameters, including accuracy (Acc), precision (Pre), sensitivity (S en), specificity (S pe), and F1 score (F1), are used. These measurements are defined as follows: The proposed method was used to conduct classification experiments in the online dataset. We conducted a comparative experiment to analyze the performance of different parts of the proposed method. We compared the classification performance in four cases. Whether the plantar pressure subregions are divided and whether DFE is carried out. If plantar pressure is not divided, then the total force is used for analysis. The experiments' results are shown in Tables 2 and 3. The results showed that the classification accuracy of the proposed method achieved 100.00% in PD diagnosis and 95.60% in severity assessment. Comparative experiments showed that each part of the proposed method improves classification accuracy. The experimental results showed that using the sub-regions signals can obtain higher classification accuracy than using the total pressure. This showed that sub-regions signal could provide more helpful information. Meanwhile, the classification accuracy can be effectively improved by the proposed DFE, whether the plantar pressure signals were divided or not.
Meanwhile, we also compared related studies on PD diagnosis and severity assessment on the same dataset from 2020. All these studies performed classification experiments for PD diagnosis, but only four performed severity assessments. The classification performance is shown in Tables 4 and 5. The accuracy of the proposed method exceeds most of the recent studies. The severity assessment experiment is applied to our dataset. The experiment's results of 10-fold cross-validation are shown in Table 6, with an accuracy achieved 96.00%. This shows that the proposed method is general, not for a specific signal acquisition device. To demonstrate the generalization ability of the proposed method, we conducted experiments using an online dataset as the training set and our dataset as the test set. The results showed an impressive accuracy of 90.00%. Conversely, when we used our dataset as the training set and the online dataset as the test set, the classification accuracy was 87.91%. The performance indicators are shown in Table 7.

Discussion
In experiments analyzing each part of the proposed method's contribution, we found that dividing plantar pressure into sub-regions effectively improved the classification accuracy. This might be because plantar pressure in different sub-regions had its specific characteristic, often masked by the total plantar pressure signals. Meanwhile, by extracting the dynamics features through DFE, the most effective part of the signals were extracted to achieve better classification accuracy.
Most previous studies had focused on the total force of plantar pressure and the data from specific sensors. The proposed method divided the plantar pressure into five sub-regions, which made this method general. The subsequent gait dynamics analysis scheme could be applied for different measuring equipment if the sub-regions are divided. The limitation of our method was that if the number of sensors were too small to be divided into 5 sub-regions, it would cause difficulties for the use of the sub-region division method. However, diagnosing PD using total plantar pressure signals achieved a classification accuracy of 92.58%.
In recent studies, deep learning and machine learning methods, such as DNN, ANN, CNN, and Random Forest, were used in PD diagnosis and severity assessment. The proposed method relied on RBFNN to model the gait dynamics of each subject and then classified the gait dynamics by comparing the output errors, which were network independent. Most of the recent works used 10-fold crossvalidation to test classification performance. Therefore, in our experiment, 10-fold cross-verification was carried out for better comparative analysis. The proposed method achieved higher classification accuracy in PD diagnosis and severity assessment.
The proposed method obtained high classification accuracy when the training and test sets were from different acquisition devices. The sub-regions division method increased its generalizability across various equipment and conditions. Additionally, the classification method based on kinematic features achieved high accuracy in PD diagnosis. However, the same method applied to severity assessment results in a significant decrease in classification accuracy. The proposed method achieved good classification accuracy in PD diagnosis and severity assessment.
The proposed method in this study achieved high accuracy in distinguishing healthy controls from patients with PD. However, a small number of patients in the early stages were included in the online dataset and our dataset. Patients in the early stages of PD were with mild motor dysfunction, and it was challenging to distinguish between healthy controls and PD in the early stages. In future studies, we intended to focus on collecting data and conducting a diagnosis of healthy controls from PD in the early stage.

Conclusions
This work proposed a plantar pressure analysis method to signify PD gait dynamics. The plantar pressure signals were divided into five sub-regions based on anatomy. Based on 3D-PSR and ICEEMDAN, a DFE was established. Gait dynamics features were extracted, and a classification mechanism based on the approximate error in RBFNN was proposed. Experiments were performed with an online dataset and our dataset. The proposed method achieved high classification accuracy in PD diagnosis and severity assessment. It had the potential for application to assess other diseases with gait disturbances related to temporal variation and instability in gait dynamics.

Use of AI tools declaration
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.