Identification and Monitoring of Parkinson’s Disease Dysgraphia Based on Fractional-Order Derivatives of Online Handwriting ^†

Abstract

Parkinson’s disease dysgraphia affects the majority of Parkinson’s disease (PD) patients and is the result of handwriting abnormalities mainly caused by motor dysfunctions. Several effective approaches to quantitative PD dysgraphia analysis, such as online handwriting processing, have been utilized. In this study, we aim to deeply explore the impact of advanced online handwriting parameterization based on fractional-order derivatives (FD) on the PD dysgraphia diagnosis and its monitoring. For this purpose, we used 33 PD patients and 36 healthy controls from the PaHaW (PD handwriting database). Partial correlation analysis (Spearman’s and Pearson’s) was performed to investigate the relationship between the newly designed features and patients’ clinical data. Next, the discrimination power of the FD features was evaluated by a binary classification analysis. Finally, regression models were trained to explore the new features’ ability to assess the progress and severity of PD. These results were compared to a baseline, which is based on conventional online handwriting features. In comparison with the conventional parameters, the FD handwriting features correlated more significantly with the patients’ clinical characteristics and provided a more accurate assessment of PD severity (error around 12%). On the other hand, the highest classification accuracy (ACC = 97.14%) was obtained by the conventional parameters. The results of this study suggest that utilization of FD in combination with properly selected tasks (continuous and/or repetitive, such as the Archimedean spiral) could improve computerized PD severity assessment.

1. Introduction

As a second most common neurodegenerative disorder, Parkinson’s disease (PD) is expected to impose an increasing social and economic burden on societies as populations age [1]. Its prevalence rate is estimated to approximately 1.5% for people aged over 65 years [2]. The risk of being affected by PD strongly increases with age, and, in the next 15 years, the incidence of PD is expected to be doubled [3,4]. The rapid degeneration of dopaminergic cells in the substantia nigra pars compacta [5] arose as the most significant biological finding associated with the disease, but the exact pathophysiological cause of PD has not yet been discovered. PD cardinal motor symptoms involve bradykinesia (slowness of movement), tremor at rest, rigidity, gait impairment, and postural instability [6,7,8]. A variety of non-motor symptoms may emerge as well—for instance, cognitive impairment, dementia, depression, sleep disorders, or anxiety [6,9,10].

Handwriting requires cognitive, perceptual, and fine motor abilities. In conjunction with motor dysfunctions in people suffering from PD, it has been proven that disrupted handwriting may be used as a significant biomarker for PD diagnosis [11,12]. Micrographia, which is associated with the progressive decrease in letters’ amplitude, is the most commonly observed handwriting abnormality in patients with PD [13,14]. Moreover, according to McLennan et al. [14], in approximately 5% of PD patients, micrographia may be observed even before the onset of the cardinal motor symptoms.

The recent advantage of new technologies coming hand-in-hand with Health 4.0 systems enables the acquisition of online handwriting signals, where temporal information is added to the x and y position. Therefore, by using a digitizing tablet, the analysis is not limited to spatial features which mainly quantify PD micrographia. In addition, we are able to quantify temporal, kinematic, and dynamic manifestations of PD dysgraphia, such as hesitations, pauses, and slow movement [7], which cannot be studied objectively using a classical paper-and-pen method. Due to this complexity, Letanneux et al. [15] started to refer to these manifestations using the generalized term PD dysgraphia.

Several research teams have explored the impact of quantitative PD dysgraphia analysis utilizing simple handwriting/drawing tasks (e.g., separate characters, a combination of two or three characters, repetitive loops, circles), as well as more complex ones (e.g., words, sentences, figures, 3D objects, and the Archimedean spiral) [8,16,17,18,19,20]. An overview of recent related works (2015–present) can be seen in

Table 1. Overview of related works focused on computerized analysis of Parkinson’s disease (PD) dysgraphia.

First Author	Year	PD/HC	Handwriting Task	Analysis	Features	Conclusions
Drotar * [17]	2015	37/38	letters, words, sentences	differential analysis (SVM)	kinematic, temporal, spatial, entropy, EMD, signal energy	The highest classification accuracy after feature selection approach was 88.13%.
Drotar * [16]	2015	37/38	letters, words, sentences	differential analysis (SVM)	kinematic, temporal, spatial, entropy, EMD, pressure	Classification performance was at its peak with on-surface features equal to AUC = 89.09%.
Heremans [24]	2015	34/10	up/down strokes at varying amplitudes	ANOVA	spatial and kinematic	Significant difference between groups was in spatial (F(2.41) = 3.97; p = 0.03).
Pereira [25]	2015	37/18	Archimedean spiral	differential an. (SVM, NB, OPF)	mean relative tremor and spatial parameters	The best results were obtained by NB classifier that provided around 79% classification accuracy.
Drotar * [8]	2016	37/38	letters, words, Archimedean spiral, sentences	differential an. (SVM, K-NN, ADA)	kinematic, temporal, spatial, entropy, EMD, pressure	Combining all exercises, SVM proved to be the best classifier with 82.5% accuracy.
Heremans [26]	2016	30/15	repetitive cursive loops	ANOVA, correlation an.	writing amplitude and velocity	PD dysgraphia is more severe in patients with freezing of gait.
Pereira [27]	2016	14/21	Archimedean spiral, meander	differential an. (CNN, OPF)	pen-based features	The best result was obtained by CNN with 87.14% classification accuracy using meander task.
Kotsavasil [21]	2017	24/20	horizontal lines	differential analysis (NB)	kinematic	Average classification accuracy was 91%.
Loconsole [18]	2017	4/7	sentence, repetitive loops	differential analysis (ANN)	temporal, kinematic, spatial	Highest classification accuracy (96.81%) was achieved using all the extracted features.
Taleb [22]	2017	16/16	letters, waves, words	differential analysis (SVM)	kinematic, stroke, pressure, entropy, energy, EMD	The highest classification accuracy was 96.88% for 12 kinematic and pressure features.
Moetesum * [23]	2018	37/38	Archimedean spiral, letters, words, sentence, loops	differential analysis (SVM)	CNN-based features	Extraction of features using CNN applied on raw handwriting data resulted in 83% classification accuracy.
Mucha * [20]	2018	30/36	Archimedean spiral	differential analysis (RF, SVM)	fractional derivatives based kinematic features	Improvement of classification accuracy by 10% (72.38%) in comparison to the baseline.
Impedovo * [28]	2018	37/38	Archimedean spiral letters, words, sentence	differential an. (RF, SVM, K-NN, NB, LDA, ADA)	kinematic, temporal, spatial, entropy, EMD, pressure	Analysis focused on PD diagnosis at earlier stages resulted in 74.76% classification accuracy.

SVM—support vector machine; EMD—empirical mode decomposition; K-NN—K-nearest neighbors; ANOVA—analysis of variance; NB—naïve Bayes classifier; OPF—optimum path forest; ANN—artificial neural network; CNN—convolutional neural network; RF—random forests; LDA—Linear Discriminant Analysis; ADA—AdaBoost; AUC—area under the receiver operating characteristics (ROC) curve; articles are sorted by the year of release and then alphabetically; * analyzes performed on the same database (Parkinson’s disease handwriting database (PaHaW) [8]).

. Most of them confirm the irreplaceability of kinematic features in PD dysgraphia analysis. Additionally, the researchers usually employ temporal, spatial, and dynamic features. Some more advanced parameters are reported too. For instance, Drotar et al. [8,16,17] demonstrated a combination of kinematic, pressure, energy, or empirical mode decomposition (EMD)-based features that resulted in a classification accuracy of up to 89% using several handwriting tasks. Kotsavasilogloua et al. [21] achieved an average prediction accuracy of 91% using simple horizontal lines and features describing the variability in the pen tip’s velocity, a deviation from the horizontal plane, and the trajectory’s entropy. Other works report even higher classification accuracies (approximately 97%), e.g., Loconsole et al. [18], who used computer vision and electromyography signal processing techniques, or Taleb et al. [22], who used a combination of features related to the correlation between kinematic and pressure characteristics (but, in this case, applied to a very small dataset). Another promising approach was published by Moetesum et al. [23], who reached an 83% classification accuracy by employing convolutional neural networks (CNN) that were used to extract discriminating visual features from handwriting data transformed into the offline mode. In 2018, Impedovo et al. reported the results of a study focused only on the early stages of PD; the best accuracy was 74.76% for a combination of three handwriting tasks. Finally, in our previous work [20], we proposed a new approach of advanced kinematic feature extraction that utilizes fractional-order derivatives (FD). This approach increased the classification accuracy by 10% (72.39%) for Archimedean spiral tasks in comparison with the baseline [20].

Although the authors of the previously mentioned studies reported high classification accuracies, further signal processing and machine learning pipeline improvements are expected to make the differential analysis even more accurate. One possible approach could involve an advanced feature extraction methodology based on fractional calculus (FC) [29,30], which enables the use of an arbitrary order of derivatives and/or integrals. Generally, FC has many applications in different fields of science [31,32,33]. For instance, it has been advantageously used during the modeling of different diseases, such as human immunodeficiency virus (HIV) [34] and malaria [35]. In addition, FC-based analytical tools have outperformed classical techniques in geology [36,37], economics and finance [38,39], etc. Moreover, in our recent paper [20], we identified a high potential for the use of FC in the kinematic analysis of PD drawings. Based on these preliminary results, we assume that FD-based handwriting features may bring improvements to PD diagnosis and assessment. In the frame of this article, we would like to go further and deeply explore the impact of FD on the PD dysgraphia diagnosis and its monitoring. More specifically, we aim to:

• investigate the relationship between newly designed FD handwriting features and a patient’s clinical data and compare these results with a baseline (i.e., results based on conventional parameters),
• evaluate the discrimination power of the FD features in terms of binary classification accuracy and compare the results to the baseline,
• use the newly designed features to establish regression models that will estimate the severity of PD and compare its performance to that of a baseline.

The rest of this paper is organized as follows: Section 2 describes the cohort of patients and the methodology, and Section 3 includes the results. A discussion is presented in Section 4, and, finally, conclusions are drawn in Section 5.

2. Materials and Methods

2.1. Dataset

For the purpose of this work, the Parkinson’s disease handwriting database (PaHaW) [8], which consists of multiple handwriting/drawing samples from 37 PD patients and 38 age- and gender-matched healthy controls (HC), was used. Since the Archimedean spiral drawing task is missing for some participants, we reduced the analyzed cohort to 33 PD patients and 36 HC. Demographic and clinical data of the participants can be found in

Table 2. Demographic and clinical data of the enrolled participants.

Gender	N	Age [years]	PD dur [years]	UPDRS V	LED [mg/day]
Parkinson’s disease patients
Females	17	71.76 ± 10.93	9.88 ± 5.27	2.18 ± 0.86	1146.03 ± 543.89
Males	16	66.50 ± 13.44	7.44 ± 4.04	2.31 ± 0.75	1673.38 ± 616.66
All	33	69.21 ± 11.10	8.70 ± 4.82	2.24 ± 0.80	1401.72 ± 630.71
Healthy controls
Females	17	61.59 ± 10.17	-	-	-
Males	19	63.32 ± 13.14	-	-	-
All	36	62.50 ± 11.70	-	-	-

PD—Parkinson’s disease; N—number of subjects; PD dur—PD duration; UPDRS V—Unified Parkinson’s disease rating scale, part V: Modified Hoehn and Yahr staging score [40]; LED—L-dopa equivalent daily dose [41].

. The participants were enrolled at the First Department of Neurology, St. Anne’s University Hospital in Brno, Czech Republic. All participants reported the Czech language as their native language and were right-handed. The patients completed their tasks approximately 1 h after their regular dopaminergic medication (L-dopa). All participants signed an informed consent form approved by the local ethics committee. Unified Parkinson’s disease rating scale, part V (UPDRS V): Modified Hoehn and Yahr staging score [40], was used to assess clinical symptoms of PD. In the frame of this work, the duration of the disease was considered as well. Descriptive visualization (histograms, regression, and residual plots) of the clinical data for the subjects participating in this study can be seen in Figure 1.

2.2. Data Acquisition

The PaHaW database [8] includes nine different handwriting tasks written in the Czech language. Their description and translation to English can be found in

Table 3. Description of the PaHaW handwriting tasks.

N	Task	Czech (Original)	English (Translation)
1	Archimedean spiral	-	-
2	repetitive loops	-	-
3	letter	l	l
4	syllable	le	le
5	word	les	forest
6	word	lektorka	lecturer
7	word	porovnat	compare
8	word	nepopadnout	not grasped
9	sentence	Tramvaj dnes už nepojede.	The tram will no longer go today.

. During all handwriting tasks, the participants were rested and seated in a comfortable position with the possibility to look at the prefilled template (see Figure 2). A digitizing tablet (Wacom Intuos 4M, Wacom, Kazo, Saitama, Japan) was overlaid with an empty paper template and participants were asked to perform all tasks using a special Wacom inking pen that gave the patients immediate visual feedback. Online handwriting signals were recorded with a sampling frequency of $f_{\mathrm{s}}=150\mathrm{Hz}$. The following time sequences were acquired: x and y coordinates ($x[t]$, $y[t]$); time-stamp (t); in-air/on-surface (on-surface movement is a movement of a pen when its tip is touching the surface, e.g., paper (i.e., it provides the information about the pen writing/drawing on the paper); vice versa, in-air movement is a movement of a pen when its tip is up to 1.5 cm above the surface [42,43]) status ($b[t]$); pressure ($p[t]$); azimuth ($az[t]$); and altitude ($al[t]$).

2.3. Feature Extraction

The main goal of this work is to compare a set of commonly used kinematic features with newly proposed FD-based features in terms of quantitative PD dysgraphia analysis. All of the handwriting features were computed using both on-surface as well as in-air movements. The two movements were quantified separately using velocity (rate at which the position of the pen changes with time [mm/s]), acceleration (rate at which the velocity of the pen changes with time [mm/s²]), jerk (rate at which the acceleration of the pen changes with time [mm/s³]), and their horizontal and vertical variants [8,44,45]. FD-based features were extracted for different values of $\alpha$. In the frame of this work, $\alpha$ ranging from 0.1 to 1.0 with a step of 0.1 was used. Subsequently, the statistical properties of the computed handwriting features were described using the mean, median, standard deviation (std), and maximum (max). Finally, all of the extracted features were divided into nine different feature sets according to the type of the movement (on-surface, in-air, and combined) and the calculation approach, i.e., the type of feature (FD-based, conventional, and combined). For more information, see

Table 4. Feature sets matrix.

Movement	FD (Count)	Conventional (Count)	Together (Count)
on-surface	4536	618	5154
in-air	2916	404	3320
together	7452	1022	8474

.

Fractional-Order Derivatives

Utilization of the FD as a substitution for the conventional differential derivative during calculation of the basic kinematic features provides a new advanced approach. The advantage of FDs is in their wide range of settings and many different approaches to approximation, e.g., Riemann–Liouville, Caputo, or Grünwald–Letnikov formulations [31,46,47]. For the purpose of this work, Jonathan Hadida’s FD Matlab implementation was used following the Grünwald–Letnikov approximation [31,48]. A direct definition of the FD $D^{\alpha }y\left(t\right)$ is based on the finite differences of an equidistant grid in $[0,\tau ]$, assuming that the function $y\left(\tau \right)$ satisfies certain smoothness conditions in every finite interval $(0,t),t\le T$. Choosing the grid [31],

$$0={\tau }_0<{\tau }_1<\dots <{\tau }_{n+1}=t=(n+1)h$$

(1)

with

$${\tau }_{k+1}-{\tau }_k=h$$

(2)

and using the notation of finite differences

$$\frac{1}{h^{\alpha }}{\Delta }_h^{\alpha }y\left(t\right)=\frac{1}{h^{\alpha }}\left(y\left({\tau }_{n+1}\right)-\sum _{v=1}^{n+1}{c}_v^{\alpha }y\left({\tau }_{n+1-v}\right)\right),$$

(3)

where

$$c_v^{\alpha }={(-1)}^{v-1}{(}_v^{\alpha }).$$

(4)

The Grünwald–Letnikov implementation is defined as

$$D^{\alpha }y\left(t\right)=\underset{h\to 0}{lim}\frac{1}{h^{\alpha }}{\Delta }_h^{\alpha }y\left(t\right),$$

(5)

where $D^{\alpha }y\left(t\right)$ denotes a derivative with order $\alpha$ of function $y\left(t\right)$, and h represents a sampling lattice.

2.4. Statistical Analysis

Prior to providing a description of the analytical setup, it is important to note that the effect of well-known confounding factors, also known as covariates, was controlled for in all of the analytical steps described below. In the frame of this work, we controlled for the effect of participants’ age, gender, and L-dopa [41] (dopaminergic medication).

To assess the strength of the relationship between the computed handwriting features and patient’s clinical data (UPDRS V and PD duration), we computed the partial Pearson’s correlation coefficient (assessment of a linear relationship), as well as the partial Spearman’s correlation coefficient (assessment of a monotonic relationship). With this approach, we aimed to identify the handwriting features that are significantly correlated with the clinical measures under focus and also to compare the FD features with conventional ones. A significance level of correlation (p) of 0.05 was selected for both of the correlation types. Only the results with a p-value below the significance level in both correlation coefficients were considered statistically significant.

Next, to evaluate and compare the power of the handwriting features to discriminate PD patients and HC, multivariate binary classification analysis was performed. For this purpose, state-of-the-art gradient boosted trees were employed. Specifically, we used the famous XGBoost algorithm [49]. The XGBoost algorithm was chosen for its ability to achieve a good performance, even for small datasets; its inherent robustness to outliers; its ability to model complex interdependencies in the data; and also its recent successes in the field of machine learning (e.g., the winning algorithm in many www.kaggle.com competitions). To train and evaluate the models, we used the following supervised learning setup: stratified 10-fold cross-validation with 20 repetitions. The performance of the trained classification models was evaluated by Matthew’s correlation coefficient (MCC) [50], classification accuracy (ACC), sensitivity (SEN), and specificity (SPE), which are defined as follows:

$$\mathrm{MCC}=\frac{\mathrm{TP}\times \mathrm{TN}-\mathrm{FP}\times \mathrm{FN}}{\sqrt{(\mathrm{TP}+\mathrm{FP})(\mathrm{TP}+\mathrm{FN})(\mathrm{TN}+\mathrm{FP})(\mathrm{TN}+\mathrm{FN})}},$$

(6)

$$\mathrm{ACC}=\frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN}}·100[\%],$$

(7)

$$\mathrm{SEN}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}·100[\%],$$

(8)

$$\mathrm{SPE}=\frac{\mathrm{TN}}{\mathrm{TN}+\mathrm{FP}}·100[\%],$$

(9)

where TP is the number of true positives, TN is the number of true negatives, FP is the number of false positives, and FN is the number false negatives.

Finally, to evaluate and compare the power of the handwriting features’ ability to predict the values of the selected clinical characteristics (UPDRS V and PD duration), multivariate regression analysis was performed. For this purpose, the same boosting tree algorithm (XGBoost) and the supervised learning setup were used. The performance of the trained regression models was evaluated by the mean absolute error (MAE), root mean square error (RMSE), and estimated error rate (EER), which are defined as follows:

$$\mathrm{MAE}=\frac{1}{n}\sum _{i=1}^n\left|y_i-\widehat{y_i}\right|,$$

(10)

$$\mathrm{RMSE}=\sqrt{\frac{1}{n}\sum _{i=1}^n{(y_i-\widehat{y_i})}^2},$$

(11)

$$\mathrm{EER}=\frac{1}{n·r}\sum _{i=1}^n\left|y_i-\widehat{y_i}\right|·100[\%],$$

(12)

where $y_i$ represents the true label of the ith observation, $\widehat{y_i}$ denotes the predicted label of the ith observation, n is the number of observations, and r is the range of the values of the predicted clinical characteristic (not the range that can be theoretically reached, but the actual range of the values in the dataset). Therefore, the EER describes a percentage of error predictions in regard to the statistical properties of the data.

3. Results

In

Table 5. Results of partial correlation analysis between handwriting features and clinical data.

UPDRS V
FD on-suface					Conventional on-surface
feature name	$\alpha$	task	$r_p$	$r_s$	$r_s$	$r_p$	task	feature name
jerk (median)	0.3	r. letters l	0.37 *	0.48 **	−0.45 *	−0.24	r. letters le	h. jerk (max)
jerk (median)	0.4	r. letters l	0.43 *	0.46 *	−0.43 *	−0.2	r. letters le	velocity (max)
h. velocity (std)	0.1	r. letters l	−0.42 *	−0.41 *	−0.42 *	0.25	r. letters l	h. jerk (max)
v. acceleration (mean)	0.7	r. loops	0.48 **	0.40 *	−0.42 *	−0.16	r. letters l	h. velocity (max)
v. velocity (std)	0.3	sentence	0.40 *	0.40 *	−0.41 *	−0.15	letter l	h. velocity (max)
FD in-air					Conventional in-air
feature name	$\alpha$	task	$r_p$	$r_s$	$r_s$	$r_p$	task	feature name
v. velocity (median)	0.9	sentence	0.44 *	0.53 **	0.43 *	0.28	r. word lektorka	acceleration (mean)
v. velocity (median)	0.8	sentence	0.40 *	0.52 **	−0.37 *	−0.31	word porovnat	h. jerk (max)
h. velocity (median)	0.5	r. letters le	−0.38 *	−0.49 **	0.36 *	0.25	r. letters l	v. velocity (median)
v. jerk (median)	0.3	r. letters le	−0.43 *	−0.49 **	0.35	0.41 *	r. letters le	h. velocity (median)
v. velocity (median)	0.7	sentence	0.37 *	0.48 **	0.35	0.19	r. word lektorka	acceleration (median)
PD Duration
FD on-surface					Conventional on-surface
feature name	$\alpha$	task	$r_p$	$r_s$	$r_s$	$r_p$	task	feature name
velocity (max)	0.1	spiral	0.54 **	0.55 **	−0.46 *	−0.40 *	r. letters l	h. velocity (max)
acceleration (max)	0.8	spiral	0.54 **	0.54 **	−0.40 *	−0.37 *	r. letters l	h. jerk (max)
acceleration (max)	0.6	spiral	0.54 **	0.54 **	−0.38 *	−0.37 *	r. letters l	velocity (max)
acceleration (max)	0.2	spiral	0.54 **	0.54 **	0.46 **	0.34	spiral	v. velocity (mean)
acceleration (max)	0.7	spiral	0.54 **	0.53 **	0.40 *	0.14	r. loops	h. acceleration (mean)
FD in-air					Conventional in-air
feature name	$\alpha$	task	$r_p$	$r_s$	$r_s$	$r_p$	task	feature name
jerk (median)	0.4	sentence	−0.37 *	−0.49 **	−0.44 *	−0.38 *	word lektorka	h. jerk (median)
jerk (max)	0.1	r. word les	0.57 **	0.46 *	0.38 *	0.40 *	word nepopad.	velocity (max)
jerk (max)	0.3	r. word les	0.57 **	0.45 *	0.37 *	0.42 *	word lektorka	h. n. jerk (mean)
velocity (max)	0.1	r. word les	0.57 **	0.45 *	−0.47 **	−0.13	r. word lektorka	h. velocity (mean)
jerk (max)	0.2	r. word les	0.57 **	0.45 *	−0.42 *	−0.13	word nepopad.	h. velocity (mean)

$\alpha$—order of FD; $r_p$—Pearson’s correlation coefficient; $r_s$—Spearman’s correlation coefficient; v.—vertical; h.—horizontal; r.—repetitive task; *—$p<0.05$; **—$p<0.01$; rows are ordered by the absolute value of Spearman’s correlation coefficient.

, the results of partial correlation analysis between the handwriting features (FD-based features, conventional features) and patients’ clinical characteristics (UPDRS V, PD duration) are summarized. The table shows the five best features according to Spearman’s correlation coefficient for each movement (on-surface, in-air).

In the case of UPDRS V (on-surface movement), the following FD-based features achieved a statistical significance of correlation: the median of jerk ($\alpha =0.3$, $\alpha =0.4$) and horizontal velocity ($\alpha$ = 0.1) for the repetitive letter l, the mean of vertical acceleration ($\alpha =0.7$) for repetitive loops, and the standard deviation of the vertical velocity ($\alpha =0.3$) for the sentence. The following conventional features achieved a statistical significance of correlation (p-value of only one of the coefficients was below the threshold): the maximum of horizontal jerk and velocity for the repetitive letters le, the maximum of horizontal jerk and horizontal velocity for the repetitive letter l, and the maximum of horizontal velocity for the letter l. Regarding UPDRS V (in-air movement), the following FD-based features achieved a statistical significance of correlation: the median of vertical velocity ($\alpha =0.9$, $\alpha =0.8$, $\alpha =0.7$) for the sentence and the median of horizontal velocity ($\alpha =0.5$) and vertical jerk ($\alpha =0.3$) for the repetitive letters le. The following conventional features achieved a statistical significance of correlation (p-value of only one of the coefficients was below the threshold): the mean of acceleration for the repetitive word lektorka, the maximum of horizontal jerk for the word porovnat, the median of the vertical velocity for the repetitive letter l, and the median of the horizontal velocity of the repetitive letters le.

For PD duration (on-surface movement), the following FD-based features achieved a statistical significance of correlation (of note: all of these features satisfied the stronger threshold for statistical significance of correlation $p<0.01$): the maximum of the velocity ($\alpha =0.1$) and acceleration ($\alpha =0.8$, $\alpha =0.7$, $\alpha =0.6$, $\alpha =0.2$) for the Archimedean spiral. The following conventional features achieved a statistical significance of correlation (p-value of only one of the coefficients was below the threshold): the maximum of horizontal velocity, horizontal jerk, and velocity for the repetitive letter l; the mean of the vertical velocity for the Archimedean spiral; and the mean of horizontal acceleration for repetitive loops. For PD duration (in-air movement), the following FD-based features achieved a statistical significance of correlation: the median of jerk ($\alpha =0.4$) for sentence, the maximum of jerk ($\alpha =0.1$, $\alpha =0.2$, $\alpha =0.3$) and velocity ($\alpha =0.1$) for repetitive word les. The following conventional features achieved a statistical significance of correlation (p-value of only one of the coefficients was below the threshold): the median and mean of horizontal jerk for the word lektorka, the maximum of the velocity for the word nepopadnout, and the mean of horizontal velocity for the repetitive word lektorka and the word nepopadnout.

The results of the multivariate binary classification analysis are summarized in

Table 6. Results of multivariate binary classification analysis (PD/HC).

Feature Set	MCC	ACC [%]	SEN [%]	SPE [%]	Feat
conventional on-surface	0.83 ± 0.18	91.19 ± 9.65	93.00 ± 15.52	70.00 ± 0.46	1
conventional in-air	0.95 ± 0.10	97.14 ± 5.71	95.50 ± 9.07	100.00 ± 0.00	1
conventional together	0.95 ± 0.11	97.14 ± 5.71	95.50 ± 9.07	100.00 ± 0.00	1
FD on-surface	0.95 ± 0.12	87.14 ± 13.48	82.00 ± 21.24	90.00 ± 30.00	1
FD in-air	0.95 ± 0.13	81.43 ± 12.86	71.50 ± 30.83	60.00 ± 48.99	3
FD together	0.95 ± 0.14	81.43 ± 15.71	69.50 ± 32.13	70.00 ± 45.83	2
all on-surface	0.95 ± 0.15	88.33 ± 14.06	89.00 ± 22.11	70.00 ± 45.83	2
all in-air	0.95 ± 0.16	97.14 ± 5.71	95.50 ± 9.07	100.00 ± 0.00	1
all together	0.95 ± 0.17	97.14 ± 5.71	95.50 ± 9.07	100.00 ± 0.00	1

MCC—Matthew’s correlation coefficient; ACC—accuracy; SEN—sensitivity; SPE—specificity; feat.—number of features important for the trained model (i.e., feature importance of the feature > 0.0); The feature importances, as well as the exact names of these features, are summarized in the text.

. In total, we built and evaluated nine different classification models. These models were selected according to the following criteria: movement type (on-surface, in-air, all), feature type (FD features, conventional features, all). We built models based on the combinations of these criteria as well. For more information, see Table 4.

With respect to the classification performance, the highest MCC achieved was 0.95 was for eight out of the total nine feature sets (with the exception being the feature set composed of conventional handwriting features computed for the on-surface movements). An interesting fact to note is that for all models based on conventional handwriting features, only a single feature was capable of providing the classification models with such a high discrimination power. In terms of the specific features important for the trained models, the following feature importances were returned by the models (feature importance quantifies the relative importance of the features in the ensemble of the trained XGBoost model [49]; therefore, the higher the value of the feature importance, the more important the feature for the prediction of the dependent variable): conventional on-surface (horizontal jerk (median) of repetitive loops), conventional in-air (horizontal velocity (median) of the sentence), conventional together (horizontal velocity (median) of the sentence), FD on-surface (jerk (max) $\alpha =0.3$ of the letters le), FD in-air (vertical acceleration (mean) $\alpha =0.6$ of the word nepopadnout ($\mathrm{FI}=0.33$), horizontal jerk (mean) $\alpha =0.9$ of the word nepopadnout ($\mathrm{FI}=0.33$), horizontal jerk (mean) $\alpha =0.2$ of the repetitive word lektorka ($\mathrm{FI}=0.33$)), FD together (jerk (max) $\alpha =0.3$ of the letters le (on-surface; $\mathrm{FI}=0.67$), horizontal jerk (mean) $\alpha =0.9$ of the word nepopadnout (in-air; $\mathrm{FI}=0.33$)), all on-surface (horizontal jerk (median) of repetitive loops ($\mathrm{FI}=0.50$), jerk (max) $\alpha =0.3$ of the letters le ($\mathrm{FI}=0.50$)), all in-air (horizontal velocity (median) of the sentence), and all together (horizontal velocity (median) of the sentence (in-air)).

The results of multivariate regression analysis are summarized in

Table 7. Results of regression analysis for clinical data.

Feature Set	MAE	RMSE	EER [%]	Feat
UPDRS V
conventional on-surface	0.59 ± 0.29	0.71 ± 0.41	13.82 ± 6.71	1
conventional in-air	0.60 ± 0.30	0.72 ± 0.42	14.01 ± 6.98	1
conventional together	0.60 ± 0.31	0.73 ± 0.42	14.05 ± 6.90	1
FD on-surface	0.60 ± 0.32	0.65 ± 0.45	12.51 ± 7.55	1
FD in-air	0.60 ± 0.33	0.68 ± 0.43	13.49 ± 7.29	1
FD together	0.60 ± 0.34	0.66 ± 0.45	13.06 ± 7.55	2
all on-surface	0.60 ± 0.35	0.65 ± 0.45	12.51 ± 7.55	1
all in-air	0.60 ± 0.36	0.71 ± 0.43	13.72 ± 7.36	1
all together	0.60 ± 0.37	0.66 ± 0.45	13.06 ± 7.55	2
PD duration
conventional on-surface	4.29 ± 0.94	5.03 ± 1.09	24.52 ± 5.39	18
conventional in-air	4.91 ± 1.38	5.56 ± 1.50	28.03 ± 7.85	16
conventional together	4.14 ± 1.32	4.85 ± 1.52	23.64 ± 7.55	16
FD on-surface	4.45 ± 0.66	5.06 ± 0.85	25.40 ± 3.75	14
FD in-air	4.79 ± 0.73	5.48 ± 0.72	27.36 ± 4.20	19
FD together	4.55 ± 0.68	5.32 ± 0.78	26.00 ± 3.88	21
all on-surface	4.48 ± 0.86	5.12 ± 0.96	25.62 ± 4.92	16 (12 F, 4 C)
all in-air	4.95 ± 1.18	5.59 ± 1.17	28.30 ± 6.75	17 (13 F, 4 C)
all together	4.70 ± 1.10	5.45 ± 1.23	26.82 ± 6.30	17 (12 F, 6 C)

UPDRS V—Unified Parkinson’s disease rating scale, part V: Modified Hoehn and Yahr staging score [40]; MAE—mean absolute error; RMSE—root mean squared error; EER—estimation error rate; F—FD-based features; C—conventional handwriting features; feat.—number of features important for the trained model (i.e., feature importance of the feature > 0.0); The feature importances, as well as the exact names of these features for models built to assess UPDRS V, are summarized in the text. In the case of PD duration, this data can be found in Table S1 provided in the Supplementary Material.

. For this purpose, we used UPDRS V and PD duration as our target variables. As in the case of binary classification, we built and evaluated nine different regression models according to the same criteria. For each of the rating scales, the table shows the results achieved using the trained models and the associated feature importance values. All obtained results are discussed in the following section.

Considering EER as our performance evaluation metric, the following results are worth pointing out. In the case of UPDRS V, the lowest EER was achieved using a single FD-based feature—specifically, the standard deviation of vertical velocity ($\alpha =0.1$) computed for the on-surface movements ($12.51\pm 7.55\%$). The same feature was selected when both FD and conventional features were considered while building the model. In general, all models achieved an EER of around 12–13%. In comparison with the conventional features, the FD-based features performed better, with a difference of about 1%. In terms of the specific features important for the trained models, the following feature importances were returned by the models: conventional on-surface (vertical normalized jerk (mean) of the repetitive word lektorka), conventional in-air (vertical velocity (mean) of the sentence), conventional together (vertical velocity (mean) of the sentence), FD on-surface (vertical velocity (std) $\alpha =0.1$ of the sentence), FD in-air (vertical velocity (median) $\alpha =0.3$ of the sentence), FD together (vertical velocity (std) $\alpha =0.1$ of the sentence (on-surface; $\mathrm{FI}=0.50$), vertical velocity (median) $\alpha =0.3$ of the sentence (in-air; $\mathrm{FI}=0.50$)), all on-surface (vertical velocity (std) $\alpha =0.1$ of the sentence), all in-air (vertical velocity (median) $\alpha =0.3$ of the sentence), and all together (vertical velocity (std) $\alpha =0.1$ of the sentence (on-surface; $\mathrm{FI}=0.50$), vertical velocity (median) $\alpha =0.3$ of the sentence (in-air; $\mathrm{FI}=0.50$)). With respect to PD duration, the lowest EER was achieved using conventional handwriting features computed for both on-surface as well as in-air movements ($23.64\pm 7.55\%$).

4. Discussion

To the best of our knowledge, except for our pilot work [20], there are no prior studies which integrate FD into a handwriting parameterization for quantitative PD dysgraphia analysis. Therefore, the results published in this paper are exploratory in nature.

In comparison with the conventional kinematic features, FD-based ones correlate more significantly with the clinical characteristics (UPDRS V and PD duration). We observed especially strong correlations for handwriting tasks based on the periodic repetition of specific movements (Archimedean spiral; repetitive letter l, syllable le, or word les). Although the levels of significance based on the conventional handwriting parameters are lower, similar handwriting tasks are involved in the most significant results. We hypothesize that this is due to their ability to highlight or better quantify the cardinal motor symptoms of PD. For example, the most significant relationship between handwriting performance and PD duration was identified in acceleration extracted from the Archimedean spiral. Rigidity combined with tremor and/or bradykinesia makes a PD patient’s handwriting/drawing less fluent (increased changes in velocity and higher acceleration). This is highlighted in a task such as the spiral, where the proper coordination of the fingers, wrist, and arm is required. Generally, the observed problems with coordination are in line with the work of Dounskaia et al. [51] and Teulings et al. [52]. To better illustrate these manifestations, Figure 3 plots the velocity profiles of repetitive loops for a healthy control and a PD patient. As can be seen, the patient introduced more changes in velocity, and their drawing became much more non-fluent. To summarize these findings, FD features in combination with properly selected tasks provide a stronger relationship with the severity and progress of PD.

On the other hand, in terms of binary classification, the conventional parameters provided the best results. The classification performance is remarkable: $\mathrm{ACC}=97.74\%$, $\mathrm{SEN}=95.50\%$, and $\mathrm{SPE}=100\%$. In fact, our results represent the highest classification accuracy that has ever been reported based on the PaHaW database (see Table 1). We hypothesize that the improvement was caused by the inclusion of the state-of-the-art XGBoost algorithm into our machine learning pipelines. As already mentioned, the result is based on one in-air feature: median horizontal velocity of a sentence. In comparison with the HC cohort, the PD patients exhibited much lower values of this measure, i.e., while writing the sentence, the PD patients were not able to perform horizontal transitions (movement between neighboring letters or words) as quickly as the HC could. This finding is in line with the work of Ma et al. [53], who observed that wrist extension stiffness in PD patients makes the handwriting in the horizontal direction more problematic. Therefore, scientists started to use the term horizontal dysgraphia [13]. Generally, vertical or horizontal dysgraphia may be considered a presymptomatic neurobehavioral biomarker of PD with possible significance in early PD diagnosis [13].

In [20], we proved that the FD features improved the accuracy of PD dysgraphia diagnosis in the Archimedean spiral drawing task by 10%. Contrary to our pilot results, in the frame of this work, these features did not lead to any improvements. After a deeper analysis, we found that this was caused by a combined task approach. Performance of the Archimedean spiral is a quasiparticle and continuous task with some repetitive patterns. It looks as though the FD features work especially well in these specific cases. Nevertheless, when combining these tasks with a complex handwriting task (such as a sentence), the measures quantifying in-air movement tend to be more discriminative (in our case, the median in-air horizontal velocity of a sentence). This brings us to the same conclusion that was given during the correlation analysis—the FD features advance the PD dysgraphia diagnosis only in some specific cases.

The best regression model, estimating the UPDRS V score with a 12.51% error, is based only on the standard deviation of on-surface vertical velocity ($\alpha =0.1$) extracted from the sentence. This FD-based parameter was selected from the feature set combining all on-surface measures; therefore, we can confirm the positive influence of FC on the regression analysis performance. In fact, the FD features outperformed the conventional ones in all scenarios. To better understand this result, we plotted vertical velocity patterns of the sentence task for different orders of FD (see Figure 4). We can observe a big difference between $\alpha =0.1$ and the rest of the orders, including the full derivative. This large distance means that we are working with completely new information that is far from that contained in the full derivative. Although it is difficult to clinically interpret this information, it is clear that FC opens new possibilities for monitoring PD severity.

Regarding the PD duration estimation results, the most successful model ($\mathrm{EER}=23.46\%$) consists of 16 conventional on-surface/in-air features (all features’ importance values can be found in Supplementary Table S1). The most frequent feature with the highest feature importance is the jerk extracted from several handwriting tasks. This probably means that as PD progresses, handwriting becomes more jerky and irregular. Vertical velocity is the second most frequent feature involved in the models, which is probably linked with micrographia. Generally, in the case of PD duration estimation, the FD-based features did not yield any improvement.

In conclusion, the FD-based features are better for modeling PD severity (in terms of UPDRS V score estimation), but they do not lead to an improvement in PD duration modeling. The progress of PD is nonlinear and very individual. This means that patients with the same PD duration can be in different stages of the disease. This fact supports our results: the estimation error of PD duration was generally much worse than the estimation error of the UPDRS V score. Since PD duration estimation is a difficult task with poor results, fine improvements based on FD parameters play no role.

5. Conclusions

This study deals with advanced approaches to PD dysgraphia diagnosis and monitoring based on FC integrated with online handwriting/drawing parameterization. To the best of our knowledge, it is the first work that performs a complex investigation into the possibilities for FC in online handwriting processing and proposes new advances in kinematic analyses based on FD. Although the conventional features provided better and very high classification accuracy, which is at the top of the state-of-the-art analyses based on the PaHaW database ($\mathrm{ACC}=97.74\%$, $\mathrm{SEN}=95.50\%$, and $\mathrm{SPE}=100\%$), the newly designed parameters were proven to work better for specific tasks (continuous and/or repetitive, such as the Archimedean spiral) and for specific applications, i.e., PD severity estimation (EER = 12.51%). However, our results need to be confirmed by subsequent scientific research.

This study has several limitations and suggestions for further improvements. Since the dataset is small, to be able to generalize the results, bigger databases should be involved. On the other hand, it is common to have such small numbers of PD patients and HC samples in PD dysgraphia analysis, e.g., see our review in Table 1. Next, we considered only the kinematic measures. To better evaluate the discrimination power of the FD features and better evaluate their ability to estimate PD severity or progress, other feature types, such as temporal, spatial, and dynamic, should be included in future comparisons. Finally, the FD-based parameters could be further explored. For instance, we can consider other approximations (e.g., Caputo) or employ FC for other measures (e.g., entropies).