Skip to main content

Effect of post-stroke spasticity on voluntary movement of the upper limb

Abstract

Background

Hemiparesis following stroke is often accompanied by spasticity. Spasticity is one factor among the multiple components of the upper motor neuron syndrome that contributes to movement impairment. However, the specific contribution of spasticity is difficult to isolate and quantify. We propose a new method of quantification and evaluation of the impact of spasticity on the quality of movement following stroke.

Methods

Spasticity was assessed using the Tonic Stretch Reflex Threshold (TSRT). TSRT was analyzed in relation to stochastic models of motion to quantify the deviation of the hemiparetic upper limb motion from the normal motion patterns during a reaching task. Specifically, we assessed the impact of spasticity in the elbow flexors on reaching motion patterns using two distinct measures of the ‘distance’ between pathological and normal movement, (a) the bidirectional Kullback–Liebler divergence (BKLD) and (b) Hellinger’s distance (HD). These measures differ in their sensitivity to different confounding variables. Motor impairment was assessed clinically by the Fugl-Meyer assessment scale for the upper extremity (FMA-UE). Forty-two first-event stroke patients in the subacute phase and 13 healthy controls of similar age participated in the study. Elbow motion was analyzed in the context of repeated reach-to-grasp movements towards four differently located targets. Log-BKLD and HD along with movement time, final elbow extension angle, mean elbow velocity, peak elbow velocity, and the number of velocity peaks of the elbow motion were computed.

Results

Upper limb kinematics in patients with lower FMA-UE scores (greater impairment) showed greater deviation from normality when the distance between impaired and normal elbow motion was analyzed either with the BKLD or HD measures. The severity of spasticity, reflected by the TSRT, was related to the distance between impaired and normal elbow motion analyzed with either distance measure. Mean elbow velocity differed between targets, however HD was not sensitive to target location. This may point at effects of spasticity on motion quality that go beyond effects on velocity.

Conclusions

The two methods for analyzing pathological movement post-stroke provide new options for studying the relationship between spasticity and movement quality under different spatiotemporal constraints.

Introduction

Stroke is one of the leading causes of long-term motor disability [1]. Most individuals with stroke present upper limb sensorimotor deficits that persist into the chronic stage (more than 6 months following the onset of stroke) [1, 2]. Spasticity, a sensorimotor disorder characterized by a velocity-dependent increase in muscle resistance stemming from hyperexcitability of the dysregulated muscle-spindle activity and stretch-reflex arc, is a prevalent sensorimotor deficit following stroke [3,4,5]. As many as 20–50% of patients develop spasticity during the first year after the event [6]. Objective and accurate quantification of spasticity and its effects on voluntary motion is important for guiding rehabilitation of the affected limbs.

While there are several clinical measures of spasticity, controversy remains about the most appropriate ones [7, 8]. Moreover, current measures are not sufficient for determining relationships between spasticity, movement deficits, and functional ability [9, 10]. To establish the effects of spasticity on voluntary motion, prior work has attempted to identify the relationship between the amount of hypertonicity measured at rest and movement disruption of voluntarily activated muscle [10,11,12]. One of the commonly used measures of spasticity is the Modified Ashworth Scale (MAS), which grades the resistance felt during passive stretching of muscles on a 6-point ordinal scale [13, 14]. A major drawback of the MAS is that the passive resistance during stretch characterizes only one aspect of the spasticity phenomenon (i.e., amount of hypertonia at rest). In addition, the scale has low resolution and poor-to-good test–retest reliability [13, 14].

Altered muscle resistance at rest may not be the only reason for disruption in voluntary movement [15]. The Tonic Stretch Reflex Threshold (TSRT) identifies where in the biomechanical joint range abnormal muscle resistance begins to contribute to disrupted muscle activation patterns and kinematics [16, 17] providing more specific information about how spasticity and movement deficits are related. TSRT can be determined objectively using the Montreal Spasticity Measure device [18]. The relationship of the TSRT angle with spasticity is based on the threshold control theory of motor control proposed by Feldman [19, 20]. According to the threshold control theory, voluntary movement is generated by regulating the spatial thresholds (of muscle length), at which muscle activation begins. The TSRT, i.e., the spatial threshold at zero velocity, is extrapolated based on the linear regression through measurements representing dynamic spatial thresholds evoked at different stretch velocities.

Upper limb recovery following damage to the brain refers to behavioral restitution (restoring premorbid movement patterns), to which spontaneous restitution is the main contributor. Improvements in upper limb function can also occur through behavioral compensation in which the system accomplishes functional tasks using altered movement patterns [21]. Standard clinical measures of upper limb function do not capture movement quality in a precise manner and therefore are inadequate for differentiating between restitution and compensation [22]. Measurement of movement kinematics and kinetics were suggested as the best way to address this problem. Although some guidelines are available regarding metrics used to characterize motor recovery [23], there is no consensus about how to identify the relationship between spasticity and motor dysfunction. Spasticity is affected by movement velocity and by multiple factors that are difficult to control, such as fatigue, secondary tasks, posture, psychological stress, and time of the day. The variability resulting from these aspects together with the inherent variability of human motion, especially during slow movement, which is typical in people with stroke, complicate the measurement of the effects of spasticity on kinematics. Therefore, a measure that can integrate the spatial and temporal aspects of motion while accounting for motion variability is required.

Stochastic models (based on random variables) offer a comprehensive yet parsimonious representation of motion data. They can capture complex, multi-dimensional, spatiotemporal phenomena while accounting for variability. These features lend stochastic models coupled with a stochastic distance measure (a distance measure between stochastic models) potential advantages over the more commonly used point measures (e.g., mean) for quantifying the effects of spasticity on voluntary motion. However, as the computation of stochastic distance measures is typically more complex, their benefits must be verified. There are several stochastic distance measures that differ according to the characteristics of the differences they capture and the ease of their computation for diverse distributions. Gaussian mixture models (GMMs) are particularly attractive stochastic models for motion representation, since they are easily adapted for spatiotemporal data representation and have standard efficient methods for parameter estimation based on maximum-likelihood estimators via the expectation–maximization algorithm [24] or Bayesian estimation [25]. Two examples of commonly used stochastic measures suitable for measuring distances between GMMs, are the bidirectional Kullback–Leibler divergence (BKLD) and Hellinger’s distance (HD) [26,27,28]. Selecting an appropriate stochastic distance measure is complex since different measures quantify different aspects of dissimilarity between distributions, and thus, are influenced differently by data attributes. HD seems particularly worth exploring in addition to KLD in the context of quantifying the influence of spasticity on voluntary motion. It offers a different perspective regarding the distance between models, and it can rectify some of the shortcomings of the KLD (and BKLD). GMMs, KLD and HD and their shortcomings are described in Appendix 1.

Davidowitz et al. [29] have recently proposed using GMMs and BKLD for quantifying the effects of spasticity (measured by the resistance to passive movement, as reflected in the MAS score) on kinematics of voluntary motion. In a cohort of 16 participants with stroke, spasticity measured by the MAS explained the BKLD of the elbow motion models of reaching movement from nearest neighbor models of healthy individuals. Deviations in individuals with stroke with higher spasticity levels were greater (larger BKLD) than those in individuals with mild spasticity. In the current study, we advanced this effort in two directions. First, we quantified the threshold angle of spasticity using TSRT. In addition, we compared two stochastic distance measures, BKLD and HD, and analyzed the advantages of each for measuring the effects of spasticity on voluntary movement. We hypothesized that both distance measures would be related to TSRT and that differences between the measured distances would highlight different aspects of spasticity, due to the variant sensitivity of BKLD and HD to different confounding variables. Preliminary results have been presented in abstract form [30, 31].

Methods

Participants

Forty-two participants with stroke (28 male, age 53.3 [10.5 SD] years, 21 with left-hemiparesis and 21 with right hemiparesis), and 13 healthy controls of similar age (9 males, mean age 60.5 [8.7 SD] years) participated in the study (Table 1). All patients were tested in the sub-acute phase (3 weeks to 6 months post stroke onset), while in a stable clinical and metabolic state. In all patients this was a first-ever ischemic or hemorrhagic stroke in middle cerebral artery territory (confirmed by MRI or CT scan of the brain) and past medical history was negative for neurological, neuromuscular, or psychiatric problems. All patients had arm paresis (Chedoke-McMaster stroke assessment 2–6/7) [32], were able to perform voluntary elbow extension/flexion movement of at least 30° per direction, had elbow flexor spasticity, and were able to provide informed consent. Individuals were excluded if they had additional neurological, neuromuscular, or orthopedic impairments, pain, difficulty comprehending instructions, or if they were taking anti-spasticity medications. Upper limb sensorimotor impairment was assessed with the Fugl-Meyer assessment scale of the upper extremity (FMA-UE) [33], a 66-point scale. Spasticity was assessed with the TSRT, calculated based on passive stretch. Data were recorded from participants in three countries: Israel, India, and Canada. Participants signed informed consent forms approved by institutional review boards of Loewenstein Rehabilitation Hospital, Raanana, Israel; Kasturba Hospital, Manipal, India and Center for Interdisciplinary Research in Rehabilitation, Montreal, Canada. In all centers the evaluators were experienced physiotherapists who underwent onsite training with the evaluation equipment and protocol.

Table 1 Mean (SD) demographic and clinical data

Procedure

Participants performed reach-to-grasp motion toward a hollow cone (6-cm diameter base) placed on a table at four target locations that required different coordination of upper limb segments (Fig. 1) [29, 34]. Target locations were at 2/3 of arm length and full arm length in the mid-sagittal plane (near-center (NC); far-center (FC)), and at full arm length ~ 30 cm to the right and left of the mid-sagittal plane (contralateral (CL)/ipsilateral (IL), depending on the hemiparetic side) (Fig. 1 bottom). For healthy participants, all targets were reachable without the need for trunk displacement. Arm length was measured with the elbow extended from the medial axillary border to the distal wrist crease. Arm motion was recorded by a wireless electromagnetic tracking system G4 (Polhemus, Colchester, VT) with 5 sensors (120 Hz) each measuring 6 degrees of freedom with respect to a base calibration frame. The static accuracy of the G4 trackers at a 2 m range is 0.75 deg root mean square (RMS) and 0.64 cm RMS.Footnote 1 Sensors were placed on the mid-sternum, midpoint of the acromial superior-lateral border, midpoint of the ventrolateral arm, dorsal forearm (1/3 forearm length proximal to ulnar head), and the index metacarpophalangeal joint. Participants sat on an armless, wooden chair, in front of the table, with feet supported and the trunk unrestricted. The initial arm posture was set with 30° elbow flexion by placement of the third fingertip on an ipsilateral seat height support (Fig. 1 top). Before motion was recorded, participants practiced reaching to each target twice (a total of 8 reaching movements). During the recording, participants performed two sets of 40 semi-randomized reach-to-grasp movements (10 trials towards 4 targets, total 80 trials) balanced in blocks across targets (so that even when not all trials were completed, e.g., due to fatigue, a balanced number of trials per target was recorded). Each trial consisted of a series of movements starting from the initial position: reaching to grasp the cone “as fast and as precisely as possible”; holding the cone for 2 s; lifting the cone toward the chin; returning it to the target position on the table; and finally returning the hand to the initial position. Only the first segment of the sequence, reaching to grasp the cone, was analyzed. The full sequence was done so that the reaching was functional (i.e., natural). Participants rested between trials and blocks as needed. Events during the trials were logged, e.g., when the participant did not succeed in grasping the cone, dropped the cone during the movement, or when the arm collided with the table.

Fig. 1
figure1

Experimental setup: Healthy participants depicted in all images. Left: Overhead schematic representation of target locations. Participants reached to four targets on the table in front of the participant. The far targets (contralateral, far-center, and ipsilateral) were at arm’s length, and the near target (near-center) was at 2/3 arm’s length. The center targets [near-center (NC) and far-center (FC)] were in the midsagittal plane. Contralateral (CL) and ipsilateral (IL) targets were 30 cm along the horizontal axis to the left (CL for right-handed participants) and right (IL for right-handed participants) of the center targets. Middle: participant with the arm in the initial position and the elbow slightly bent. Right: participant grasping a target cone (NC target)

Analysis

Trials were discarded in case of a recording error, e.g., lack of communication between the sensor hub and the computer, or task failure, e.g., motion started prior to the cue. Sensor data were filtered using a standard 2-way (zero lag), low-pass, third-order Butterworth filter with a 6-Hz cutoff. The onset and offset of the first movement were defined as the times at which the forearm tangential velocity exceeded and remained above or decreased and remained below 10% of the peak forearm tangential velocity for at least 0.1 s. Tangential velocity was computed by differentiating position samples. Joint centers were calculated [35, 36] and joint angles were reconstructed using the Cardan angle convention [37].

The spatial and temporal dimensions of the movement trajectories differ in magnitude considerably. Using different value ranges in stochastic modeling would result in a construction of a biased model, where the dimension with larger values would have a larger weight [38]. To avoid this, joint angles were linearly scaled to the range of [− 1, 1], similar to the average task duration:

$${x}_{i,new}=2\frac{{x}_{i}-\mathrm{min}(x)}{\mathrm{max}\left(x\right)-\mathrm{min}(x)}-1$$

where x is the original joint angle vector, xi is angle value and xi,new is the scaled value.

Spatiotemporal GMMs (i.e., GMM parameters) were estimated for each participant per target based on the scaled movements. Model parameters were estimated using the expectation–maximization algorithm initialized using K-means [38]. The expectation–maximization algorithm requires input vectors of equal length. To create equal length vectors a function representing each trajectory was approximated using a general regression neural network [39]. The approximated function was sampled at a constant rate determined for each model based on the average trajectory length (per participant per target). Models with 2 to 25 Gaussian components were tested and Bayesian Information Criteria [40] were used in order to choose the best fit.

BKLD and HD distances between the estimated GMMs were calculated for each target (see Appendix 1 for more detail). BKLD values were computed using variational approximation [41] and HD values were computed using the unscented transform [42]. For the stroke group, distances between GMMs of individuals with stroke and control participants (stroke–control) were calculated. For the control group, distances between the GMMs of control participants and the other control participants were calculated (control–control). The control–control distances serve as a baseline for similar models. This is especially important for BKLD values, which do not have an upper bound. Since movement patterns may differ between healthy individuals, control–control distances indicate the range of acceptable differences between typical motion models. Distances that were considerably larger were considered to represent abnormal movement patterns. In both cases (stroke–control and control–control), the final BKLD or HD scores for each participant were determined as the minimal score based on the nearest-neighbor methodology [43].

Averages were calculated for each participant per target of the following common upper limb kinematic measures of reach-to-grasp tasks [44] which quantify both spatial and temporal motion characteristics: movement time, final elbow extension angle, mean elbow velocity, peak elbow velocity, and the number of elbow velocity peaks (related to lack of smoothness). Average final angle and average mean elbow velocity were calculated using circular mean computation (mean computation suited for a circular quantity). Elbow velocity was computed by differentiating measured joint angles. The peak elbow velocity was computed over the velocity vector. The mean velocity was computed for all values higher than 0.1*maximal velocity to handle cases of segmented motion. Movement time was calculated as the difference between movement offset and onset (of the forearm tangential velocity defined above). The number of elbow velocity peaks was calculated by differentiating the velocity profile and counting the number of times the acceleration profile became and stayed positive for more than 100 ms.

Statistical analysis

Statistical analysis was performed using R Studio IDE for R (version 3.4.2). Analysis was performed using Linear Mixed Models (LMMs) with restricted maximum likelihood (REML) criterion for convergence [45]. LMMs were preferred over the more traditional analysis of variance (ANOVA) since they yield asymptotically efficient estimators, i.e., tend toward being optimal (minimal variance) as the sample size increases for both balanced and unbalanced research designs, while ANOVA produces an optimal estimator only for balanced designs. REML produces variance component estimators closer to ANOVA with a smaller bias than maximum likelihood [46]. All LMMs included participants as random effect,

$${y}_{i,j}=\beta {X}_{i,j}+{\alpha }_{j}+{\epsilon }_{i,j}$$

where yi,j is the outcome value of measurement i, for participant j. β is the fixed effects vector, Xi,j is the observation vector, αj ~ N(0,σα2) is the random effect of participant j, and εi,j ~ N(0,σε2) is the random error.

The Wald test χ2 statistic was calculated for testing the fixed effects of the LMMs. Conditional R2 (Rc2) and marginal R2 (Rm2) values were evaluated for all models [47]. The Rc2 represents the variance explained by both fixed and random factors, and thus indicates how the model fits the participants. The marginal Rm2 represents the variance explained by fixed factors only, and thus indicates how the model fits the general population of people affected by stroke.

Stroke–control and control–control HD and BKLD values were analyzed using LMMs with HD or BKLD as the outcome value, with the group measured divergence (stroke–control, control–control), target [near-center (NC), far-center (FC), contralateral (CL), ipsilateral (IL)], and their interaction as factors. Since BKLD was right-skewed, log BKLD (log-BKLD) was analyzed to adhere to the model assumption of normality. The influence of spasticity on voluntary motion was examined using LMMs with the different measures (stroke–control HD, stroke–control BKLD, movement time, final elbow extension angle, mean elbow velocity, peak elbow velocity, and the number of velocity peaks) as outcome values and TSRT, FMA-UE, target and interactions with the target, were defined as factors.

Results

Healthy controls completed 97% (8% SD) of trials, whereas participants with stroke completed 79% (22% SD) of trials. For controls, 6% (0.7% task failure) and for participants with stroke 22% (16% task failure) of all trials were discarded. The movements made by the control group were faster, smoother, and less variable than the movements made by the stroke group (Fig. 2, Table 2).

Fig. 2
figure2

An example of elbow joint motion traces and respective GMMs. The y-axis shows elbow extension angles scaled to [− 1, 1]. The x-axis shows the time scaled to the average movement time. The traces are overlaid with an ellipsoid-based representation [48] of the respective Gaussian mixture models for a control participant (left) and a participant with stroke (right, FMA = 36, TSRT = 128.3°) for reaches toward each target. Each ellipsoid depicts one distribution component. Ellipsoids are located at the component mean, with axes along the directions of the covariance matrix eigenvectors and show the solutions to \({\left(x-\mu \right)}^{^{\prime}}{\Sigma }^{-1}\left(x-\mu \right)\) of a respective component. NC near-center, FC far-center, CL contralateral, IL ipsilateral

Table 2 Mean (SD) estimates of kinematic measures

The GMMs of the control group had fewer Gaussian components than the models of the stroke group [control 5.77 (0.90), stroke 11.04 (1.31); χ2 = 601.89, p < 0.001] (Fig. 2). The number of components differed between targets (χ2 = 19.62, p < 0.001). In the stroke group, the number of GMM components of models of reaching to the NC target was lower than the amount in movements toward the three far targets, FC, CL and IL targets (p < 0.001), with no difference between the three far targets.

Stroke–control and control–control distances

For both divergences, stroke–control average HD and log-BKLD, were larger than control–control group values (HD: χ2 = 17.96, p < 0.001, Rm2 = 0.13, Rc2 = 0.88; log-BKLD: χ2 = 24.27, p < 0.001, Rm2 = 0.16, Rc2 = 0.89) (Fig. 3). The control–control BKLD values were an order of magnitude smaller. In addition, the stroke–control values had large standard deviations, while the standard deviations of the control–control were low. For HD, there was no effect of target location or interaction between group and target location. For log-BKLD, there was an effect of target location (χ2 = 33.64, p < 0.001) and no interaction between target and group. Participants from both groups had lower log-BKLD values for reaches to the near-center target than for all other targets (p < 0.001).

Fig. 3
figure3

Boxplots of HD (left) and log-BKLD (right) for elbow extension reaching the four targets. Control–control (red) and stroke–control (blue). HD Hellinger’s distance, BKLD bidirectional Kullback–Liebler divergence, NC near-center, FC far-center, CL contralateral, IL ipsilateral. Dots indicate outliers smaller/larger than 1.5 times the inter-quartile range (IQR)

Effects of spasticity on voluntary motion kinematics

Participant’s age, days since stroke onset, country (Israel, India, Canada), gender, and their interactions with target were not significant in linear mixed models (LMMs) for any of the measures. Therefore, these factors were not included in the final LMMs. The LMMs are presented in Table 3. All measures had very high Rc2 (0.88–0.96). HD had the highest Rm2 (0.49). log-BKLD had the second highest Rm2 (0.45). The Rm2 values of the rest of the measures were much lower (0.22–0.33).

Table 3 Wald Chi-square values for the linear mixed models (LMMs)

All measures except HD, peak elbow velocity, and movement time were strongly related to target location (log-BKLD: χ2 = 55.02, p < 0.001; final elbow extension angle: χ2 = 101.72, p < 0.001; mean velocity: χ2 = 27.8, p < 0.001; number of elbow velocity peaks: χ2 = 36.36, p < 0.001) and there was no interaction between target and TSRT (Figs. 3, 4). Individuals with stroke had lower log-BKLD values (i.e., were more similar to controls) for reaches to the near-center target than for all other (far) targets (p < 0.001 for each of the 3 far targets). They had lower final elbow extension angles for the near-center target comparing to all far targets (p < 0.001 for each of the 3 far targets). The average elbow velocity of individuals with stroke for movement to the near-center target was higher than to the ipsilateral target (p < 0.001) and marginally higher than that to the far-center target (p = 0.05). The average elbow velocity toward the ipsilateral target was the lowest (near-center: p < 0.001; far-center: p < 0.01; contralateral: p < 0.001). Average number of elbow velocity peaks towards the near-center target was similar to that toward the ipsilateral but lower than that for the far-center target (p < 0.01) and the contralateral target (p < 0.001).

Fig. 4
figure4

Boxplots by target for the common kinematic measures. Box plots with line over means for a final elbow extension angle; b movement time; c mean elbow velocity; d peak elbow velocity; e number of elbow velocity peaks, per target for healthy control (red) and stroke (blue) groups. NC near-center, FC far-center, CL contralateral, IL ipsilateral. Dots indicate outliers smaller/larger than 1.5 times IQR

All the measures were strongly related to the FMA-UE (HD: χ2 = 41.10, p < 0.001; log-BKLD: χ2 = 32.88, p < 0.001; movement time: χ2 = 14.14, p < 0.001; final elbow angle: χ2 = 10.14, p < 0.01; mean elbow velocity: χ2 = 11.20, p < 0.001; peak elbow velocity: χ2 = 8.03, p < 0.01; number of elbow velocity peaks: χ2 = 10.47, p < 0.01). Participants with lower FMA-UE scores had higher HD and log-BKLD values (larger divergence from healthy models) (Fig. 5). For all measures, there was no interaction between FMA-UE and TSRT. HD, log-BKLD, mean elbow velocity, and peak elbow velocity were related to TSRT (HD: χ2 = 8.02, p < 0.01, log-BKLD: χ2 = 5.11, p < 0.05; mean elbow velocity: χ2 = 11.00, p < 0.001; peak elbow velocity: χ2 = 12.53, p < 0.001). There is no agreed-upon division of TSRT values to sub-groups reflecting different degrees of severity, and visualizing the relationship of these measures (HD, log-BKLD, mean elbow velocity) to TSRT is not trivial (Fig. 6). In order to examine if the relationship between HD, log-BKLD, mean elbow velocity and TSRT is better estimated by a linear or quadratic regression fit, we added a quadratic term of TRST to the regression model and then the new model was compared to the linear model using the F test. The F-tests showed that a quadratic regression fit for TSRT is significantly better than a linear fit for all measures, when all the data points are considered or when data of participants with non-severe (mild and moderate) impairment are considered (HD, all data: F164,1 = 15.7 p < 0.001, non-severe data: F88,1 = 8.5 p < 0.001; log-BKLD, all data: F164,1 = 6.9 p < 0.01, non-severe data: F88,1 = 13.4 p < 0.01; mean elbow velocity, all data: F164,1 = 5.3 p < 0.05, non-severe data: F88,1 = 4.0 p < 0.05). While participants with more severe impairment have higher HD, higher log-BKLD, or lower mean elbow velocity, the quadratic regression fit indicates that participants with mid-range TSRT values (~ 100–130°) have higher HD, higher log-BKLD, or lower mean elbow velocity values than participants with either lower or higher TSRT values. One subject had both severe impairment and very low TSRT but low HD and log-BKLD scores, which was very different from all other participants with severe impairment. When excluding the data of this subject from the data of the severe impairment group, there was no advantage for the quadratic regression fit over the linear fit. Moreover, the linear fit was nearly horizontal for both HD and log-BKLD.

Fig. 5
figure5

Boxplots of Hellinger’s distance (left) and log-BKLD (right) for elbow extension while reaching to each target, in stroke patients with mild, moderate and severe motor impairment, as reflected in FMA-UE scores. FMA-UE scores between 0 and 30, 31 and 50, and 51 and 66, represent severe, moderate, and mild motor impairment, respectively [49, 50]. NC near-center, FC far-center, CL contralateral, IL ipsilateral, FMA-UE Fugl-Meyer assessment scale of the upper extremity, BKLD bidirectional Kullback–Liebler divergence. Dots indicate outliers smaller/larger than 1.5 times IQR

Fig. 6
figure6

HD, log-BKLD, and mean elbow velocity by TSRT. Top: Quadratic regression fit to all the data; Middle: Quadratic regression fit to data from non-severe subjects (FMA-UE scores above 31); Bottom: Linear regression fit to severe data (FMA-UE scores between 0 and 30). Significance of differences between linear and quadratic regression marked below figures. Dots represent individual participant data. Red dots represent severe impairment (FMA-UE scores between 0 and 30), green dots represent moderate impairment (FMA-UE scores between 31 and 50), and blue dots represent mild impairment (FMA-UE scores between 51 and 66) [41, 42]. In the bottom figures, data of a subject with severe impairment and very low TSRT but low HD and log-BKLD scores marked in black. The data of this subject was not used in the regression fit of the severe data. HD Hellinger’s distance, BKLD bidirectional Kullback–Liebler divergence, TSRT tonic stretch reflex threshold, FMA-UE Fugl-Meyer assessment scale of the upper extremity. F-test for quadratic vs linear regression fit presented below each graph; Significance levels: *p < 0.05; **p < 0.01; ***p < 0.001

Discussion

Spasticity measured by the TSRT was strongly related to both model distances (HD and log-BKLD) and to the mean elbow extension velocity. This suggests that spasticity has a major effect on the quality of voluntary movement in the hemiparetic upper limb. For explaining the distance measures (HD and log-BKLD) and the velocity, both TSRT and FMA-UE were required. This suggests that the effects of spasticity (measured by TSRT) go beyond the effects captured by overall clinical assessments of motor impairment (i.e., FMA-UE). The LMM for HD had the highest Rm2, suggesting that the results for HD are more generalizable to the population of people affected by stroke than the results of all other measures.

The different LMMs that characterized the two distances measured, HD and log-BKLD, shed important light on the influence of spasticity on voluntary motion. Log-BKLD along with several other kinematic measures (final elbow extension angle, mean velocity, and the number of velocity peaks) was related to target location, while HD, movement time, and peak elbow velocity were not related to it. Individuals with stroke made much slower upper limb reaching movements (in the current experiment on average 37–38% slower depending on target) than control participants. It is well documented that slow motions are less smooth and more variable than fast short movements [51]. Large motion variability leads to a large dispersion of spatiotemporal data points and consequently to GMMs with more components and larger variances. When such models are approximated by models of condensed data, e.g., a model fit to motion data of healthy participants, much information may be lost (i.e., the distance measured by BKLD). This explains the large differences between control–control and stroke–control log-BKLDs. It can also explain the influence of target location on log-BKLD since elbow motion velocity differed between targets. In contrast, HD quantifies the separability of the models and not the information loss. Therefore, it is less affected by motion variability. Indeed, our results show that HD was related to spasticity and not to target location. This points to an effect of spasticity on motion quality that goes beyond the effects on velocity.

We found that in general, individuals with stroke who had severe impairment (low FMA-UE) exhibited motion that was less similar to the motion of healthy subjects (high HD) regardless of their TSRT. Individuals with stroke who had mild or moderate impairment (high FMA-UE) and either low level of spasticity (very high TSRT) or high level of spasticity (very low TSRT) exhibited motion that was more similar to the motion of healthy subjects (low HD). Recently, Turpin et al. [17] found that TSRT measured during active motion differed from TSRT measured during passive motion, and that the subjects with lower motor impairment severity (higher FMA-UE scores) had greater modulation of TSRT thresholds during passive and active movement. Individuals with severe impairment had smaller changes between active and passive TSRT, i.e., had a lower capability of modulating the stretch reflex threshold. This phenomenon and the differences between the active and passive stretch reflex thresholds may explain the current results. The dissimilarity of the motion of individuals with low FMA-UE regardless of their TSRT could be related to their inability to modulate the stretch reflex threshold.

There were a few outliers (smaller than Q1− 1.5 IQR or larger than Q3+ 1.5 IQR) for all outcome measures in all divisions tested (healthy vs with stroke, by severity, by target). There were no extreme outliers (smaller than Q1− 3IQR or larger than Q3+ 3IQR). The outliers were not removed for the result analysis. The outliers mainly occurred towards the higher range for the final elbow angle, movement time, mean elbow velocity, peak elbow velocity, and the number of elbow velocity peaks, and towards the lower range for both Hellinger’s distance and for log-BKLD, suggesting a skewed data distribution. The REML convergence criterion used is robust against skewness in terms of bias [52].

Methods for assisting post-stroke upper limb recovery have received much attention. However, the identification of effective rehabilitation interventions has been unsatisfactory [1, 53]. The practices currently dominating stroke rehabilitation are focused largely on achieving better performance during activities of daily living through training protocols constrained by the period allocated for in-patient rehabilitation. While this approach is usually successful in regaining the capacity for independent gait, the results of functional rehabilitation of the hemiparetic upper limb are less effective. Thus, widespread efforts have been dedicated to developing novel intervention protocols aimed at restoring sensorimotor function of the hemiparetic upper limb, while minimizing non-adaptive motor compensations. Refined tools to differentiate between true recovery vs. compensation and to quantify and analyze each component of the upper motor neuron syndrome are a critical necessity for impairment-oriented therapies [21, 23]. An advanced understanding regarding the role of spasticity and its effect on voluntary motion, along with the ability to quantify this effect may facilitate improved guidance of recovery efforts.

The main study limitation is that only spasticity in the elbow joint was measured by TSRT, and it is possible that spasticity in muscles spanning adjacent joints may have affected the elbow movement.

Conclusions

The advanced analyses of movement kinematics used in the current study point to a major effect of spasticity (measured by TSRT) on the quality of voluntary movement in the hemiparetic upper limb. This effect goes beyond the effects captured using clinical tools for the assessment of motor impairment (e.g., FMA-UE). Spasticity and impaired movement quality are core components of the upper motor neuron syndrome. Advanced methodologies for their quantitative assessment, enabling better understanding of the relationship between the two, are crucial for segregation of dynamics reflecting true recovery from amelioration of overall arm function based on compensation. Such methodology is also likely to facilitate the development of novel rehabilitation modalities aimed to promote motor recovery.

Availability of data and materials

The datasets generated and/or analyzed during the current study will not be publicly available due to patient confidentiality rules, but anonymized data is available from the corresponding author on reasonable request.

Notes

  1. 1.

    https://polhemus.com/_assets/img/G4_User_Manual_URM10PH238-D.pdf.

Abbreviations

ANOVA:

Analysis of variance

BKLD:

Bidirectional Kullback–Leibler divergence

CL:

Contralateral

CT:

Computerized tomography

FC:

Far-center

FMA-UE:

Fugl-Meyer assessment scale of the upper extremity

GMM:

Gaussian mixture model

HD:

Hellinger’s distance

IL:

Ipsilateral

IQR:

Inter-quartile range

LMM:

Linear Mixed Model

MAS:

Modified Ashworth Scale

MRI:

Magnetic resonance imaging

NC:

Near-center

SD:

Standard deviation

TSRT:

Tonic Stretch Reflex Threshold

References

  1. 1.

    Langhorne P, Coupar F, Pollock A. Motor recovery after stroke: a systematic review. Lancet Neurol. 2009;8(8):741–54.

    PubMed  Article  PubMed Central  Google Scholar 

  2. 2.

    Broeks JG, Lankhorst GJ, Rumping K, Prevo AJ. The long-term outcome of arm function after stroke: results of a follow-up study. Disabil Rehabil. 1999;21:357–64.

    CAS  Article  Google Scholar 

  3. 3.

    Lance JW. Pathophysiology of spasticity and clinical experience with baclofen. In: Feldman RG, Young RR, Koella WP, editors. Spasticity: disordered motor control. Chicago: Year Book Medical Publisher; 1980. p. 185–220.

    Google Scholar 

  4. 4.

    Pandyan AD, Gregoric M, Barnes MP, Wood D, Van Wijck F, Burridge J, Hermens H, Johnson GR. Spasticity: clinical perceptions, neurological realities and meaningful measurement. Disabil Rehabil. 2005;27:2–6.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  5. 5.

    Baude M, Nielsen JB, Gracies J-M. The neurophysiology of deforming spastic paresis: a revised taxonomy. Ann Phys Rehabil Med. 2019;62:426–30.

    PubMed  Article  PubMed Central  Google Scholar 

  6. 6.

    Wissel J, Manack A, Brainin M. Toward an epidemiology of poststroke spasticity. Neurology. 2013;80(3 suppl2):S13–9.

    PubMed  Article  PubMed Central  Google Scholar 

  7. 7.

    Malhotra S, Pandyan AD, Day CR, Jones PW, Hermens H. Spasticity, an impairment that is poorly defined and poorly measured. Clin Rehabil. 2009;23(7):651–8.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  8. 8.

    Burridge JH, Wood DE, Hermens HJ, Voerman GE, Johnson GR, Van Wijck F, Platz T, Gregoric M, Hitchcock R, Pandyan AD. Theoretical and methodological considerations in the measurement of spasticity. Disabil Rehabil. 2005;27(1–2):69–80.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  9. 9.

    Elovic EP, Simone LK, Zafonte R. Outcome assessment for spasticity management in the patient with traumatic brain injury: the state of the art. J Head Trauma Rehabil. 2004;19:155–77.

    PubMed  Article  PubMed Central  Google Scholar 

  10. 10.

    Bensmail D, Robertson JVG, Fermanian C, Roby-Brami A. Botulinum toxin to treat upper-limb spasticity in hemiparetic patients: analysis of function and kinematics of reaching movements. Neurorehabil Neural Repair. 2010;24(3):273–81.

    PubMed  Article  PubMed Central  Google Scholar 

  11. 11.

    Mochizuki G, Centen A, Resnick M, Lowrey C, Dukelow SP, Scott SH. Movement kinematics and proprioception in post-stroke spasticity: assessment using the Kinarm robotic exoskeleton. J Neuroeng Rehabil. 2019;16(1):146.

    PubMed  PubMed Central  Article  Google Scholar 

  12. 12.

    Wood DE, Burridge JH, van Wijck FM, McFadden C, Hitchcock RA, Pandyan AD, Haugh A, Salazar-Torres JJ, Swain ID. Biomechanical approaches applied to the lower and upper limb for the measurement of spasticity: a systematic review of the literature. Disabil Rehabil. 2005;27(1–2):19–32.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  13. 13.

    Bohannon RW, Smith MB. Interrater reliability of a modified Ashworth scale of muscle spasticity. Phys Ther. 1987;67(2):206–7.

    CAS  Article  Google Scholar 

  14. 14.

    Blackburn M, van Vliet P, Mockett SP. Reliability of measurements obtained with the modified Ashworth scale in the lower extremities of people with stroke. Phys Ther. 2002;82:25–34.

    PubMed  Article  PubMed Central  Google Scholar 

  15. 15.

    Krakauer JW, Carmichael ST. Chapter 2.10 Broken movement. The neurobiology of motor recovery after stroke. Cambridge: MIT Press; 2017.

  16. 16.

    Calota A, Feldman AG, Levin MF. Spasticity measurement based on tonic stretch reflex threshold in stroke using a portable device. Clin Neurophysiol. 2008;119(10):2329–37.

    PubMed  Article  Google Scholar 

  17. 17.

    Turpin N, Feldman AG, Levin MF. Stretch-reflex threshold modulation during active elbow movements in post-stroke survivors with spasticity. Clin Neurophysiol. 2017;128(10):1891–7.

    PubMed  Article  Google Scholar 

  18. 18.

    Calota A, Levin MF. Tonic stretch reflex threshold as a measure of spasticity: implications for clinical practice. Top Stroke Rehabil. 2009;16(3):177–88.

    PubMed  Article  PubMed Central  Google Scholar 

  19. 19.

    Feldman A. Once more on the equilibrium-point hypothesis (λ model) for motor control. J Mot Behav. 1986;18(1):17–54.

    CAS  PubMed  Article  Google Scholar 

  20. 20.

    Feldman AG. Referent control of action and perception. New York: Springer; 2015.

    Book  Google Scholar 

  21. 21.

    Levin MF, Kleim JA, Wolf SL. What do motor “recovery” and “compensation” mean in patients following stroke? Neurorehabil Neural Repair. 2009;23(4):313–9.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  22. 22.

    Demers M, Levin MF. Do activity level outcome measures commonly used in neurological practice assess upper-limb movement quality? Neurorehabil Neural Repair. 2017;31(7):623–37.

    PubMed  Article  PubMed Central  Google Scholar 

  23. 23.

    Kwakkel G, Van Wegen EEH, Burridge JH, Winstein CJ, Van Dokkum LEH, Alt Murphy M, Levin MF, Krakauer JW, on behalf of the ADVISORY Group. Standardized measurement of quality of upper limb movement after stroke: consensus-based core recommendations from the second stroke recovery and rehabilitation roundtable. Int J Stroke. 2019;14(8):783–91.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  24. 24

    Dempster AP, Laird NM, Rubin DB. Maximum likelihood from incomplete data via the EM algorithm. J R Stat Soc Ser B (Methodol). 1977;39:1–38.

    Google Scholar 

  25. 25.

    Rasmussen CE. The infinite Gaussian mixture model. In: Solla SA, Leen TK, Muller K-R, editors. Advances in neural information processing systems 12. Massachusetts: MIT Press; 2005. p. 554–60.

    Google Scholar 

  26. 26.

    Jensen JH, Ellis DP, Christensen MG, Jensen SH. Evaluation distance measures between gaussian mixture models of MFCCS. In: Int. Conf. music information retrieval, Vienna, Austria, September 23–27, 2007.

  27. 27.

    Goldberger J, Aronowitz H. A distance measure between GMMs based on the unscented transform and its application to speaker recognition. In: EU. conf. speech comm. tech. Lisbon, Portugal, September 4–8, 2005.

  28. 28.

    Cutler A, Cordero-Brana OI. Minimum Hellinger’s distance estimation for finite mixture models. J Am Stat Assoc. 1996;91(436):1716–23.

    Article  Google Scholar 

  29. 29.

    Davidowitz I, Parmet Y, Frenkel-Toledo S, Banina MC, Soroker N, Solomon JM, Liebermann DG, Levin MF, Berman S. Relationship between spasticity and upper-limb movement disorders in individuals with subacute stroke using stochastic spatiotemporal modeling. Neurorehabil Neural Repair. 2019;33(2):141–52.

    PubMed  Article  PubMed Central  Google Scholar 

  30. 30.

    Lackritz H, Parmet Y, Frenkel-Toledo S, Banina MC, Soroker N, Solomon JM, Liebermann DG, Levin MF, Berman S. Quantifying the effects of spasticity on reaching movement patterns using stochastic spatiotemporal modeling. In: Prog. motor control XII, The Netherlands, Amsterdam, July 7–10, 2019.

  31. 31.

    Lackritz H, Parmet Y, Frenkel-Toledo S, Banina MC, Soroker N, Solomon JM, Liebermann DG, Levin MF, Berman S. Stochastic motion modeling, implemented for measuring the effects of spasticity on kinematics. In: Karniel computational motor control workshop, Beer-Sheva, March 24–25, 2019.

  32. 32.

    Gowland C, Stratford P, Ward M, Moreland J, Torresin W, Van Hullenaar S, Plews N. Measuring physical impairment and disability with the Chedoke-McMaster Stroke Assessment. Stroke. 1993;24(1):58–63.

    CAS  Article  Google Scholar 

  33. 33.

    Fugl-Meyer AR, Jaasko L, Leyman I, Olsson S, Steglind S. The post-stroke hemiplegic patient. 1. A method for evaluation of physical performance. Scand J Rehabil Med. 1975;7(1):13–31.

    CAS  PubMed  PubMed Central  Google Scholar 

  34. 34.

    Levin MF, Banina MC, Frenkel-Toledo S, Berman S, Soroker N, Solomon JM, Liebermann DG. Personalized upper limb training combined with anodal-tDCS for sensorimotor recovery in spastic hemiparesis: study protocol for a randomized controlled trial. Trials. 2018. https://0-doi-org.brum.beds.ac.uk/10.1186/s13063-017-2377-6.

    Article  PubMed  PubMed Central  Google Scholar 

  35. 35.

    O'Brien JF, Bodenheimer Jr RE, Brostow GJ, Hodgins JK. Automatic joint parameter estimation from magnetic motion capture data. In: Graphics Interface Conf. Montreal, Quebec, Canada, May 15–17, 2000. pp. 53–60.

  36. 36.

    Prokopenko RA, Frolov AA, Biryukova EV, Roby-Brami A. Assessment of the accuracy of a human arm model with seven degrees of freedom. J Biomech. 2001;34(2):177–85.

    CAS  PubMed  Article  Google Scholar 

  37. 37.

    Diebel J. Representing attitude: Euler angles, unit quaternions, and rotation vectors. Matrix. 2006;58(15–16):1–35.

    Google Scholar 

  38. 38.

    Bishop CM. Pattern recognition and machine learning. New York: Springer; 2006.

    Google Scholar 

  39. 39.

    Specht DF. A general regression neural network. IEEE Trans Neural Netw. 1991;2(6):568–76.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  40. 40.

    Cohn DA, Ghahramani Z, Michael IJ. Active learning with statistical models. J Artif Intell Res. 1996;4:129–45.

    Article  Google Scholar 

  41. 41.

    Hershey JR, Olsen PA, Rennie SJ. Variational Kullback-Leibler divergence for hidden Markov models. In: IEEE workshop automatic speech recognition understanding, Kyoto, Japan, December 9–13, 2007. pp. 323–28.

  42. 42.

    Tamura RN, Boos DD. Minimum Hellinger distance estimation for multivariate location and covariance. J Am Stat Assoc. 1986;81(393):223–9.

    Article  Google Scholar 

  43. 43.

    Altman NS. An introduction to kernel and nearest-neighbor nonparametric regression. Am Stat. 1992;46(3):175–85.

    Google Scholar 

  44. 44.

    Schwarz A, Kanzler CM, Lambercy O, Luft AR, Veerbeek JM. Systematic review on kinematic assessments of upper limb movements after stroke. Stroke. 2019;50(3):718–27.

    Article  Google Scholar 

  45. 45.

    Satterthwaite FE. An approximate distribution of estimates of variance components. Biom Bull. 1946;2(6):110–4.

    CAS  Article  Google Scholar 

  46. 46.

    Cheung MWL. Implementing restricted maximum likelihood estimation in structural equation models. Struct Equ Modeling. 2013;20(1):157–67.

    Article  Google Scholar 

  47. 47.

    Nakagawa S, Schielzeth H. A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods Ecol Evol. 2013;4(2):133–42.

    Article  Google Scholar 

  48. 48.

    Friendly M, Monette G, Fox J. Elliptical insights: understanding statistical methods through elliptical geometry. Stat Sci. 2013;28(1):1–39.

    Article  Google Scholar 

  49. 49.

    Duncan PW, Goldstein LB, Horner RD, Landsman PB, Samsa GP, Matchar DB. Similar motor recovery of upper and lower extremities after stroke. Stroke. 1994;25(6):1181–8.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  50. 50.

    Subramanian SK, Yamanaka J, Chilingaryan G, Levin MF. Validity of movement pattern kinematics as measures of arm motor impairment poststroke. Stroke. 2010;41:2303–8.

    PubMed  Article  PubMed Central  Google Scholar 

  51. 51.

    Park SW, Marino H, Charles SK, Sternad D, Hogan N. Moving slowly is hard for humans: limitations of dynamic primitives. J Neurophysiol. 2017;118(1):69–83.

    PubMed  PubMed Central  Article  Google Scholar 

  52. 52.

    Banks BD, Mao IL, Walter JP. Robustness of the restricted maximum likelihood estimator derived under normality as applied to data with skewed distributions. J Dairy Sci. 1985;68(7):1785–92.

    Article  Google Scholar 

  53. 53.

    Bernhardt J, Hayward KS, Kwakkel G, Ward NS, Wolf SL, Borschmann K, Krakauer JW, Boyd LA, Carmichael ST, Corbett D, Cramer SC. Agreed definitions and a shared vision for new standards in stroke recovery research: the stroke recovery and rehabilitation roundtable taskforce. Neurorehabil Neural Rep. 2017;31(9):793–9.

    Article  Google Scholar 

  54. 54.

    Ali SM, Silvey SD. A general class of coefficients of divergence of one distribution from another. J R Stat Soc Ser B. 1966;28(1):131–42.

    Google Scholar 

  55. 55.

    Csiszar I. Information-type measures of difference of probability distributions and indirect. Stud Sci Math Hung. 1967;2:299–318.

    Google Scholar 

  56. 56.

    Qiao Y, Minematsu N. A study on invariance of f-divergence and its application to speech recognition. IEEE TRANs Sig Proc. 2010;58(7):3884–990.

    Article  Google Scholar 

  57. 57.

    Hershey JR, Olsen PA. Approximating the Kullback Leibler divergence between gaussian mixture models. In: IEEE Int. conf. acoustics, speech and signal proc. Honolulu, Hawaii, USA, April 15–20, 2007.

  58. 58.

    Julier S, Uhlmann J. A general method for approximating non-linear transformations of probability distributions. Technical Report, Department of Engineering Science, University of Oxford. 1996.

  59. 59.

    Kristan M, Leonardis A, Skocaj D. Multivariate online kernel density estimation with Gaussian kernels. Pattern Recogn. 2011;44(10–11):2630–42.

    Article  Google Scholar 

Download references

Acknowledgements

The authors acknowledge Rhona Guberek, Maureen McMahon, Franceen Kaizer, Marie-Therese Laramée, Arel Shasha, Tal Galinka, Akash Shah, Réjean Prévost, and Subramanian Durairaj for their invaluable contributions to the success of this project.

Funding

This project is supported by the Canada-Israel Health Research Program (MFL and DGL), a program that is jointly funded by: The Canadian Institutes of Health Research, Azrieli Foundation, International Development Research Center and Israel Science Foundation; IDRC Grant number 108186-001, ISF Grant number 2392, and by the Helmsley Charitable Trust through the Agricultural, Biological and Cognitive Robotics Center of Ben-Gurion University of the Negev (SB).

Author information

Affiliations

Authors

Contributions

SB, HL, YP, MFL, DGL, and NS conceptualized the project and the methodology. MFL, DG, and SB acquired the funding and secured the resources. SFT, MCB, and JMS coordinated the data acquisition. SB and HL performed the data analysis including the writing of data analysis programs. YP provided consultation about statistical analysis. SB and HL wrote the initial draft of the manuscript. All authors contributed to the review and editing of the final draft of the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Sigal Berman.

Ethics declarations

Ethics approval and consent to participate

Ethical approval was obtained from the appropriate ethics committees of each study site prior to the beginning of the study. The Centre de Recherche Interdisciplinaire en Réadaptation du Montréal Métropolitain (Montréal, Canada) approved the study for the Jewish Rehabilitation Hospital and the Institut de Réadaptation Gingras-Lindsay de Montréal (CRIR1112-1115). The Institutional Ethics Review Board of Loewenstein Rehabilitation Hospital (Ra’anana, Israel) approved the study for the Loewenstein Rehabilitation Hospital (000-11-15-LOE) and the Institutional Ethics Committee of the Tel Aviv University (Tel Aviv, Israel). The Institutional Ethics Committee of Kasturba Hospital (Manipal, India) approved the study for the Kasturba Medical Hospital (IEC-32/2016). Written informed consent was obtained from each participant.

Consent for publication

Not applicable.

Competing interests

MFL held a US patent for the Montreal Stretch Reflex Threshold analysis. The remaining authors declare that they have no competing interests or any conflicts of interest in the authorship or publication of this study.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Appendix 1: Distance measures between GMMs

Appendix 1: Distance measures between GMMs

Gaussian mixture model

A GMM is defined by a convex sum of K Gaussian probability densities:

$$p\left(x\right)=\sum_{i=1}^{K}{w}_{i}{g}_{i}\left(x|{\mu }_{i},{\Sigma }_{i}\right)$$

where wi are convex weights (sum to 1 and are all positive) and gi are Gaussian density functions with expectation μi and variance Σi.

The trajectory of a single joint angle has two dimensions (time and space). Therefore, it can be represented by a mixture of spatiotemporal Gaussian density functions,

$${g}_{i}\left(x|{\mu }_{i},{\Sigma }_{i}\right)=\frac{1}{{\left(2\pi \right)}^\frac{2}{2}\sqrt{\left|{\Sigma }_{i}\right|}}exp\left\{-0.5{\left(x-{\mu }_{i}\right)}^{T}{\left({\Sigma }_{i}\right)}^{-1}\left(x-{\mu }_{i}\right)\right\}$$
$${\mu }_{i}={\left({\mu }_{t,i}^{T}, {\mu }_{s,i}^{T}\right)}^{T}, {\Sigma }_{i}=\left(\begin{array}{cc}{\Sigma }_{tt,i}& {\Sigma }_{ts,i}\\ {\Sigma }_{ts,i}& {\Sigma }_{ss,i}\end{array}\right)$$

where μs,i is the spatial expectation vector, μt,i is the temporal expectation vector of the ith component, Σtt,i and Σss,i are the variance matrices of temporal and spatial elements, and Σts,i is the covariance matrix between them for the ith component. T is the transpose operation.

A GMM can be learned based on data from several repetitions of a motion trajectory (using parameter estimation methods such as expectation–maximization or Bayesian estimation [19, 20]). Motion quality can be assessed by comparing a motion model performed by an individual with stroke to models of motion performed by healthy individuals of a similar age using a stochastic distance measure. Due to the multivariate nature of the GMM, the classical maximum-likelihood based goodness-of-fit measures (e.g., the χ2 test) cannot be used. An alternative method for measuring the dissimilarity (divergence) between two probability distributions is based on using a measure from the family of invariant measures called f-divergences [54,55,56]. An f-divergence, Df(P||Q) quantifies the divergence of the probability distribution Q from the probability distribution P. For continuous random variables,

$${D}_{f}\left(P||Q\right)={\int }_{\Omega }f\left(\frac{p\left(x\right)}{q\left(x\right)}\right)q\left(x\right)dx$$

where p(x) and q(x) are the probability density functions of P and Q (respectively) on the measurable space Ω. f:(0…∞) → R is a real and convex function and f(1) = 0 (so the f-divergence between two identical distributions is zero). The f-divergences differ based on the choice of the f function.

Kullback–Leibler divergence

KLD [26, 27] is a commonly used f-divergence for which the f function is,

$${f}_{KL}=tln\left(t\right).$$

Therefore, for continuous random variables, KLD is defined by substituting t by p(x)/q(x) and computing the expectation,

$${D}_{KL}\left(P||Q\right)={\int }_{-\infty }^{\infty }p\left(x\right)ln\left(\frac{p\left(x\right)}{q\left(x\right)}\right)dx.$$

KLD quantifies the divergence of Q from P, by taking the expectation with probability P of the log likelihood ratio between probabilities P and Q. The KLD is high when the difference between P and Q is large in regions for which P is high. KLD measures how much information is lost when P is approximated by Q and is thus appropriately called relative entropy. KLD is non-negative DKL(P||Q) ≥ 0, where DKL(P||Q) = 0 indicates that no information is lost when P is approximated by Q. Larger KLD values indicate more information loss, i.e., larger divergence. However, the interpretation of the KLD values is not straightforward as KLD does not have an upper bound.

KLD is not a distance measure as it is not symmetric, and it does not satisfy the triangle inequality. A symmetric variant of the KLD, the BKLD can be computed and used as a distance measure,

$$BKLD=\frac{{D}_{KL}\left(P||Q\right)+{D}_{KL}\left(Q||P\right)}{2}.$$

However, although it is symmetric, the BKLD is not a proper distance measure since it does not satisfy the triangle inequality. This may lead to an increased computation load in cases of multiple comparisons. KLD between GMMs cannot be analytically computed. However, there are several suitable approximation methods, including variational approximation, Monte Carlo simulation, Gaussian approximation, and lower-bound approximation [57]. Among these, variational approximation provides a good balance between accuracy and computation time [41].

Hellinger’s distance

HD is also a commonly used f-divergence. The f function of the squared HD is,

$${f}_{HD}={\left(\sqrt{t}-1\right)}^{2}.$$

Therefore, for continuous random variables, the squared HD is defined by substituting t by p(x)/q(x) and taking the expectation, divided by 2 for simplifying equation,

$${D}_{HD}^{2}\left(P||Q\right)=\frac{1}{2}{\int }_{-\infty }^{\infty }{\left(\sqrt{p(x)}-\sqrt{q(x)}\right)}^{2}dx=1-{\int }_{-\infty }^{\infty }\sqrt{p\left(x\right)q\left(x\right)}dx.$$

HD is a proper distance metric. It is symmetric and it satisfies the triangle inequality. HD is bounded between 0 and 1, 0 ≤ DHD(P||Q) ≤ 1. HD is 0 if (and only if) P and Q are identical and HD is 1 if (and only if) P and Q are mutually exclusive. HD is a measure of the separability between samples from the two distributions. HD is larger when P and Q differ (and thus the integral of \(\sqrt{p\left(x\right)q(x)}\) is smaller). HD has been shown to be invariant, consistent, and robust [50]. HD between GMMs cannot be computed analytically, however, the unscented HD [42, 58] provides a highly accurate estimate [59].

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Lackritz, H., Parmet, Y., Frenkel-Toledo, S. et al. Effect of post-stroke spasticity on voluntary movement of the upper limb. J NeuroEngineering Rehabil 18, 81 (2021). https://0-doi-org.brum.beds.ac.uk/10.1186/s12984-021-00876-6

Download citation

Keywords

  • Stroke
  • Spasticity
  • Hemiparesis
  • Kinematics
  • Stochastic model
  • Gaussian mixture model
  • Hellinger’s distance
  • Kullback–Liebler divergence