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

1.   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.

2.   Methods

2.1.   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:

Figure 1.  The schematic diagram of sub-regions division.

DownLoad: Full-Size Img PowerPoint

● Toes (Sub-region 1, R1): The region of phalanges in front of the phalangeal metatarsal joint.

● Forefoot (Sub-region 2, R2): The area around the phalangeal metatarsal joint.

● Midfoot (Sub-region 3, R3): The area between Sub-region2 and Sub-region5, divided equally.

● Rearfoot (Sub-region 4, R4): The area between Sub-region2 and Sub-region5, divided equally.

● Heel (Sub-region 5, R5): The region of the calcaneus.

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.

Figure 2.  The plantar pressure of a PD patient and a healthy subject. Above: The measured total plantar pressure during ten seconds. Bottom: The five sub-regions' plantar pressure.

DownLoad: Full-Size Img PowerPoint

2.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.

DownLoad: Full-Size Img PowerPoint

2.2.1.   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 = \left\{ {{q_1}, {q_2}, {q_3}, \ldots, {q_K}} \right\}$, the phase space can be reconstructed as follows:

where $\tau$ is the time delay, $i = 1, 2, \ldots, K - (d - 1)\tau$, $K$ is the length of the series and $d$ is the embedding dimension.

The distance vector $D = \left\{ {D_1}, {D_2}, \ldots, {D_i}, \right\}$ 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 $\tau = 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.

Figure 4.  3D-PSR in five sub-regions' plantar pressure.

DownLoad: Full-Size Img PowerPoint

Figure 5.  The ED curves of five sub-regions in a PD patient and a healthy subject.

DownLoad: Full-Size Img PowerPoint

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.

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 $imf_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.

Figure 6.  The power spectrum of IMFs in five sub-regions of a PD patient and a healthy subject.

DownLoad: Full-Size Img PowerPoint

2.2.2.   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).

The $\hat x = [{{\hat x}_1}, \ldots, {{\hat x}_n}]$ is the state vector, $A = diag[{a_1}, \ldots, {a_n}]$ is a diagonal matrix, ${a_i} > 0$ is a constant. Moreover, neural weight is updated by the Eq (2.9).

where ${{\tilde x}_i} = {{\hat x}_i} - {x_i}$, and ${{\tilde W}_i} = {{\hat W}_i} - W_i^*$ is the constant weight vector.

With Eqs (2.4), (2.8) and (2.9), the derivative of the estimation error ${\tilde x}_i$ can be calculated by Eq (2.10).

The overall system is shown in Eqs (2.11) and (2.12).

and

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 RBFNN $\bar WS(x)$. The $M$ estimators of training gait features are constructed by using the learning knowledge acquired during the training stage, as shown below:

where $j = 1, \ldots, M$, $x_j = {{\rm{[}}{{x}_{j1}}, \ldots, {{x}_{jn}}{\rm{]}}^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 = {{\rm{[}}{{e}_{j1}}, \ldots, {{e}_{jn}}{\rm{]}}^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_{tj}}$ 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 output $\dot 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 output $\tilde{\dot x}$ 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.

3.   Experiments and results

3.1.   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 \times 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.

Table 1.  Demographics of the PD patients and healthy controls in the online dataset ^[28] and our dataset (Values are expressed as mean $\pm$ standard deviation).

Clinical characteristics	HC (Online [28])	PD (Online [28])	PD (Ours)
Age /years	63.65 $\pm$ 8.70	66.23 $\pm$ 9.59	63.84 $\pm$ 9.18
Weight /kg	72.76 $\pm$ 12.35	72.39 $\pm$ 11.96	68.18 $\pm$ 13.15
Height /m	1.68 $\pm$ 0.09	1.67 $\pm$ 0.08	1.68 $\pm$ 0.09
Gender	40 male, 31 female	57 male, 34 female	33 male, 17 female
UPDRS	-	31.62 $\pm$ 11.83	36.32 $\pm$ 16.29
Severity	-	55 mild, 36 moderate	23 mild, 27 moderate

---

 | Show Table

DownLoad: CSV

3.2.   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 $\dot x$. The size of $x$ and $\dot x$ 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 ($Sen$), specificity ($Spe$), 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 sub-regions 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.

Table 2.  Performance of the proposed method in PD diagnosis.

Evaluation metrics	Total +RBFNN	Total +DFE+RBFNN	Sub-region +RBFNN	Sub-region +DFE+RBFNN (Proposed)
Accuracy	87.57%	93.89%	95.63%	100.00%
Precision	87.91%	95.56%	95.07%	100.00%
Sensitivity	91.11%	93.56%	97.78%	100.00%
Specificity	83.04%	94.29%	92.86%	100.00%
F1 score	89.13%	94.44%	96.25%	100.00%

---

 | Show Table

DownLoad: CSV

Table 3.  Performance of the proposed method in severity assessment of PD.

Evaluation metrics	Total +RBFNN	Total +DFE+RBFNN	Sub-region +RBFNN	Sub-region +DFE+RBFNN (Proposed)
Accuracy	63.78%	69.89%	79.03%	95.89%
Precision	72.48%	74.39%	82.79%	98.33%
Sensitivity	66.00%	78.00%	83.67%	95.00%
Specificity	60.00%	57.50%	72.50%	97.50%
F1 score	68.09%	75.87%	82.61%	96.18%

---

 | Show Table

DownLoad: CSV

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.

Table 4.  PD diagnosis classification results in the online dataset of comparisons with existing studies.

Reference	Accuracy	Sensitivity	Specificity
Maachi.(2020) [30]	98.70%	98.10%	100.00%
Veeraragavan.(2020) [16]	97.40%	97.70%	97.05%
Yurdakul.(2020) [15]	98.30%	98.60%	97.50%
Setiawan.(2021) [31]	96.63%	94.46%	97.69%
Ferreira.(2022) [32]	84.60%	80.00%	-
Proposed method	100.00%	100.00%	100.00%

---

 | Show Table

DownLoad: CSV

Table 5.  The results of PD severity assessment in the online dataset of comparisons with existing studies.

Reference	Accuracy	Sensitivity	Specificity
Maachi.(2020) [30]	85.30%	85.30%	87.30%
Veeraragavan.(2020) [16]	87.10%	-	-
Setiawan.(2021) [31]	94.58%	92.08%	95.60%
Ferreira.(2022) [32]	78.60%	72.30%	72.30%
Proposed method	95.89%	95.00%	97.50%

---

 | Show Table

DownLoad: CSV

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.

Table 6.  Performance of the severity assessment by the proposed method in our dataset.

Evaluation metrics	Total +RBFNN	Total +DFE+RBFNN	Sub-region +RBFNN	Sub-region +DFE+RBFNN (Proposed)
Accuracy	66.00%	68.00%	80.00%	96.00%
Precision	67.00%	63.33%	78.33%	96.67%
Sensitivity	76.67%	86.67%	88.33%	96.67%
Specificity	50.00%	55.00%	76.67%	96.67%
F1 score	67.38%	71.33%	80.33%	96.00%

---

 | Show Table

DownLoad: CSV

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.

Table 7.  Validation of the generalization ability.

Evaluation Metrics	Training set: online dataset. Test set: our dataset.	Training set: our dataset. Test set: online dataset.
Accuracy	87.91%	90.00%
Precision	92.31%	84.62%
Sensitivity	87.27%	95.65%
Specificity	88.89%	85.19%
F1 score	89.72%	89.80%

---

 | Show Table

DownLoad: CSV

4.   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 cross-validation 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.

5.   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.

Acknowledgments

This research was funded by the National Natural Science Foundation of China (U1913208), the Science and Technology Program of Tianjin (21JCZDJC00170), Tianjin Key Medical Discipline (Specialty) Construction Project (TJYXZDXK-052B), Science and Technology Program of Tianjin (21JCYBJC01440) and Tianjin Health Research Project (TJWJ2022MS033, TJWJ2022MS031).

Conflict of interest

The authors declare there is no conflict of interest.