IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 61, NO. 8, AUGUST 2014

2235

Position-Dependent Characterization of Passive Wrist Stiffness Autumn L. Pando, Hyunglae Lee, Will B. Drake, Neville Hogan, and Steven K. Charles∗

Abstract—Because the dynamics of wrist rotations are dominated by stiffness, understanding wrist rotations requires a thorough characterization of wrist stiffness in multiple degrees of freedom. The only prior measurement of multivariable wrist stiffness was confined to approximately one-seventh of the wrist range of motion (ROM). Here, we present a precise nonlinear characterization of passive wrist joint stiffness over a range three times greater, which covers approximately 70% of the functional ROM of the wrist. We measured the torque–displacement vector field in 24 directions and fit the data using thin-plate spline smoothing optimized with generalized cross validation. To assess anisotropy and nonlinearity, we subsequently derived several different approximations of the stiffness due to this multivariable vector field. The directional variation of stiffness was more pronounced than reported previously. A linear approximation (obtained by multiple linear regression over the entire field) was significantly more anisotropic (eigenvalue ratio of 2.69 ± 0.52 versus 1.58 ± 0.39; p < 0.001) though less misaligned with the anatomical wrist axes (12.1 ± 4.6◦ versus 21.2 ± 9.2◦ ; p < 0.001). We also found that stiffness over this range exhibited considerable nonlinearity—the error associated with a linear approximation was 20–30%. The nonlinear characterization over this greater range confirmed significantly greater stiffness in radial deviation compared to ulnar deviation. This study provides a characterization of passive wrist stiffness better suited to investigations of natural wrist rotations, which cover much of the wrist’s ROM. It also provides a baseline for the study of neurological and/or orthopedic disorders that result in abnormal wrist stiffness. Index Terms—Impedance, nonlinear, passive, resistance, stiffness, wrist.

Manuscript received November 21, 2013; revised January 22, 2014 and March 18, 2014; accepted March 20, 2014. Date of publication March 25, 2014; date of current version July 15, 2014. The work of H. Lee was supported in part by a Samsung Fellowship. The work of N. Hogan was supported in part by the Eric P. and Evelyn E. Newman Fund. A. L. Pando and H. Lee contributed equally to this work. Asterisk indicates corresponding author. A. L. Pando was with the Department of Mechanical Engineering, Brigham Young University, Provo, UT 84602 USA. She is now with Quicken Loans Inc. Scottsdale, AZ 85260 USA (e-mail: [email protected]). H. Lee was with the Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139 USA. He is now with the Sensory Motor Performance Program, Rehabilitation Institute of Chicago, Chicago, IL 60611 USA (e-mail: [email protected]). W. B. Drake was with the Department of Mechanical Engineering, Brigham Young University, Provo, UT 84602 USA. He is now with B-K Manufacturing Arab, AL 35016 USA (e-mail: [email protected]). N. Hogan is with the Department of Mechanical Engineering and the Department of Brain and Cognitive Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139 USA (e-mail: [email protected]). ∗ S. K. Charles is with the Department of Mechanical Engineering and the Neuroscience Center, Brigham Young University, Provo, UT 84602 USA (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TBME.2014.2313532

I. INTRODUCTION ECADES of upper limb motor control studies have investigated how the neuromuscular system controls the inertial dynamics of shoulder and elbow movements (reaching), recent studies have shown that the dynamics of wrist movements are dominated by stiffness, not inertia [1]. Passive wrist stiffness is greater in radial–ulnar deviation (RUD) than in flexion– extension (FE) and does not align with the anatomical axes [2], creating nontrivial dynamics for which the neuromuscular system must compensate in order to make coordinated wrist movements. The effects of stiffness are observable in wrist movement behavior: incomplete compensation for the passive stiffness of the wrist creates the pattern of path curvature observed in wrist movements [3], [4]. Despite its dominant role in wrist movements, the passive stiffness of the wrist has only been measured in a small portion of the range of motion (ROM) of the wrist joint. Although the ROM of the wrist joint ranges from approximately 77◦ of flexion (FLX) to 67◦ of extension (EXT) and from 21◦ of radial deviation (RD) to 36◦ of ulnar deviation (UD) [5], the only previous in vivo measurement of passive wrist stiffness in both degrees of freedom [2] was limited to 17◦ movements in FE and RUD, which is approximately one-seventh of the wrist ROM and does not permit reliable assessment of nonlinearity [see Fig. 1(a)]. In contrast, natural wrist movements generally involve large portions of the wrist ROM [6], [7]. The purpose of this study was to characterize the passive stiffness of the wrist joint over a much greater range than previously measured to facilitate investigations of natural wrist behavior. We measured and analyzed passive wrist stiffness1 in 15 subjects from 37◦ FLX to 36◦ EXT and from 16◦ RD to 28◦ UD [see Fig. 1(a)], which covers approximately 42% of the ROM of the wrist and 70% of wrist motion during activities of daily living [6]. We found the stiffness over this range to exhibit considerable nonlinearity, with significantly greater stiffness in RD than UD. The directional variation of stiffness differed significantly from that measured over the smaller range [2], being both more anisotropic and less misaligned with respect to the anatomical axes (p < 0.001 in both cases). These data—the nonlinear vector field and its linear approximations—provide a more accurate and detailed characterization of the wrist to

D

1 Because stiffness is a locally linear approximation of the torquedisplacement relationship, a nonlinear torque-displacement relationship is more correctly termed “nonlinear static mechanical impedance” than “nonlinear stiffness” (static mechanical impedance refers to the static component of the relationship between displacement and the torque it evokes). However, because “stiffness” is a simpler and more familiar term, we use it throughout this paper to describe both linear and nonlinear torque–displacement relationships.

0018-9294 © 2014 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications standards/publications/rights/index.html for more information.

2236

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 61, NO. 8, AUGUST 2014

USA). This robot is capable of moving the wrist independently in FE, RUD, and forearm pronation–supination (PS) [8]. The robot sensed displacement in each degree of freedom (DOF) at 200 Hz with a resolution of 0.0006◦ and was able to apply 1.95 Nm of torque in FE and RUD. Each subject was seated with their right arm attached to the robot in a parasagittal plane. The shoulder was in approximately 0◦ of abduction, 25◦ of FLX, and 0◦ of internal rotation, with the elbow flexed at approximately 65◦ to place the forearm in the horizontal plane. The subject’s distal forearm was fastened with a strap to a custommade attachment plate [indicated by “1” in Fig. 1(b)] fixed to the PS stage of the robot (which remained stationary during the experiment). The strap was attached to the PS stage [at “4” in Fig. 1(c)] and this plate [“1” in Fig. 1(b)] and pressed the dorsal aspect of the distal forearm (more specifically the dorsal tubercle of the radius and the dorsal-most protuberance of the ulnar head) against this plate. We found previously [3] that this arrangement places the forearm in a neutral position compatible with the ISB standard [9] and minimizes movement in PS. The hand was attached at the metacarpals to the robot end-effector via a custom-made interface [see “2” in Fig. 1(b)]. This interface attached to the linear bearing of the robot [see “3” in Fig. 1(b)] and was machined to be concave on the side on which it interfaced with the hand to wrap around the back of the hand. Importantly, this interface allowed the fingers to be free and the wrist and finger muscles to relax as much as possible. To record the amount of muscle activity during the stiffness measurement, surface EMG sensors (Trigno by Delsys, Boston, MA, USA) were placed over the flexor carpi radialis, flexor carpi ulnaris, extensor carpi radialis, and extensor carpi ulnaris muscles and set to sample 200 times per second. Fig. 1. (a) Range of motion over which stiffness was measured in this study (“Current range,” filled black circles), compared to the only previous 2-DOF measurement of wrist stiffness [2] (“Prior range,” dashed circle). Also shown are the approximate full range of motion of the wrist (“Wrist ROM,” solid black line, adapted from [5]) and the robot (“Robot ROM,” dashed rectangle). The actual displacement (filled black circles) differed from the programmed targets (“Target,” open circles) because of the control system and torque limit of the robot motors (see Methods). Negative extension and radial deviation represent flexion and ulnar deviation, respectively. (b–c) Photo of setup. Each subject’s distal forearm was strapped to the pronation-supination (PS) stage of the robot (which remained stationary during the experiment) via a strap attached to the PS stage (4) and a custom-made plate (1). The hand was attached at the metacarpals (2) to the end-effector of the wrist robot (3).

support studies of wrist motor control, and also a baseline for evaluating wrist pathologies. II. METHODS A. Subjects Fifteen right-handed individuals (seven male, eight female; age range 20–27 years; BMI range 16.8–24.7) with no selfreported neuromuscular or biomechanical disorders participated in this study. Following procedures approved by Brigham Young University’s Institutional Review Board, informed consent was obtained from all subjects. B. Experimental Setup Stiffness was measured using an InMotion Wrist Rehabilitation Robot (Interactive Motion Technologies, Watertown, MA,

C. Protocol The robot rotated each subject’s wrist between a center target at neutral wrist position (defined below) and 24 peripheral targets requiring FE, RUD, or combinations [see Fig. 1(a)]. The movement cycle began in pure EXT and proceeded counterclockwise, with five repetitions per target before moving counterclockwise to the next target. Targets were placed on the periphery of a rectangle close to the ROM limits of the wrist robot. The robot was commanded to reach each target via proportional-derivative control (with proportional and derivative gains of 10 Nm/rad and 0.1 Nms/rad, respectively). The robot proceeded toward each target until it either reached the target (within the error allowed by the controller gains) or the torque limit of the robot motors (1.95 Nm in FE and RUD) [see Fig. 1(a)]. The torques applied by the robot motors were estimated from the voltages sent to the servo-amplifiers (essentially voltage-controlled current sources for the motors). From the estimated motor torques, we calculated the torques at the robot end-effector. The relationships between the voltage sent to the servo-amplifiers, the current sent to the motors, and the resulting torques were previously characterized in detail for the similar prototype robot [10]. FE and RUD were defined to be in neutral position when the long axis of the third metacarpal aligned with the long axis of the

PANDO et al.: POSITION-DEPENDENT CHARACTERIZATION OF PASSIVE WRIST STIFFNESS

forearm [9]. More specifically, the wrist was in neutral FE and RUD when the elbow joint center (EJC), the wrist joint center (WJC), and the center of the head of the third metacarpal were aligned. The EJC was defined midway between the medial and lateral epicondyles. The WJC was defined proximodistally midway between the distal end of the radius and the proximal end of the third metacarpal, and mediolaterally and ventrodorsally midway between the medial and lateral aspects and midway between the ventral and dorsal aspects of the distal forearm, respectively. When moving a joint from rest, the initial stiffness is often significantly larger than the stiffness of the remaining portion of the movement [11], [12]. To enable the later removal of this short-range stiffness effect, we commanded the robot to start each outbound movement 2◦ “behind” the center target (in the direction opposite of the peripheral target) and make one continuous movement from there through the center target to the peripheral target. Inbound movements were commanded to start from the peripheral target and travel through the center target to 2◦ on the other side of the center target, from where the next repetition immediately followed. Subjects were asked to relax their arm throughout the experiment. To avoid evoking reflex action, the robot was programmed to rotate the wrist smoothly (following a minimum jerk profile) and at a low average velocity of 5.2◦ /s. The actual average velocity of the robot (different from the programmed velocity because of the controller dynamics) was 4.3◦ /s. The activity of the four prime wrist muscles was monitored to ensure that the measured stiffness was indeed passive. The entire stiffness measurement lasted 35 min. After the robot measured stiffness, we measured subjects’ maximum voluntary contraction (MVC) three times in each of the four instrumented wrist muscles by asking subjects to flex or extend their wrist as hard as they could (close to neutral FE) while we resisted by pressing against their hand and supporting their distal forearm. D. Data Processing To remove the effects of short-range stiffness from the data (see above), the displacement and torque data of the first 2◦ were removed from each movement. The recorded torque included not only the torque required to overcome the stiffness of the wrist joint, but also the torque necessary to overcome the dynamics of the robot and the gravitational effects of the hand. Because these nonstiffness-related torques are position dependent, they would affect the estimated stiffness unless removed. These torques are subject specific, so we measured them separately for each subject by rerunning the entire protocol with a weight (in the place of the subject’s hand) under similar conditions: the weight was matched to the subject’s hand mass and placed at the center of mass of the subject’s hand (estimated using regression equations from [13]) relative to the robot, and the robot FE and RUD offsets were the same as in the subject’s original measurement. The nonstiffness-related torques were then subtracted from the total recorded torque data. We also observed effects due to hysteresis, which required greater torques for outbound movements (loading responses) than inbound move-

2237

ments (unloading responses). For inbound movements, hysteresis created a discontinuity around the neutral position, which we removed by shifting the inbound torque by the mean distance between outbound and inbound torque (this was done separately for each movement direction). E. Analysis 1) EMG: After subtracting the mean EMG from the raw EMG (separately for each movement), we rectified and low-pass filtered the remaining EMG signal using a second-order Butterworth filter with cut-off frequency at 3 Hz. We then averaged the EMG over the duration of each movement and normalized by the MVC for each muscle. We then compared this average to the baseline level, estimated as the average EMG over the 1-s period at the beginning of each movement before the wrist had displaced more than 0.5% of its total displacement (we excluded movements whose displacement after 1 s exceeded 0.5% of the total displacement). 2) Stiffness Estimation: The multivariable torque–angle relation of the passive wrist was represented as a vector field V : (τFE , τRUD ) = V (θFE , θRUD )

(1)

where θFE and θRUD are angular displacements in the FE and RUD directions, respectively, and τFE and τRUD are the corresponding applied torques. The components of the vector field were estimated as scalar functions φ1 and φ2 based on thin-plate spline (TPS) smoothing [14] with generalized cross validation (GCV) [15]: τFE = φ1 (θFE , θRUD )

(2)

τRUD = φ2 (θFE , θRUD ) .

(3)

This method provides an optimal compromise between fidelity to the data and roughness of the solution, and the method is sufficiently robust to eliminate the effect of measurement noise, even when the amount of noise is unknown [15]. A more detailed explanation of the vector field decomposition is provided in Appendix A of [16]. Each of the five measurement repetitions was approximated separately by TPS smoothing with GCV (described above) and averaged into a single vector field to estimate a continuous field (separately for outbound and inbound data). This continuous vector field approximation (CVFA) is a complete characterization of wrist stiffness. To simplify further analysis, we used it to evaluate local stiffness at any point of interest as a linear approximation to the torque–displacement vector field. In particular, we computed a “directional stiffness,” Knon−linear for each of the 24 directions, calculated as the slope of the least squares linear fit between displacement in that direction and the component of the restoring torque in that direction.2 Included in the fit were all points (from the CVFA) in a given direction between neutral position and the edge of the CVFA, which was bounded by a box with sides at 14◦ in RD, 29◦ in UD, and 37◦ in FLX and EXT (the values for the bounding box were chosen based 2 Note that this directional stiffness ignores torque components orthogonal to the displacement and hence cannot be used to assess nonconservative behavior.

2238

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 61, NO. 8, AUGUST 2014

on actual displacements). To quantify how well the linearized fit approximated the nonlinear torque data in that direction, we calculated the square of the Pearson correlation coefficient R2 for each movement direction. Although this directional stiffness was estimated as a linear fit in each direction, over the entire range it is a nonlinear estimate because stiffness is not the same in opposite directions. Finally, we also computed a linear estimate of stiffness about the origin by applying multiple linear regression to the entire range (as opposed to each direction), producing a single 2 × 2 stiffness matrix [2]. We performed Jarque–Bera tests (MATLAB’s jbtest function) of normality of the data and a one-way ANOVA to examine differences in stiffness between FLX, EXT, RD, and UD. Furthermore, we used Tukey’s honestly significant difference (HSD) test to find pairs that differed significantly from each other. 3) Nonlinearity of Wrist Stiffness: To evaluate the degree to which passive wrist stiffness is nonlinear, we compared our nonlinear estimate to our linear estimate by calculating for each direction the difference in stiffness, normalized by the directional stiffness in that direction Kerror = |Klinear − Knon -linear | /Knon -linear

(4)

where Knon -linear is the directional stiffness (defined above) and Klinear is the linear stiffness in each direction. To provide a fair comparison with Knon -linear (which only includes stiffness along the displacement, as described above), we derived Klinear separately for each direction from the conservative portion of the stiffness matrix as the ratio between the component of the restoring torque in the direction of displacement and the corresponding displacement. 4) Spring-Like Property of the Wrist: In general, a vector field can be decomposed according to Helmholtz’s theorem into two components: conservative and nonconservative. The conservative component is the gradient of a scalar function, and there is no work associated with a closed-loop integral within a conservative vector field. In contrast, the nonconservative component exhibits a rotational field (aka curl field), and closed-loop displacements add or remove energy from the system. Energetically passive systems, such as purely elastic springs, are perfectly conservative. However, due to the possible presence of reflex feedback coupling motions of one DOF to muscles acting on another, the passive wrist may be nonconservative, in particular if the gains of coupling feedback are unequal [17]. To quantify the extent to which our measurements of passive wrist stiffness were nonconservative, we decomposed the stiffness matrix K into its conservative component (Ks , the symmetric part of K) and nonconservative component (Ka , the antisymmetric part of K) and compared their relative magnitudes. Specifically, we calculated the ratio of the square roots of the determinants of each component   ρ = det (Ka )/ det (Ks ). (5) For a system to be fundamentally spring-like requires that the nonconservative component of the stiffness be zero, i.e., det (Ka ) = 0 and ρ = 0. A more detailed discussion of springlike behavior is provided in Appendixes B and C of [16].

To determine whether any nonconservative components were significant, we compared them to any nonconservative components in the data taken during the measurement with the mass in place of the human subject’s hand (described above). Because the mass was completely passive, any nonconservative components identified were clearly artifacts, possibly due to noise from sensors and actuators. Nonconservative stiffness components measured with a human subject that were greater than these artifactual components were regarded as noteworthy. One subject’s value of ρ was more than three standard deviations from the mean of all subjects. This subject’s data were excluded from further analysis. 5) Representations of Wrist Stiffness: The nonlinear and linear stiffness estimates were represented as follows. Nonlinear stiffness was graphed as a polar plot representing stiffness in each movement direction, the angle corresponding to the direction of displacement, and the radius corresponding to stiffness in that direction. To allow graphical comparison of our linear stiffness estimate to both our nonlinear estimate and to previous measurements of linear stiffness, we represented our linear stiffness estimate as a polar plot and as a stiffness ellipse.3 The stiffness ellipse marks the locus of the restoring torque due to a unit displacement in each direction and is characterized by its shape, orientation, and area [18]. The lengths and directions of the semimajor axes of the stiffness ellipse are given by the eigenvalues (λ1 and λ2 , with λ1 ≥ λ2 ) and eigenvectors of the associated stiffness matrix, respectively. We represented the shape of the ellipse by the ratio of the semimajor axis lengths (computed as λ1 /λ2 ), the orientation by the angle between the major axis (given by the eigenvector associated with λ1 ), and the RUD axis, and we computed the area of the ellipse as πλ1 λ2 . III. RESULTS Raw displacement, torque, and EMG data for a representative movement involving FE and RUD are shown in Fig. 2. Movement in PS was minimal, with a total range of less than 2◦ for all subjects (not shown). A. EMG Analysis In general, subjects were able to maintain their wrist muscles in a relaxed condition during the experiment: when averaged over all subjects, the activation levels of flexors (FCR and FCU) and extensors (ECR and ECU) were less than 1.4 and 3.2% of MVC, respectively (see Table I). There was no significant difference between this muscle activity and the estimated noise levels (p > 0.05 for all muscles; Table I). B. Stiffness Estimation A representative vector field approximation for a single repetition of a subject’s outbound data is shown in Fig. 3. This is an 3 It is common to only plot the conservative portion of the stiffness matrix because only that portion is guaranteed to have real eigenvalues, from which the shape, orientation, and area of the ellipse are calculated. Because the nonconservative portion of passive systems is small in general (it would be zero in the absence of experimental error) and in our experiment (see Results), the ellipses of the total and the conservative portion of the stiffness matrix are similar.

PANDO et al.: POSITION-DEPENDENT CHARACTERIZATION OF PASSIVE WRIST STIFFNESS

2239

Fig. 2. Raw displacement, torque, and EMG data, all from the same representative movement involving both FE and RUD. (a) Displacements that occurred during the measurement with the subject are shown in black, whereas the displacements without the subject (with a mass in place of the subject’s hand) are in gray. The data preceding the circles were considered to be in the range of short-range stiffness and were excluded from further analysis. (b) The torques that produced the displacements shown in A. Negative torque in extension and radial deviation represent torque in flexion and ulnar deviation, respectively. Black and gray colors indicate measurements with and without the subject, as in A. (c) Raw EMG data in each of the four measured wrist muscles (distributed vertically for visualization). The bar on the bottom right shows 10% MVC. (d) Detrended, rectified, and low-pass filtered EMG signals, as well as their averages over the duration of a single movement. The order of muscles (top to bottom) is the same as in c.

TABLE I MUSCLE ACTIVITY AND ESTIMATED NOISE LEVEL, EXPRESSED AS A PERCENTAGE OF MVC

extremely precise characterization: the mean deviation between the “raw” torque measurements and the scalar function estimates was less than 0.0005 Nm for both τFE and τRUD in all subjects, which is about 200 times smaller than the measurement error of the apparatus (about ±0.1 Nm). This precision emerges from the averaging provided by the thin-plate spline smoothing we employed. The five repetitions were approximated separately for both outbound and inbound data. No statistically significant difference was found between them. Averaged over all subjects, p = 0.70 ± 0.07, 0.67 ± 0.08, 0.65 ± 0.09, and 0.60 ± 0.09 (mean ± SE) for φ1 outbound, φ1 inbound, φ2 outbound, and φ2 inbound, respectively (with all subjects satisfying p > 0.05 in all cases and power ≈ 1.0), confirming the repeatability of our measurements.

We averaged the repetitions into a single vector field which provides the most complete representation of passive wrist stiffness, and we identified linear and nonlinear (directional) stiffness estimates from this single vector field. Clear directional variation of wrist stiffness was observed in FE–RUD space (F3 , 52 = 27.29, p ≈ 0.0 when KRD , KUD , KFLX , and KFLX were compared): stiffness in the FE direction was consistently lower than in the RUD direction for both nonlinear and linear stiffness estimates, resulting in an anisotropic shape [see Fig. 4(a)]. In addition, the nonlinear stiffness estimate showed differences in FLX versus EXT and RD versus UD which were qualitatively similar across subjects: 11 out of 14 subjects showed a general trend of KRD > KUD > KEXT > KFLX , while two subjects showed a pattern of KRD > KUD > KFLX > KEXT and 1 subject the pattern of KUD > KRD > KEXT > KFLX . Considering all subjects together, KRD was significantly higher than KUD (F1,26 = 6.26, p = 0.019), while KFLX and KEXT were not significantly different (F1,26 = 2.36, p = 0.137), that is KRD > KUD > KEXT ≈ KFLX . Directional stiffness values, estimated as separate linear fits in each movement direction, with their corresponding R2 values, are given in Table II. High R2 values in all movement directions demonstrated that estimating stiffness for each movement direction as a linear fit was a good approximation of the nonlinear field. The linear approximation of the entire CVFA is presented as the conservative portion of the stiffness matrix and as the

2240

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 61, NO. 8, AUGUST 2014

IV. DISCUSSION Because stiffness effects dominate the dynamics of wrist rotations, understanding how the neuromuscular system controls wrist rotations requires a knowledge of wrist stiffness throughout the ROM of the wrist. Here, we provide an extremely precise characterization of relaxed wrist stiffness over the central 42% of the ROM of the wrist, which covers approximately 70% of wrist activity during activities of daily living [6]. While our findings agree qualitatively with a previous study, which covered 14% of wrist ROM [2], we found stiffness over this three-times greater range to be significantly more anisotropic but less misaligned with the anatomical axes, and we provide both linear and nonlinear estimates of passive wrist stiffness over this greater range. A. Comparison to Previous Studies

Fig. 3: Representative continuous vector field approximation for a single subject, shown as (a) the torque in extension versus extension and radial deviation and (b) the torque in radial deviation versus extension and radial deviation (negative extension and radial deviation represent flexion and ulnar deviation, respectively). The red data points are the torque compensated measurements, which were approximated by the surfaces representing scalar functions (φ 1 ,φ 2 ) obtained by thin-plate spline smoothing with generalized cross validation. EXT and RD represent extension and radial deviation, respectively.

anisotropy, angle, and area of the corresponding stiffness ellipse (see Table III), and displayed as a polar plot [see Fig. 4(b)] and as a stiffness ellipse [see Fig. 4(c)]. We compared the linear stiffness estimate to the nonlinear stiffness estimate using Kerror , which represents the normalized difference between the polar plots of the nonlinear stiffness [see Fig. 4(a)] and linear stiffness [see Fig. 4(b)] in each direction. Averaged over all directions, Kerror was 0.33 ± 0.04 (mean ± SE of all subjects), 0.20 ± 0.03, and 0.25 ± 0.03 for outbound, inbound, and combined data, respectively. Errors substantially increased when only the half of the data points nearest to the limit of the ROM were used: Kerror was 0.57 ± 0.11, 0.27 ± 0.04, and 0.34 ± 0.05 for outbound, inbound, and combined data, respectively. C. Spring-Like Property of the Passive Wrist A representative decomposition of one subject’s outbound torque–displacement vector field into conservative and nonconservative components is shown in Fig. 5(a). As with this subject, in general, the nonconservative component was substantially smaller than the conservative components: averaged over all subjects, ρ (4) was 0.131, 0.071, 0.102 for outbound, inbound, and combined data, respectively. Three subjects showed nonconservative components which were modest compared to conservative components but statistically significantly different from zero. An example subject’s data exhibiting nonconservative components is shown in Fig. 5(b).

In addition to providing stiffness over a range three times greater than previous studies, the current study involved several innovations designed to provide a more accurate measurement of wrist stiffness. These included 1) a custom-built attachment device that encouraged relaxation of finger and wrist muscles, 2) placement of the wrist in standard neutral position to be consistent between subjects and allow comparison with future studies, 3) removal of short-range stiffness effects, 4) compensation for gravity and robot dynamics on a subject-by-subject basis, and 5) a nonlinear multi-DOF data fit using thin-plate spline smoothing optimized with generalized cross validation. To relate our 2-DOF stiffness measurement to prior studies, which are almost all limited to a single DOF, we compared the diagonal elements of our stiffness matrix (which represent stiffness in pure FE and RUD) to the prior measurements. Note that our matrix elements were obtained through regression over the entire vector field (in both DOF) and therefore represent stiffness over a much larger range than 1-DOF measurements. Previous 1-DOF measurements ranged from less than 0.15 Nm/rad [19] and 0.32–0.7 Nm/rad [20] to 0.88 Nm/rad [2], 1 Nm/rad [11], 2.2 Nm/rad [21], and 3 Nm/rad [22]. The 20-fold range of these previous measurements likely reflects different experimental procedures; nevertheless our average stiffness in FE (0.84 Nm/rad) falls in the middle of it. Our measurement of stiffness in RUD (see K22 = 2.02 Nm/rad in Table III) is higher than the two previous single-DOF measurements of 1.45 [23] and 1.63 Nm/rad [2]. Since stiffness increases toward the limits of the wrist’s ROM, the increased RUD stiffness seen in our study is likely due to the larger range over which we measured stiffness (the range over which we measured FE stiffness was also larger than in prior studies but did not approach the limit of the wrist’s ROM in FE). To relate the conservative portion of the stiffness matrix (including off-diagonal elements) to prior work, we compared our measurement of stiffness to the only previously published measurement of coupled wrist stiffness [2] and found similarities and differences. We compared the linear stiffness estimates from the two studies by comparing the shape, orientation, and size of the stiffness ellipses (see Table III). While our findings are qualitatively similar (the shape is anisotropic, with greater stiffness in RUD than in FE; the orientation of the ellipse is slightly

PANDO et al.: POSITION-DEPENDENT CHARACTERIZATION OF PASSIVE WRIST STIFFNESS

2241

Fig. 4. Stiffness measurements averaged over all subjects (left column) and for individual subjects (right column). (a) Polar plots of directional stiffness. The angle corresponds to the direction of displacement, and the radius corresponds to stiffness in that direction (see Methods for details). (b) Polar plots of the (linear) stiffness matrix. (c) Stiffness ellipses of the conservative portion of the linear stiffness, which mark the locus of restoring torque due to unit displacements in all directions. Negative torque in extension (EXT) and radial deviation (RD) represent torque in flexion and ulnar deviation, respectively. In all plots in the left column, the thick solid line represents the mean stiffness, and the gray bands represent the mean ± 1 standard error. The solid thin line and the dotted thin line in a represent the mean stiffness for outbound and inbound movements, respectively. The outbound and inbound means are also included in b and c but are barely distinguishable from the overall mean. All plots in the right column depict the average of outbound and inbound stiffness for each subject.

tilted counterclockwise relative to the anatomical axes; and the size of ellipse is highly variable among subjects), we found that our stiffness measurements were significantly more anisotropic (p < 0.001; eigenvalue ratio of 2.69 ± 0.52 versus 1.58 ± 0.39) and less tilted (p < 0.001; 12.1 ± 4.6◦ versus 21.2 ± 9.2◦ ), with a marginally significant difference in the area of the stiffness ellipse (p = 0.054; 5.61 ± 3.42 versus 7.19 ± 3.41).

We believe that these differences may be due to the following (deliberate) methodological differences. First, the current measurement covers a range that is three times greater. Second, in the previous study, the origin of the stiffness measurement was chosen as the position for grasping a handle, while we placed the origin in standard neutral wrist position [9] to allow for widespread comparison with future studies. Compared to

2242

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 61, NO. 8, AUGUST 2014

TABLE II DIRECTIONAL STIFFNESS AND CORRESPONDING R2 VALUES

As in the previous 2-DOF study [2], intersubject differences in the shape and orientation of wrist stiffness were small (for all subjects, stiffness is greater in RUD than FE, with the direction of greatest stiffness slightly pronated from pure RUD), while differences in the magnitude of wrist stiffness were relatively large (see Fig. 4 and Table III). This is to be expected: the shape and orientation reflect the anatomical geometry of the wrist bones and tendon routing, which is relatively stereotyped between subjects, whereas stiffness magnitude increases with the size of muscles and ligaments (as do springs arranged in parallel), which varies greatly between subjects. B. Comparison Between Nonlinear and Linear Estimates

TABLE III LINEAR STIFFNESS

The simplified linear and nonlinear estimates of stiffness in this study were compared to evaluate the degree to which passive wrist stiffness is nonlinear over this range. While both methods revealed that wrist stiffness is anisotropic and pronated relative to the anatomical axes, the estimation error (measured as Kerror ) was substantial for both outbound, inbound, and combined data (higher than 30%, 20%, and 25%, respectively). Furthermore, the finding that this error increased substantially toward the limit of the ROM of the wrist implies that the accuracy of a linear representation of stiffness (from multiple linear regression) is best for smaller ranges and decreases as the size of the movement range increases. Furthermore, a linear analysis is unable to discern differences in stiffness in opposite directions, whereas the nonlinear estimate revealed substantial differences KRD > KUD > KEXT > KFLX ). C. Spring-Like Behavior

our study, subjects in the previous study were in approximately 20◦ of EXT and 2◦ of UD, with the anatomical axes pronated approximately 11◦ relative to the robot. Third, in the previous study, the fingers were flexed around a handle. In our study, we deliberately chose to leave the fingers unconstrained to allow the wrist to be in standard neutral position and to encourage finger and wrist muscles to relax fully. Fourth, the previous study accounted for gravity and robot dynamics through an average measurement of these effects, while we measured and removed these effects separately for each subject.

An ideal mechanical spring is energetically passive. Comparison of the conservative and nonconservative components of the torque field provided a quantification of the energetic passivity of the passive wrist joint. Most subjects’ nonconservative component was statistically indistinguishable from zero, indicating that the wrist of young, healthy subjects is predominantly springlike in a relaxed condition. Identification of energetic passivity is important, since it ensures stable interaction with a dynamically passive environment [24]. Although the nonconservative component was small in all subjects, it was statistically nonzero in three subjects. Unlike the linear approximation, the CVFA allows one to determine the location (in the FE-RUD space) where the nonconservative component is nonzero [see Fig. 5(b)]. As observed previously in the human ankle, this information has clinical relevance because significant nonconservative stiffness in a specific region implies that intermuscular feedback between muscles acting on different DOFs is unbalanced in that region [25]. D. Limitations of This Study The stiffness reported here is only valid for the range over which it was measured [see black circles in Fig. 1(a)]. The robot targets (open circles) were chosen to maximize the measurement range within the constraints of the ROM of the wrist and robot (solid line and dashed rectangle, respectively; the slight

PANDO et al.: POSITION-DEPENDENT CHARACTERIZATION OF PASSIVE WRIST STIFFNESS

2243

Fig. 5. Comparison of conservative (left column) and non-conservative components (right column) of the restoring torque. Torque is represented by an arrow (units shown in the center of the figure) with its tail at the tip of the displacement vector. The axes in all subplots are the same as in the bottom left subplot. (a) Representative subject. (b) Subject showing non-conservative components (non-zero curl). Statistically significant non-zero curl is highlighted with a box.

disagreement between the targets and the overlap of the wrist and robot ROM was caused by differences in neutral position between subjects). However, in practice, the robot did not reach the targets (compare black and open circles) because of the steady-state error of the proportional-derivative controller and the torque limit of the robot motors (see Methods). This resulted in an anisotropic measurement range and, by EXT, measurements closer to the limit of the wrist ROM in some directions than in others. Because stiffness generally increases as one approaches the limit of the ROM, stiffness measurements can be expected to be higher in the directions in which the measurement range approached more closely the limit of the ROM. Indeed, our finding that KRD > KUD > KEXT ≈ KFLX matches the percentage of the ROM measured [see Fig. 1(a)]. However, this simply reflects the fact that stiffness and ROM are two sides of the same phenomenon: stiffness is greater in RD than in UD, FLX, and EXT—even when the nominal measurement range was isotropic [2]—at least in part because its ROM is smaller than in UD, FLX, and EXT. In other words, given that the ROM is the locus at which the stiffness torque reaches the same, high value in each direction (an “isotorque curve”), the order of stiffness (KRD > KUD > KEXT ≈ KFLX ) is the opposite of the order of ROM (ROMRD < ROMUD < ROMEXT ≈ ROMFLX [5]). This is true for the small, isotropic measurement range covering 14% of the wrist ROM in [2],4 the larger, anisotropic measurement range covering 42% of the wrist ROM in the current study, and would likely be true for the entire wrist ROM as well. 4 Formica et al. found K R D > K U D > K E X T > K F L X (i.e., a difference between K E X T and K F L X ), but the origin was placed in 20◦ of extension relative to neutral position, so the order of ROM was ROMR D < ROMU D < ROME X T < ROMF L X (i.e., there was a difference between ROME X T and ROMF L X ), validating the statement that the orders of stiffness and ROM are inversely related.

Likewise, the reported stiffness is, strictly speaking, only valid for the conditions under which it was measured. These conditions include a relaxed-muscle state, pseudostatic conditions (average speed was 4.3◦ /s), unconstrained fingers, and straight movements to and from a neutral posture. That said, to the extent that we measured static behavior (which was the reason for the procedures used), the stiffness depends only on the current location and not on the starting location or direction of the movement (ignoring short-range stiffness and hysteresis). Also, studies of shoulder–elbow movements suggest that while muscle activity clearly increases stiffness magnitude, the change in the shape and orientation of the stiffness ellipse is quite limited [26], [27]. Accordingly, relaxed wrist stiffness measured under conditions similar to those in this study [2] has successfully been used to explain wrist behavior under different conditions, including fast movements [4]. Finally, the displacement and torque applied to the wrist joint were not measured directly at the wrist joint but estimated from other measured variables. Care was taken to align the axes of the wrist and robot as closely as possible, but the alignment was not perfect. In particular, the location of the robot’s differential gear mechanism [8] forces the FE axis of the robot to be approximately 2 cm proximal to the FE axis of the wrist [see Fig. 1(b)–(c)]. Also, while the robot is designed to have an offset between the RUD axes of the robot and the wrist, this offset causes a slightly nonlinear relationship between the RUD angles of the wrist and the robot [10]. Further, the torque applied to the wrist was assumed to be equal to the torque at the robot end-effector, which was estimated from the voltages sent to the servo-amplifiers [10]. The combined effect of these limitations (which have been present in all prior 2-DOF measurements of wrist stiffness) on the estimate of wrist stiffness is believed to

2244

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 61, NO. 8, AUGUST 2014

be relatively small, but further experimentation and analysis are required to quantify it. E. Relevance to Motor Control and Stiffness Pathologies The anisotropic nature of wrist stiffness causes path curvature [3] unless the neuromuscular system specifically compensates for the anisotropy [4]. Our stiffness estimates are more anisotropic than previous measurements, suggesting that stiffness anisotropy may be an even greater control challenge than previously recognized. Another consequence of this anisotropy is that rotating the wrist in some directions requires less torque than in others, providing the opportunity for the neuromuscular system to choose “paths of least resistance” when it has the choice. It was previously observed [2] that the direction of least stiffness (from EXT and slight RD to FLX and slight UD) corresponds to the direction of the “dart thrower’s motion,” a commonly used wrist movement direction [7], [28]. The nonlinear estimate shows a dramatic decrease in stiffness in this direction, even compared to directions that are close-by [seen as a “pinching” of the stiffness trace in Fig. 4(a)], indicating an even greater anisotropy than indicated by the stiffness ellipse [Fig. 4(c)]. Therefore, following this “path of least resistance” would be more advantageous than previously believed. Abnormal joint viscoelasticity is one of the most common symptoms of neuromuscular and biomechanical disorders and demands more consistent and quantitative analysis. Here we provide a quantitative characterization of relaxed healthy wrist stiffness over the central 42% of the wrist ROM and about 70% of wrist activity during activities of daily living. Following recent investigations of ankle stiffness [29]–[31], the wrist stiffness measurements presented here can serve as a baseline for evaluating stiffness pathologies, tracking changes with injury and rehabilitation, and assessing the effectiveness of orthotic, implant, and arthroplasty designs. REFERENCES [1] S. K. Charles and N. Hogan, “Dynamics of wrist rotations,” J. Biomechanics, vol. 44, pp. 614–621, 2011. [2] D. Formica, S. K. Charles, L. Zollo, E. Guglielmelli, N. Hogan, and H. I. Krebs, “The Passive stiffness of the wrist and forearm,” J. Neurophysiol., vol. 108, pp. 1158–1166, 2012. [3] S. K. Charles and N. Hogan, “The curvature and variability of wrist and arm movements,” Exp. Brain Res., vol. 203, pp. 63–73, May 2010. [4] S. K. Charles and N. Hogan, “Stiffness, not inertial coupling, determines path curvature of wrist motions,” J. Neurophysiol., vol. 107, pp. 1230– 1240, 2012. [5] K. N. An, R. A. Berger, and W. P. I. Cooney, Biomechanics of the Wrist Joint. New York, NY, USA: Springer-Verlag, 1991. [6] J. Ryu, W. P. Cooney, L. J. Askew, K. N. An, and E. Y. S. Chao, “Functional ranges of motion of the wrist joint,” J. Hand Surgery-Amer. Vol., vol. 16A, pp. 409–419, May 1991. [7] A. K. Palmer, F. W. Werner, D. Murphy, and R. Glisson, “Functional wrist motion—A biomechanical study,” J. Hand Surgery-Amer. Volume, vol. 10 A, pp. 39–46, 1985. [8] H. Krebs, B. T. Volpe, D. Williams, J. Celestino, S. Charles, D. Lynch, and N. Hogan, “Robot-aided neurorehabilitation: A robot for the wrist rehabilitation,” IEEE Trans. Neural Syst. Rehabil. Eng., vol. 15, no. 3, pp. 327–335, Sep. 2007.

[9] G. Wu, F. C. T. Van Der Helm, H. E. J. Veeger, M. Makhsous, P. Van Roy, C. Anglin, J. Nagels, A. R Karduna, K. McQuade, X. Wang, F. W. Werner, and B. Buchholz, “ISB recommendation on definitions of joint coordinate systems of various joints for the reporting of human joint motion—Part II: Shoulder, elbow, wrist and hand,” J. Biomechanics, vol. 38, pp. 981–992, May 2005. [10] J. R. Celestino, “Characterization and control of a robot for wrist rehabilitation,” M.S. Thesis, Massachusetts Institute of Technology, Cambridge, MA, USA, 2003. [11] H. W. Axelson and K. E. Hagbarth, “Human motor control consequences of thixotropic changes in muscular short-range stiffness,” J. Physiol.London, vol. 535, pp. 279–288, Aug. 2001. [12] K. S. Campbell, “Short-range mechanical properties of skeletal and cardiac muscles,” in Muscle Biophysics: From Molecules to Cells (Advance in Experimental Medicine and Biology Series), D. E. Rassier, Ed., New York: Springer 2010, vol. 682, pp. 223–246. [13] P. De Leva, “Adjustments to Zatsiorsky-Seluyanov’s segment inertia parameters,” J. Biomech., vol. 29, pp. 1223–1230, Sep. 1996. [14] F. L. Bookstein, “Principal warps - thin-plate splines and the decomposition of deformations,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 11, no. 6, pp. 567–585, Jun. 1989. [15] G. Wahba, Spline Models for Observational Data. Philadelphia, PA, USA: SIAM, 1990. [16] H. Lee, P. Ho, M. A. Rastgaar, H. I. Krebs, and N. Hogan, “Multivariable static ankle mechanical impedance with relaxed muscles,” J. Biomech., vol. 44, pp. 1901–1908, Jul. 2011. [17] N. Hogan, “The mechanics of multi-joint posture and movement control,” Biological Cybern., vol. 52, pp. 315–331, 1985. [18] J. M. Dolan, M. B. Friedman, and M. L. Nagurka, “Dynamic and loaded impedance components in the maintenance of human arm posture,” IEEE Trans. Syst. Man Cybern., vol. 23, no. 3, pp. 698–709, May./Jun. 1993. [19] S. L. Lehman and B. M. Calhoun, “An identified model for human wrist movements,” Exp. Brain Res., vol. 81, pp. 199–208, 1990. [20] C. Gielen and J. C. Houk, “Nonlinear viscosity of human wrist,” J. Neurophysiol., vol. 52, pp. 553–569, 1984. [21] A. B. Leger and T. E. Milner, “Passive and active wrist joint stiffness following eccentric exercise,” Eur. J. Appl. Physiol., vol. 82, pp. 472–479, Aug. 2000. [22] S. J. De Serres and T. E. Milner, “Wrist muscle activation patterns and stiffness associated with stable and unstable mechanical loads,” Exp. Brain Res., vol. 86, pp. 451–458, 1991. [23] N. Rijnveld and H. I. Krebs, “Passive wrist joint impedance in flexionextension and abduction-adduction,” presented at the Int. Conf. Rehabil. Robot., Noordwijk, The Netherlands, 2007. [24] N. Hogan, “On the stability of manipulators performing contact tasks,” IEEE J. Robot. Autom., vol. 4, no. 6, pp. 677–686, Dec. 1988. [25] H. Lee, T. Patterson, J. Ahn, D. Klenk, A. Lo, H. I. Krebs, and N. Hogan, “Static ankle impedance in stroke and multiple sclerosis: A feasibility study,” in Proc. IEEE Annu. Int. Conf. Eng. Med. Biol. Soc., 2011, pp. 8523–8526. [26] F. A. Mussa-Ivaldi, N. Hogan, and E. Bizzi, “Neural, mechanical, and geometric factors subserving arm posture in humans,” J. Neuroscience, vol. 5, pp. 2732–2743, 1985. [27] E. J. Perreault, R. F. Kirsch, and P. E. Crago, “Voluntary control of static endpoint stiffness during force regulation tasks,” J. Neurophysiol., vol. 87, pp. 2808–2816, Jun. 2002. [28] Z. M. Li, L. Kuxhaus, J. A. Fisk, and T. H. Christophel, “Coupling between wrist flexion-extension and radial-ulnar deviation,” Clinical Biomech., vol. 20, pp. 177–183, Feb. 2005. [29] A. Roy, L. W. Forrester, R. F. Macko, and H. I. Krebs, “Changes in passive ankle stiffness and its effects on gait function in people with chronic stroke,” J. Rehabil. Res. Develop., vol. 50, pp. 555–571, 2013. [30] L. W. Forrester, A. Roy, H. I. Krebs, and R. F. Macko, “Ankle training with a robotic device improves hemiparetic gait after a stroke,” Neurorehabil. Neural Repair, vol. 25, pp. 369–377, 2011. [31] A. Roy, H. I. Krebs, C. T. Bever, L. W. Forrester, R. F. Macko, and N. Hogan, “Measurement of passive ankle stiffness in subjects with chronic hemiparesis using a novel ankle robot,” J. Neurophysiol., vol. 105, pp. 2132–2149, 2011. Authors’ photographs and biographies not available at the time of publication.