 Research
 Open Access
 Published:
Effect of poststroke spasticity on voluntary movement of the upper limb
Journal of NeuroEngineering and Rehabilitation volume 18, Article number: 81 (2021)
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 FuglMeyer assessment scale for the upper extremity (FMAUE). Fortytwo firstevent 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 reachtograsp movements towards four differently located targets. LogBKLD 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 FMAUE 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 poststroke 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 longterm 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 velocitydependent increase in muscle resistance stemming from hyperexcitability of the dysregulated musclespindle activity and stretchreflex 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 6point 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 poortogood 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, multidimensional, 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 maximumlikelihood 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
Fortytwo participants with stroke (28 male, age 53.3 [10.5 SD] years, 21 with lefthemiparesis 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 subacute phase (3 weeks to 6 months post stroke onset), while in a stable clinical and metabolic state. In all patients this was a firstever 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 (ChedokeMcMaster 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 antispasticity medications. Upper limb sensorimotor impairment was assessed with the FuglMeyer assessment scale of the upper extremity (FMAUE) [33], a 66point 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.
Procedure
Participants performed reachtograsp motion toward a hollow cone (6cm 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 midsagittal plane (nearcenter (NC); farcenter (FC)), and at full arm length ~ 30 cm to the right and left of the midsagittal 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 midsternum, midpoint of the acromial superiorlateral 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 semirandomized reachtograsp 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.
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 2way (zero lag), lowpass, thirdorder Butterworth filter with a 6Hz 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:
where x is the original joint angle vector, x_{i} is angle value and x_{i,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 Kmeans [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 nearestneighbor methodology [43].
Averages were calculated for each participant per target of the following common upper limb kinematic measures of reachtograsp 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,
where y_{i,j} is the outcome value of measurement i, for participant j. β is the fixed effects vector, X_{i,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 R^{2} (R_{c}^{2}) and marginal R^{2} (R_{m}^{2}) values were evaluated for all models [47]. The R_{c}^{2} represents the variance explained by both fixed and random factors, and thus indicates how the model fits the participants. The marginal R_{m}^{2} 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 [nearcenter (NC), farcenter (FC), contralateral (CL), ipsilateral (IL)], and their interaction as factors. Since BKLD was rightskewed, log BKLD (logBKLD) 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, FMAUE, 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).
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 logBKLD, were larger than control–control group values (HD: χ^{2} = 17.96, p < 0.001, R_{m}^{2} = 0.13, R_{c}^{2} = 0.88; logBKLD: χ^{2} = 24.27, p < 0.001, R_{m}^{2} = 0.16, R_{c}^{2} = 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 logBKLD, 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 logBKLD values for reaches to the nearcenter target than for all other targets (p < 0.001).
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 R_{c}^{2} (0.88–0.96). HD had the highest R_{m}^{2} (0.49). logBKLD had the second highest R_{m}^{2} (0.45). The R_{m}^{2} values of the rest of the measures were much lower (0.22–0.33).
All measures except HD, peak elbow velocity, and movement time were strongly related to target location (logBKLD: χ^{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 logBKLD values (i.e., were more similar to controls) for reaches to the nearcenter 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 nearcenter 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 nearcenter target was higher than to the ipsilateral target (p < 0.001) and marginally higher than that to the farcenter target (p = 0.05). The average elbow velocity toward the ipsilateral target was the lowest (nearcenter: p < 0.001; farcenter: p < 0.01; contralateral: p < 0.001). Average number of elbow velocity peaks towards the nearcenter target was similar to that toward the ipsilateral but lower than that for the farcenter target (p < 0.01) and the contralateral target (p < 0.001).
All the measures were strongly related to the FMAUE (HD: χ^{2} = 41.10, p < 0.001; logBKLD: χ^{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 FMAUE scores had higher HD and logBKLD values (larger divergence from healthy models) (Fig. 5). For all measures, there was no interaction between FMAUE and TSRT. HD, logBKLD, mean elbow velocity, and peak elbow velocity were related to TSRT (HD: χ^{2} = 8.02, p < 0.01, logBKLD: χ^{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 agreedupon division of TSRT values to subgroups reflecting different degrees of severity, and visualizing the relationship of these measures (HD, logBKLD, mean elbow velocity) to TSRT is not trivial (Fig. 6). In order to examine if the relationship between HD, logBKLD, 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 Ftests 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 nonsevere (mild and moderate) impairment are considered (HD, all data: F_{164,1} = 15.7 p < 0.001, nonsevere data: F_{88,1} = 8.5 p < 0.001; logBKLD, all data: F_{164,1} = 6.9 p < 0.01, nonsevere data: F_{88,1} = 13.4 p < 0.01; mean elbow velocity, all data: F_{164,1} = 5.3 p < 0.05, nonsevere data: F_{88,1} = 4.0 p < 0.05). While participants with more severe impairment have higher HD, higher logBKLD, or lower mean elbow velocity, the quadratic regression fit indicates that participants with midrange TSRT values (~ 100–130°) have higher HD, higher logBKLD, 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 logBKLD 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 logBKLD.
Discussion
Spasticity measured by the TSRT was strongly related to both model distances (HD and logBKLD) 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 logBKLD) and the velocity, both TSRT and FMAUE 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., FMAUE). The LMM for HD had the highest R_{m}^{2}, 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 logBKLD, shed important light on the influence of spasticity on voluntary motion. LogBKLD 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 logBKLDs. It can also explain the influence of target location on logBKLD 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 FMAUE) 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 FMAUE) 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 FMAUE 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 FMAUE 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 logBKLD, suggesting a skewed data distribution. The REML convergence criterion used is robust against skewness in terms of bias [52].
Methods for assisting poststroke 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 inpatient 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 nonadaptive 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 impairmentoriented 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., FMAUE). 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.
Abbreviations
 ANOVA:

Analysis of variance
 BKLD:

Bidirectional Kullback–Leibler divergence
 CL:

Contralateral
 CT:

Computerized tomography
 FC:

Farcenter
 FMAUE:

FuglMeyer assessment scale of the upper extremity
 GMM:

Gaussian mixture model
 HD:

Hellinger’s distance
 IL:

Ipsilateral
 IQR:

Interquartile range
 LMM:

Linear Mixed Model
 MAS:

Modified Ashworth Scale
 MRI:

Magnetic resonance imaging
 NC:

Nearcenter
 SD:

Standard deviation
 TSRT:

Tonic Stretch Reflex Threshold
References
 1.
Langhorne P, Coupar F, Pollock A. Motor recovery after stroke: a systematic review. Lancet Neurol. 2009;8(8):741–54.
 2.
Broeks JG, Lankhorst GJ, Rumping K, Prevo AJ. The longterm outcome of arm function after stroke: results of a followup study. Disabil Rehabil. 1999;21:357–64.
 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.
 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.
 5.
Baude M, Nielsen JB, Gracies JM. The neurophysiology of deforming spastic paresis: a revised taxonomy. Ann Phys Rehabil Med. 2019;62:426–30.
 6.
Wissel J, Manack A, Brainin M. Toward an epidemiology of poststroke spasticity. Neurology. 2013;80(3 suppl2):S13–9.
 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.
 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.
 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.
 10.
Bensmail D, Robertson JVG, Fermanian C, RobyBrami A. Botulinum toxin to treat upperlimb spasticity in hemiparetic patients: analysis of function and kinematics of reaching movements. Neurorehabil Neural Repair. 2010;24(3):273–81.
 11.
Mochizuki G, Centen A, Resnick M, Lowrey C, Dukelow SP, Scott SH. Movement kinematics and proprioception in poststroke spasticity: assessment using the Kinarm robotic exoskeleton. J Neuroeng Rehabil. 2019;16(1):146.
 12.
Wood DE, Burridge JH, van Wijck FM, McFadden C, Hitchcock RA, Pandyan AD, Haugh A, SalazarTorres 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.
 13.
Bohannon RW, Smith MB. Interrater reliability of a modified Ashworth scale of muscle spasticity. Phys Ther. 1987;67(2):206–7.
 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.
 15.
Krakauer JW, Carmichael ST. Chapter 2.10 Broken movement. The neurobiology of motor recovery after stroke. Cambridge: MIT Press; 2017.
 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.
 17.
Turpin N, Feldman AG, Levin MF. Stretchreflex threshold modulation during active elbow movements in poststroke survivors with spasticity. Clin Neurophysiol. 2017;128(10):1891–7.
 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.
 19.
Feldman A. Once more on the equilibriumpoint hypothesis (λ model) for motor control. J Mot Behav. 1986;18(1):17–54.
 20.
Feldman AG. Referent control of action and perception. New York: Springer; 2015.
 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.
 22.
Demers M, Levin MF. Do activity level outcome measures commonly used in neurological practice assess upperlimb movement quality? Neurorehabil Neural Repair. 2017;31(7):623–37.
 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: consensusbased core recommendations from the second stroke recovery and rehabilitation roundtable. Int J Stroke. 2019;14(8):783–91.
 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.
 25.
Rasmussen CE. The infinite Gaussian mixture model. In: Solla SA, Leen TK, Muller KR, editors. Advances in neural information processing systems 12. Massachusetts: MIT Press; 2005. p. 554–60.
 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.
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.
Cutler A, CorderoBrana OI. Minimum Hellinger’s distance estimation for finite mixture models. J Am Stat Assoc. 1996;91(436):1716–23.
 29.
Davidowitz I, Parmet Y, FrenkelToledo S, Banina MC, Soroker N, Solomon JM, Liebermann DG, Levin MF, Berman S. Relationship between spasticity and upperlimb movement disorders in individuals with subacute stroke using stochastic spatiotemporal modeling. Neurorehabil Neural Repair. 2019;33(2):141–52.
 30.
Lackritz H, Parmet Y, FrenkelToledo 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.
Lackritz H, Parmet Y, FrenkelToledo 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, BeerSheva, March 24–25, 2019.
 32.
Gowland C, Stratford P, Ward M, Moreland J, Torresin W, Van Hullenaar S, Plews N. Measuring physical impairment and disability with the ChedokeMcMaster Stroke Assessment. Stroke. 1993;24(1):58–63.
 33.
FuglMeyer AR, Jaasko L, Leyman I, Olsson S, Steglind S. The poststroke hemiplegic patient. 1. A method for evaluation of physical performance. Scand J Rehabil Med. 1975;7(1):13–31.
 34.
Levin MF, Banina MC, FrenkelToledo S, Berman S, Soroker N, Solomon JM, Liebermann DG. Personalized upper limb training combined with anodaltDCS for sensorimotor recovery in spastic hemiparesis: study protocol for a randomized controlled trial. Trials. 2018. https://doi.org/10.1186/s1306301723776.
 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.
Prokopenko RA, Frolov AA, Biryukova EV, RobyBrami A. Assessment of the accuracy of a human arm model with seven degrees of freedom. J Biomech. 2001;34(2):177–85.
 37.
Diebel J. Representing attitude: Euler angles, unit quaternions, and rotation vectors. Matrix. 2006;58(15–16):1–35.
 38.
Bishop CM. Pattern recognition and machine learning. New York: Springer; 2006.
 39.
Specht DF. A general regression neural network. IEEE Trans Neural Netw. 1991;2(6):568–76.
 40.
Cohn DA, Ghahramani Z, Michael IJ. Active learning with statistical models. J Artif Intell Res. 1996;4:129–45.
 41.
Hershey JR, Olsen PA, Rennie SJ. Variational KullbackLeibler divergence for hidden Markov models. In: IEEE workshop automatic speech recognition understanding, Kyoto, Japan, December 9–13, 2007. pp. 323–28.
 42.
Tamura RN, Boos DD. Minimum Hellinger distance estimation for multivariate location and covariance. J Am Stat Assoc. 1986;81(393):223–9.
 43.
Altman NS. An introduction to kernel and nearestneighbor nonparametric regression. Am Stat. 1992;46(3):175–85.
 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.
 45.
Satterthwaite FE. An approximate distribution of estimates of variance components. Biom Bull. 1946;2(6):110–4.
 46.
Cheung MWL. Implementing restricted maximum likelihood estimation in structural equation models. Struct Equ Modeling. 2013;20(1):157–67.
 47.
Nakagawa S, Schielzeth H. A general and simple method for obtaining R2 from generalized linear mixedeffects models. Methods Ecol Evol. 2013;4(2):133–42.
 48.
Friendly M, Monette G, Fox J. Elliptical insights: understanding statistical methods through elliptical geometry. Stat Sci. 2013;28(1):1–39.
 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.
 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.
 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.
 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.
 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.
 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.
 55.
Csiszar I. Informationtype measures of difference of probability distributions and indirect. Stud Sci Math Hung. 1967;2:299–318.
 56.
Qiao Y, Minematsu N. A study on invariance of fdivergence and its application to speech recognition. IEEE TRANs Sig Proc. 2010;58(7):3884–990.
 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.
Julier S, Uhlmann J. A general method for approximating nonlinear transformations of probability distributions. Technical Report, Department of Engineering Science, University of Oxford. 1996.
 59.
Kristan M, Leonardis A, Skocaj D. Multivariate online kernel density estimation with Gaussian kernels. Pattern Recogn. 2011;44(10–11):2630–42.
Acknowledgements
The authors acknowledge Rhona Guberek, Maureen McMahon, Franceen Kaizer, MarieTherese 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 CanadaIsrael 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 108186001, ISF Grant number 2392, and by the Helmsley Charitable Trust through the Agricultural, Biological and Cognitive Robotics Center of BenGurion University of the Negev (SB).
Author information
Affiliations
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
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 GingrasLindsay de Montréal (CRIR11121115). The Institutional Ethics Review Board of Loewenstein Rehabilitation Hospital (Ra’anana, Israel) approved the study for the Loewenstein Rehabilitation Hospital (0001115LOE) 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 (IEC32/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:
where w_{i} are convex weights (sum to 1 and are all positive) and g_{i} 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,
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 maximumlikelihood based goodnessoffit 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 fdivergences [54,55,56]. An fdivergence, D_{f}(PQ) quantifies the divergence of the probability distribution Q from the probability distribution P. For continuous random variables,
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 fdivergence between two identical distributions is zero). The fdivergences differ based on the choice of the f function.
Kullback–Leibler divergence
KLD [26, 27] is a commonly used fdivergence for which the f function is,
Therefore, for continuous random variables, KLD is defined by substituting t by p(x)/q(x) and computing the expectation,
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 nonnegative D_{KL}(PQ) ≥ 0, where D_{KL}(PQ) = 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,
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 lowerbound 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 fdivergence. The f function of the squared HD is,
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,
HD is a proper distance metric. It is symmetric and it satisfies the triangle inequality. HD is bounded between 0 and 1, 0 ≤ D_{HD}(PQ) ≤ 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.
About this article
Cite this article
Lackritz, H., Parmet, Y., FrenkelToledo, S. et al. Effect of poststroke spasticity on voluntary movement of the upper limb. J NeuroEngineering Rehabil 18, 81 (2021). https://doi.org/10.1186/s12984021008766
Received:
Accepted:
Published:
Keywords
 Stroke
 Spasticity
 Hemiparesis
 Kinematics
 Stochastic model
 Gaussian mixture model
 Hellinger’s distance
 Kullback–Liebler divergence