Experimental measure of arm stiffness during single reaching movements with a time-frequency analysis.

J Neurophysiol 110: 2484 –2496, 2013. First published August 14, 2013; doi:10.1152/jn.01013.2012.

Experimental measure of arm stiffness during single reaching movements with a time-frequency analysis Davide Piovesan,1,2,3 Alberto Pierobon,3 Paul DiZio,3,4 and James R. Lackner3,4 1 Sensory Motor Performance Program (SMPP), Rehabilitation Institute of Chicago, Chicago, Illinois; 2Department of Physical Medicine and Rehabilitation, Northwestern University, Chicago, Illinois; 3Ashton Graybiel Spatial Orientation Laboratory, Brandeis University, Waltham, Massachusetts; and 4Volen Center for Complex Systems, Brandeis University, Waltham, Massachusetts

Submitted 21 November 2012; accepted in final form 11 August 2013

joint stiffness measure; time-frequency analysis; reflex modulation

of the human arm can be represented as a linear combination of three force components that are dependent on acceleration, velocity, and position. The proportionality coefficients of such components are inertia, damping, and stiffness, respectively. Inertia depends on the mass distribution and geometry of the limb segments. The inertia of each segment about its center of gravity is invariant with respect to joint rotation. By contrast, stiffness and damping can be modulated by neural activity during movement (Mussa-Ivaldi et al. 1985; Nichols and Houk 1976; Perreault et al. 2001; Tsuji et al. 1995). Hence, human motor control studies often focus on the estimation of joint stiffness and damping (Hondori and Tech 2011; Hondori et al. 2012; Masia et al. 2011). Joint stiffness is mostly influenced by muscle activity, segmental reflexes, and tissue biomechanics (Osu and Gomi 1999; Perreault et al. 2002). While many experiments have studied stiffness during postural tasks, the modulation of stiffness during unconstrained movements is still poorly understood and its measurement has produced contrasting results, even under seemingly similar experimental conditions (Burdet et al. 2000;

THE RESISTANCE TO MOTION

Address for reprint requests and other correspondence: D. Piovesan, Sensory Motor Performance Program (SMPP), Rehabilitation Institute of Chicago, 345 E. Superior, Suite 1396, Chicago, IL 60611-2654 (e-mail: [email protected]). 2484

Darainy et al. 2007; Frolov et al. 2006; Gomi and Kawato 1997; Mah 2001). The development of the virtual trajectory theory (Hogan 1984) led to the prediction that stiffness would be higher during ongoing movements and less during static postural tasks (Flash 1987). Indeed, Won and Hogan (1995) found that significant restoring forces were generated during perturbed movements to maintain stability, supporting the hypothesis that high joint stiffness is used to control movements. However, measures of stiffness made during movements (Burdet et al. 2000; Gomi and Kawato 1997) also provided values lower than those obtained during postural tasks such as supporting the limb against gravity (Darainy et al. 2004; Mussa-Ivaldi et al. 1985; Tsuji et al. 1995). Segmental reflexes play a crucial role in the modulation of stiffness: as reflex gain increases, so does joint stiffness (Houk 1979). Reflexive activity is inhibited or reduced during voluntary movement (Bawa and Sinkjaer 1999; Seki et al. 2003), suggesting that stiffness might be lower in active motor tasks than in postural tasks. Joint stiffness also depends on the mechanical properties of muscles. Displacements introduced by dynamic estimation procedures, whether directly or as a consequence of force perturbations, usually modify the trajectory beyond the natural variability of unperturbed baseline movements. The elastic forces generated at the joints are assumed by these techniques to vary linearly across the range of displacement. However, Kearney and Hunter (1982) and Loram and colleagues (2007a, 2007b, 2009) demonstrated that a strong nonlinearity exists in the elastic force field around an unperturbed movement trajectory, which may result in measured stiffness appearing higher for small displacements. This factor could partly explain the incongruence between the low stiffness reported in dynamic estimation studies, where perturbations were ⬃1 cm, and the high stiffness predicted by virtual trajectory theories, which implies very small displacements. Many techniques have been proposed to estimate stiffness and damping during postural tasks. By contrast, methods for assessing stiffness during limb motion are less common, mainly because of the technical difficulties associated with requiring a repeatable task and the long duration of the experiments. Stiffness is often estimated by using perturbations to elicit measurable deviations from baseline trajectories (Burdet et al. 2000; Darainy et al. 2007; Gomi and Kawato 1997). Each resulting trajectory reflects the displacement induced on the trajectory of the normally unperturbed voluntary movement. The voluntary and induced components of a perturbed trajectory are difficult to separate in the time domain, especially when many move-

0022-3077/13 Copyright © 2013 the American Physiological Society

www.jn.org

Downloaded from http://jn.physiology.org/ by 10.220.33.4 on November 24, 2017

Piovesan D, Pierobon A, DiZio P, Lackner JR. Experimental measure of arm stiffness during single reaching movements with a time-frequency analysis. J Neurophysiol 110: 2484 –2496, 2013. First published August 14, 2013; doi:10.1152/jn.01013.2012.—We tested an innovative method to estimate joint stiffness and damping during multijoint unfettered arm movements. The technique employs impulsive perturbations and a time-frequency analysis to estimate the arm’s mechanical properties along a reaching trajectory. Each single impulsive perturbation provides a continuous estimation on a single-reach basis, making our method ideal to investigate motor adaptation in the presence of force fields and to study the control of movement in impaired individuals with limited kinematic repeatability. In contrast with previous dynamic stiffness studies, we found that stiffness varies during movement, achieving levels higher than during static postural control. High stiffness was associated with elevated reflexive activity. We observed a decrease in stiffness and a marked reduction in long-latency reflexes around the reaching movement velocity peak. This pattern could partly explain the difference between the high stiffness reported in postural studies and the low stiffness measured in dynamic estimation studies, where perturbations are typically applied near the peak velocity point.

ARM STIFFNESS MEASURE WITH TIME-FREQUENCY ANALYSIS

MATERIALS AND METHODS

Subjects Thirteen right-handed adults (10 men, 3 women; age: 32 ⫾ 14 yr, mass: 78 ⫾ 15 kg, height: 1.75 ⫾ 0.08 m) participated after giving informed consent to a protocol approved by the Brandeis Human Subjects Committee. Protocol The subject sat on an adjustable stool at a custom-built table with a start position and a target position specified on its surface as described below (Fig. 1). The table was designed to ensure that its lowest resonant frequency (30 Hz) was above the bandwidth of interest for the stiffness estimations. The height of the table surface was individually leveled to the sixth rib of each subject in order to keep reaching trajectories on a horizontal plane passing through the subject’s gleno-humeral joint (see Fig. 1). Forward reaches of the subject’s right arm were aligned with the median plane, and contact of the arm with the table was allowed only at the starting position and the target location. The planarity of the movements and displacement of the shoulder centroid for this type of movement are reported in Piovesan et al. (2011c). The y-axis of the Cartesian reference frame was aligned with the starting and target positions, oriented in the direction of the movement; the x-axis pointed laterally and the z-axis vertically in a right-handed reference frame. The angular configuration of the reaching arm was kept consistent across subjects. Considering both shoulder and elbow angles to be at

Fig. 1. Experimental setup. The reaching distance d was adjusted for each subject to normalize the angular configuration of the joints at the start and target positions across subjects. After lifting their elbow to support their arm against gravity, subjects reached to the target in one single fluid motion. The trajectory was linear, parallel to the reference system’s y-axis. The x-axis is oriented laterally. The z-axis is oriented up, in a right-hand reference system. The position of the arm was recorded with active markers numbered from 1 to 8. A light manipulandum randomly perturbed the hand, applying low-magnitude forces along directions A, B, and C at points 1/4d, 1/2d, and 3/4d along the trajectory, within a tolerance of ⫾1%.

0° with the arm completely extended laterally, the shoulder and elbow joint angles were 35 ⫾ 4° and 115 ⫾ 6° with the index finger at the starting point and 60 ⫾ 6° and 55 ⫾ 8° at the target, respectively. The joint angles were computed from the positions of a set of active optical markers projected onto the horizontal plane (Fig. 1). Subjects were instructed to place their finger at the starting point and selfsupport the weight of their arm against gravity prior to the initiation of each reaching movement. Surface electromyographic (SEMG) activity was recorded and monitored in real time on an external video screen to ensure a consistent initial motor drive across subjects. This baseline SEMG activity was used as reference to normalize the SEMG signal recorded during each movement. The experiment was performed in complete darkness to eliminate any visual feedback. Two colored light-emitting diodes (LEDs) visible through a Plexiglas panel on the table marked the start and target positions. They were bright enough to detect though insufficient to illuminate the subject’s hand during the movement. Subjects were instructed to keep the tip of their index finger at the illuminated starting position, wait until the target position lighted up, and then reach in a single motion. Verbal feedback was given to the subject to ensure correct movement timing (0.500 ⫾ 0.035 s, calculated between 10% and 90% of the maximum joint angular displacements). If subjects reported any surface contact before the end of the reaching movement, that trial was discarded. After the completion of each movement, the target LED was turned off and the subject returned to the starting position to await the onset of the next trial. The subject wore a wrist guard rigid enough to transmit mechanical perturbations to the limb. It constrained wrist flexion and offered a solid attaching point near the styloid process for mounting a light (25 g) instrumented aluminum scaffold. A low-mass PHANToM robot (SensAble, Wilmington, MA) was connected to the scaffold through a gimbal (see Fig. 1). Each subject completed 12 sets of 15 movements. Within each set, six movements were unperturbed and nine were perturbed by a randomized force pulse delivered by the robot to the limb. Perturbations were applied along the trajectory at three distinct points, located at one-fourth, one-half, and three-fourths of the total distance, d, between the starting and ending points measured along the y-axis. For

J Neurophysiol • doi:10.1152/jn.01013.2012 • www.jn.org


ment repetitions are required and movement repeatability is compromised by fatigue or motor noise (Darainy et al. 2007). Some variability can be eliminated by using autoregressive models or lookup tables (Burdet et al. 2000; Darainy et al. 2007); however, multiple trials are still necessary to estimate the most probable unperturbed trajectory, as well as the stiffness at each estimation location along the trajectory. We have developed a minimally invasive technique for single-trial stiffness estimation using a single brief perturbation during a reaching movement. Our method is based on an analysis of the vibrational energy of the moving limb as a function of time and frequency (Piovesan et al. 2009, 2012c). Using a short-time Fourier transform (STFT) spectrogram and its reassignment, we can identify the spectrum of the excited natural frequencies. This allows us to measure stiffness and damping in both Cartesian and joint space with a standard modal analysis, without the need for a baseline trajectory. We demonstrate here that our method uses perturbations small enough to avoid interference with the course of the voluntary movement, nonetheless providing valid dynamical estimates of the stiffness associated with displacements within the natural range of trajectory variability. Stiffness measured during some segments of a movement were indeed higher than during posture, in contrast with the aforementioned dynamic stiffness studies, and were consistently associated with higher reflexive activity. Our technique does not require the assumption of stationarity or repeatability of motor execution. Hence, arm mechanics can be estimated on a movement-by-movement basis, thereby eliminating the need to establish a baseline trajectory. Estimation of stiffness and damping requires computation of arm inertia. We used eight regression equations and a direct measure based on water displacement to assess the influence of body segment parameters (BSPs) on our estimations (Piovesan et al. 2011c).

2485

2486


each perturbed trial, a single 0.050-s force pulse was applied at one of the aforementioned points, choosing from a set of three possible directions and three associated magnitudes (A ⫽ 3 N at 165°, B ⫽ 4 N at 45°, and C ⫽ 5 N at 285°, with 0° aligned to the x-axis). The starting and ending angular positions of the arm were the same across subjects, and the relationship between hand displacement and joint angles along the trajectory depended upon the mechanical configuration of the limb; therefore each perturbation was applied at the same joint angle configuration. Instrumentation

→

→

¨ , X˙, X, u, t)→ M(X, t) ⭸ ¨ X (t) ⫹ C(X˙, X, u, t) ⭸ X˙ (t) ⫹ K(X ⭸ X (t) ⫽ 0 (1) →

where ⭸X is the vector of deviation from a planned trajectory whose components are x and y on the plane of movement, the two degrees of → freedom (DOF) of the system in the Cartesian reference frame. ⭸X˙ →

¨ are the first and second derivative of the position vector with and ⭸X respect to time, and u is the vector of muscle activations. M, C, and K are the Cartesian representations of the matrices of inertia, damping, and stiffness, respectively. The eigenvalues and eigenvectors of Eq. 1 are the squared natural frequencies of the system and its vibrational modes, respectively. To solve the eigen-problem, Eq. 1 is decoupled and normalized through the following steps (Inman 1989; Piovesan et al. 2012c): 1

Sm ⫽ [M(X)]⫺ 2 P where matrix Sm is such that T Sm · M · Sm ⫽ In T Sm · (C) · Sm ⫽ diag [2⌫ j(t)]

BSP Estimations

(2)

(3)

T · (K) · Sm ⫽ diag [␩2j (t)] Sm

Our stiffness measurement method relies on estimates of BSPs, so it is critical to have an adequate model of the inertial characteristics of the arm. We used eight methods from the literature that provide the inertial parameters by means of regression equations of limb geometrical measures and a method based on water displacement. The methods used different techniques, including geometrical estimation [Hanavan 1964 (HV)], cadaver-based regressions [Dempster 1955 (DE); Chandler et al. 1975 (CH); Clauser et al. 1969 (CL)], in vivo mass scans [Zatsiorsky and Seluyanov 1983 (Z1); Zatsiorsky 2002 (Z2); de Leva 1996 (DL)], and photogrammetry [McConville et al. 1980 (MC)]. With the water displacement method, BSPs were calculated from the geometry of the arm and the volume of each of its subsections [Piovesan et al. 2011c (PI)]. The weight of each arm component was computed assuming the density to be uniform and the percentage distribution of various tissue types as a given (Clarys et al. 1986a, 1986b). The moments of inertia were computed about the rotation axes with the parallel axis theorem. Biomechanical Model and Data Processing

1

␦Y ⫽ M 2 · ␦X P is the eigenvector matrix, In is the identity matrix, and ␩2j (t) and ⌫j(t) are the natural frequency and modal damping ratio of the jth mode (Piovesan et al. 2012c). Equation 1 can be recast as the following homogeneous diagonal system:

␦¨ Y ⫹ diag [2⌫ j(t)]␦Y˙ ⫹ diag [␩2j (t)]␦Y ⫽ 0

(4)

Our technique estimates the solution of Eq. 4 based on the analysis of the system’s response to an impulsive perturbation in the timefrequency domain. At each instant t, the spectrogram of the impulse response of each DOF i presents as many amplitude peaks Aj ⫽ S(␻j, t) (local maxima) as the total number of DOFs. Each maximum occurs at frequency ␻j(t), which represents the instantaneous resonant frequency of the jth mode. Natural frequency and resonant frequency are variables that coincide only when the system is undamped. The natural frequency ␩j(t) is a function of the resonant frequency ␻j(t) and damping ratio ⌫j(t). It can be shown that (Galleani and Cohen 2004)

The model. In the following section, we provide a brief description and justification of our stiffness estimation technique; a complete theoretical and procedural analysis is presented in Piovesan et al. J Neurophysiol • doi:10.1152/jn.01013.2012 • www.jn.org

⌫j ⫽ ⫺␣j ⫺

␻˙ j 2␻ j

(5)


The force applied to the subjects’ wrist was measured by a six-axis load cell (Nano17—SI-25– 0.25, ATI Industrial Automation, Apex, NC). The transducer was located between the gimbal and the wrist guard, on the aluminum scaffold. Three single-axis accelerometers (Kistler 8352A-10M3) were also mounted on the scaffold and oriented along a right-hand reference frame centered and aligned with both the scaffold and the load cell. All signals were sampled at 4,000 Hz after antialiasing filtering at 1,024 Hz with a four-pole Butterworth filter. The trajectory position signals were measured with a three-sensor Optotrak motion analysis system (Northern Digital, Waterloo, ON, Canada) using eight active markers placed on the subjects’ right arm and sternum (Fig. 1), sampled at 200 Hz. The robot was positioned on top of a Kistler force platform to measure the reaction force. This allowed us to measure the force in both the mobile and stationary Cartesian reference systems, eliminating the need for coordinate transformations. SEMG was measured by four bipolar, AC-coupled surface electrodes positioned on the long head of biceps brachii (BLH), the clavicular head of pectoralis (PCH), the triceps brachii lateral head (TLH), and the deltoid posterior head (DPH), respectively. SEMG signals were collected with a wireless system (MIE, Leeds, UK), amplified (8,000⫻), low-pass filtered (2,000 Hz) for antialiasing, notch filtered (4-pole bidirectional Butterworth filter) between 59 Hz and 61 Hz to limit power line interference, and band-pass filtered (4-pole bidirectional Butterworth filter) between 35 and 500 Hz, to suppress movement artifacts and high-frequency noise. The signals were then full-wave rectified. SEMG, acceleration, and force signals were digitized with an OUDA unit (Northern Digital).

(2012c). Classical methods for estimating stiffness assume the system to be time invariant (stationary). The typical model for a timeinvariant system is a set of differential equations whose solution can be found in either the time or frequency domain. Classical Fourier transforms can be used to approach the problem in the frequency domain. When a system is stationary, the frequency content of its solution is a constant in the time domain. Our method does not require the assumption of stationarity or require the frequency content to be constant in time. The STFT provides a method for analyzing a time-varying system with the same mathematical tools applied to a stationary system by restricting observations to small temporal windows within which the stationarity assumption still applies. The signal is convoluted with a set of elementary “window” functions that are nonzero within a subset of the time domain. A spectrogram, a time-frequency representation of the magnitude of the convolution of such window functions with the signal, is computed at subsequent times. Each time point in the spectrogram represents the average of all STFTs enclosing that instant in their windows function. A second-order linear biomechanical system can be described by the following differential equation:


␩2j (t) ⫽ ␻2j ⫹ ␣2j ⫹

␴ j␻˙ j ␻j

⫺ ␣˙ j

(6)

where

␣ j(t) ⫽

d dt

共ln A ji(t)兲 ⫽

(7)

A ji(t)

n

→ sj 兺 j⫽1

⭸X⫽

Aj C˙ j Cj

⫽ ⫽

␣ j · B j · e␣ j·t B j · e␣ j·t

␣ j · B(a j ⫹ i␻ j)2 · e␣ j·t B(a j ⫹ i␻ j)2 · e␣ j·t

p j · B j · e(␣ ⫹i␻ )t 兺 j⫽1 →

j

共 H 兲T H H T CH ⫽ RH CH 共 RH 兲 '

ⱍ ⫽

n

→ p j · 共B j · e␣ ·t兲ei␻ ·t 兺 j⫽1 j

j

n

→ p j · A j · ei␻ ·t 兺 j⫽1 j

→

where pj is the eigenvector associated with the jth vibrational mode. The system is second order; therefore the frequency content of → → ¨ and ⭸X˙ are the same, as can be seen from calculating the second ⭸X derivative of Eq. 8 with respect to time: →

¨⫽ ⭸X

n

→ s¨ j 兺 j⫽1

ⱍ ⫽

n

→ p j · B j共␣ j ⫹ i␻ j兲2 · e(␣ ⫹i␻ )t 兺 j⫽1 j

j

ⱍ ⫽

(9) n

→ p j · 共B j共␣ j ⫹ i␻ j兲2 · e␣ ·t兲ei␻ ·t 兺 j⫽1 j

j

ⱍ ⫽

n

→ p j · C j · ei␻ ·t 兺 j⫽1 j

The oscillating terms in both Eqs. 8 and 9 are in the form ei␻ j·t → → ¨ oscillate at the same resonant demonstrating that both ⭸X˙ and ⭸X

(11)

'

⭸ JT ⭸␪

F (12)

C ␪ ⫽ J TC HJ

(8)

ⱍ ⫽

'

The position of the hand along the trajectory is needed to calculate H RH= (t). Its Cartesian trajectory is obtained from the Optotrak data (see Protocol). The transformation between the end-point stiffness, damping, and inertia and their respective joint equivalents is accomplished by means of the configuration-dependent Jacobian matrix J(␪):

M ␪ ⫽ J T M HJ

j

⫽ ␣j

H

K H ⫽ R H'K H' R H'

K ␪ ⫽ J TK HJ ⫹ n

(10)

A practical consequence of this relationship is that ␣j(t) can be obtained with simple accelerometers, and there is no need to employ sensitive motion tracking systems to precisely measure arm vibration. The eigenvectors of the vibrational modes are invariant with respect to the double derivation; therefore, a procedure to estimate the eigenvector matrix P, previously described in Piovesan et al. (2012c), can be applied using acceleration in place of position time signals. The acceleration signal used in the spectrogram is measured in the mobile reference frame centered on the accelerometer H=, which rotates and translates solidly with the wrist. The estimated stiffness expressed in the moving accelerometer reference frame is then projected onto the inertial H frame H by means of a time-varying rotation matrix RH= (t), namely:

ⱍ ⫽

⫽ ␣j

The estimation of the Jacobian matrix is subject specific and relying on individual anthropometric measurements and the kinematic configuration of the limb (Piovesan et al. 2011c). Because the modal ⭸JT analysis follows a free response, the term F ⫽ 0. ⭸␪ Resolution issue. With the proposed method, it is possible to map the behavior of the system in both the time and frequency domains. However, the resolutions of the spectrogram in time and frequency are constrained by the “bandwidth theorem,” which states that the product of the frequency and time resolutions is a constant. In addition, averaging across windows in the construction of the spectrogram results in a “smear” of energy density across time and frequency that limits the ability to identify precisely local amplitude maxima. We partially overcome this limitation by postprocessing the simple STFT with a reassignment technique (Fulop and Fitz 2006a, 2006b; Nelson 2001). Because the STFT is a complex function, phase information is available and is stationary around the local maximum in the amplitude spectrum. By computing the derivatives of the STFT’s phase with respect to time and frequency, we can identify time and frequency shifts of magnitude smaller than the original STFT resolution. This feature allows us to “reassign” the location of the maximum energy density in the time-frequency domain (Nelson 2001) corresponding to the resonant frequency ␻j(t). Our reassignment technique is less sensitive to noise than the commonly used Wigner-Ville transform (Boashash and O’Shea 1993), especially when it is applied to the energy distribution information of a multi-DOF system. The reassigned STFT has a constant time-frequency resolution and uniformly averaged-out noise. The Wigner-Ville transform has inherently high



and Aji is the instantaneous maximum of the spectrogram of the jth mode calculated for the impulse response of the ith DOF. The unknown variables of Eq. 4, ␩2j (t) and ⌫j(t), can be identified with Eqs. 5–7, by extracting the resonant frequency ␻j and the corresponding magnitude Aj from the spectrogram of the impulse response of each DOF i independently (Piovesan et al. 2012c). Our approach allows a continuous estimation of the resonant frequencies and damping of the system for as long as the perturbing effects of the force impulse are identifiable in the spectrogram. We tested the practical applicability of this theoretical observation by applying perturbations at three consecutive positions along the trajectory of the arm, comparing the stiffness and damping characteristics estimated from one perturbation to those stemming from the effects of perturbations applied at previous points. We show in RESULTS that stiffness and damping estimations could be made before perturbations had elicited any voluntary muscle activity to correct a movement trajectory. The resonant frequencies ␻j and corresponding magnitudes Aj can be derived from the spectrogram of either the position or the acceleration of the perturbed arm trajectory. The general solution of Eq. 1 → is the superposition of each damped mode of vibration sj (Piovesan et al. 2012c), which can be written in the following form: →

frequency ␻j. Consequently, from Eq. 7, the term ␣j(t) can be calculated as the ratio of the first derivative with respect to time of either the instantaneous amplitude of the acceleration or the instantaneous amplitude of the position, namely: A˙ j

A˙ ji(t)

2487

2488


resolutions in both time and frequency across different regions of the time-frequency plane. However, it is very sensitive to noise because of the nonlocal, nonlinear process by which higher accuracy is achieved and is more likely to generate artifacts than the STFT reassignment procedure. EMG Processing To identify neural correlates of stiffness modulation we analyzed the SEMG signals of four muscles, one flexor and one extensor for each joint involved in the movement. The level of cocontraction was assessed by measuring the electrical activity during nonperturbed trials. We wanted to separate the effect of the intrinsic properties of the muscle-tendon system on stiffness from that of neural control, including reflexive signals. Applying perturbations of three different magnitudes along three different directions allowed us to test whether reflexes were elicited and whether the perturbations induced any potentially voluntary activity to correct a movement trajectory. Before using the SEMG signals for analyzing reflexive responses, we verified post hoc that the background SEMG activity was statistically the same before the onset of each perturbation, as well as before each movement onset. We then analyzed EMG activity within 0.035-s time windows around each perturbation (t ⫽ 0 at the onset of the perturbation). The background activity (BG) was defined for each muscle as the average EMG signal within the interval ⫺0.035 to 0 s. Spinal activity reflects a fast feedback loop and appears as short-latency (SL) reflexes that are measurable within the 0.015– 0.050 s interval (Johnson et al. 1993; Tatton and Lee 1975). EMG activity 0.050 s after the perturbation onset is likely to be conditioned by transcortical feedback loops that give rise to medium-latency (ML) and long-latency (LL) reflexes (Mutha et al. 2008; Pruszynski et al. 2011). On the basis of visual inspection of the EMG data, we identified ML reflexes within the interval 0.050 – 0.085 s and LL reflexes within 0.085– 0.120 s. Variations in the EMG activity between 0.120 s and 0.165 s cannot be easily attributed to a defined feedback loop, as both voluntary and reflexive components could coexist within this interval, and the onset of voluntary activity cannot be uniquely isolated. We therefore assumed the voluntary component of muscular activity to be fully developed in the interval Vol ⫽ 0.165– 0.200 s, following the perturbation onset.

Stiffness Analysis We obtain a stiffness profile within a wide temporal window following each perturbation, allowing for comparison of results within overlapping time windows. Figure 2 shows how the time courses of the stiffness estimated from different perturbations follow similar trends. Thus we can compare stiffness estimated at 1/2d and 3/4d produced by perturbations applied at different previous instants along the trajectory. Statistical analysis shows that the measures obtained with different perturbations are comparable as reported in Table 1, in which subject is a random factor. The left side of the table compares the estimates calculated at 0.01 s after the perturbation at 1/2d and at 0.105 s after the perturbation at 1/4d, which are equivalent

RESULTS

Time-Frequency Analysis The magnitude of hand acceleration measured by the accelerometer ensemble during typical perturbed trials is shown in

Fig. 3. Time profile of force perturbation. Perturbations were measured at the table (with the force plate) and at the hand (with the load cell) in order to estimate the transmissibility of vibration along the robot.



Fig. 2. Examples of measured accelerations following a perturbation. Left: acceleration magnitude at the wrist of 1 subject perturbed with a type A perturbation (3 N at 165°) applied at 1/4d, 1/2d, and 3/4d. Center: stiffness estimations from each perturbation on left. Solid, dashed, and dotted lines represent elbow stiffness, shoulder stiffness, and interjoint stiffness, respectively. Right: reassigned spectrograms of the acceleration signal for each perturbation. Color code: red, 1/4d; green, 1/2d; light blue, 3/4d.

Fig. 2. Acceleration data for the complete movement are presented for trials in which a type A perturbation (3 N at 165°) is applied at three points along the trajectory. Oscillations following the application of each perturbation are evident. The reassigned spectrograms of the filtered signals in Fig. 2 show the time-varying loci of the two resonant frequencies elicited by perturbations applied at the three different positions along the trajectory. At any point along the trajectory, the lower resonant frequency is slightly higher than the frequency proper of the voluntary movement [which is in agreement with results obtained by Bennett and colleagues (1992)]. Typical force profiles measured with both the load cell and the force plate are shown in Fig. 3. The easily identifiable oscillation visible in the force plate signal (⬃30 Hz) represents the lowest resonance frequency of the table. The table’s vibration is transmitted to the robot’s base and, although very attenuated, through the robotic arm to the load cell, and it can be seen as a small residual oscillation in the force component along the y-axis. Our technique allows us to filter out the vibration generated by robot and table by partitioning the frequency band around the arm’s resonance frequencies, because the experimental apparatus vibrates at frequencies outside the band of interest for stiffness measures. Joint angular displacement and velocity profiles for unperturbed movements are presented in Fig. 4, showing movement duration to be ⬃0.5 s. Acceleration signals were high-pass filtered with a cutoff of 1.8 Hz to attenuate the carrier frequency proper of the movement without affecting higher frequencies in the measured spectrum.


relative points along the arm trajectory. The estimated elbow stiffness is statistically different across perturbation epochs and perturbation directions, while the other stiffness parameters are comparable across trials. This pattern relates to the fact that the relative direction of the perturbation with respect to the arm Table 1. Statistical effects of factors influencing each stiffness term

segment direction varied at each perturbation epoch. As a consequence, the perturbations were more or less aligned with the direction of the local vibrational eigenvector, and segmental reflexes of different intensity were elicited by different perturbations (discussed below in detail). The right side of Table 1 illustrates the statistical compatibility of joint stiffness estimated at the relative points along the trajectory located 0.01 s after the perturbation point at 3/4d and across perturbations applied at all three positions. In this analysis, the interjoint stiffness is statistically different across epochs but not across perturbation directions. This means that the orientation of the stiffness ellipses, but not their magnitude, is influenced by the point of application of the perturbation. This analysis verifies the feasibility of using a single perturbation to estimate the stiffness within a wide window of the trajectory. The value of each stiffness coefficient follows a very similar time profile, even though their estimation results from perturbations applied at different instants and along different directions relative to the arm. The length of the useful time window available for the estimation is limited only by the dissipation of elastic energy due to damping. To compare our results to related studies (Burdet et al. 2000; Darainy et al. 2007; Gomi and Kawato 1997), we computed a set of stiffness, damping, and inertial estimates at 0.135 s after the onset of each perturbation. Since the force signal follows the SEMG signal with a delay of ⬃0.050 s (Carter et al. 1993; Frolov et al. 2000), each set of parameters so computed corresponded to the instant at which force generation was affected by the approximate peak of a ML reflex, which is the prevalent reflex seen in this experiment. Our method requires an independent estimation of the limb’s inertia. Methods in the literature can give different values, and the intrasubject variability of the inertial parameter estimation across methods influences stiffness and damping calculations. We used the method proposed by Zatsiorsky (2002) as our benchmark because it provides estimations that are the closest to the average across several methods, as we found in a previous study (Piovesan et al. 2011c). For each subject, inertial end-point ellipses at the position of the estimations are presented in Fig. 5. The dispersion across trials of the end-point inertial matrix, shown in its elliptical representation, is due to the variability of the kinematic configuration at each point of estimation. Figure 5 also depicts the end-point stiffness and damping ellipses in a Cartesian reference frame for each

at 3/4d (1/4d at 0.200 s—1/2d at 0.105 s—3/4d at 0.01 s)

at 1/2d (1/4d at 0.105 s—1/2d at 0.01 s) Source

d.f.

KSS

KSE

KEE

KSS

KSE

KEE

Subject Inertial method Perturbation epoch Set Perturbation direction

12 8 1 11 2

0.008 ⬍0.001 0.13 0.58 0.18

⬍0.001 0.64 0.98 0.55 0.66

⬍0.001 ⬍0.001 0.009 0.39 0.01

0.004 ⬍0.001 0.07 0.70 0.36

⬍0.001 ⬍0.001 0.005 0.297 0.834

⬍0.001 ⬍0.001 0.06 0.53 0.56

P values are shown. Left side of the table compares estimates of each component of the joint stiffness matrix K. Subscript S refers to the shoulder, E to the elbow. Components are calculated at 0.01 s after the perturbation at 1/2d and at 0.105 s after the perturbation at 1/4d, which are equivalent relative points along the arm trajectory of length d. Right side illustrates statistical analysis of the joint stiffness components estimated at 0.01 s after the perturbation point at 3/4d, 0.105 s after the perturbation point at 1/2d, and 0.200 s after the perturbation point at 1/4d. The position in space coincides for the 3 estimations. Fixed factors are the inertial methods used to estimate the stiffness, the perturbation epochs, the sets of movement repetitions, and the perturbation direction. Subject is a random factor. d.f., Degrees of freedom. P values indicating a statistical influence of the perturbation epoch and perturbation direction are highlighted in bold. J Neurophysiol • doi:10.1152/jn.01013.2012 • www.jn.org


Fig. 4. Average EMG and corresponding movement kinematics. Top: average surface EMG (SEMG) across subjects normalized to baseline muscle activation. au, Arbitrary unit. Middle: average angular velocity at the joints. Bottom: average joint angles. Dashed line with light gray [standard deviation (STD)] dispersion represents shoulder joint kinematic; solid line with dark gray STD dispersion represents elbow joint kinematic. Vertical dashed lines correspond to start and finish of the movement, defined as 10% and 90%, respectively, of the total angular excursion. Vertical dotted lines indicate when perturbations were delivered.

2489

2490



Fig. 5. Ellipses of impedance components estimated 0.135 s after perturbation onset, represented in the Cartesian space, for subjects S1–S13. Dashed lines represent 95% confidence interval. Color code: red, 1/4d; green, 1/2d; light blue, 3/4d.

subject, calculated with the inertial properties represented in the same figure. The variability of the stiffness estimates originates mostly from the uncertainty in the calculation of the resonant frequencies and, to a lesser extent, as a consequence

of imprecision in the estimation of the other parameters (Piovesan et al. 2012c). From Fig. 5 it can be seen that subjects with higher inertia (e.g., S1, S10, and S13) have higher stiffness and subjects with low inertia (e.g., S2, S7, S11) have lower



2491

stiffness. This pattern suggests that stiffness scales with the inertia of the subject’s arm. Since the squared natural frequency of a system is proportional to the ratio between stiffness and inertia and the movement duration was similar across subjects, this result supports the hypothesis that one goal of stiffness modulation is to maintain the resonant frequency of the limb close to the frequency of the movement (Bennett et al. 1992). The orientations of the stiffness ellipses are consistent with those obtained by Gomi and Kawato (1997); however, compared with their results, arm stiffness estimates in the direction of movement are significantly higher at the first estimation point and lower at the last. These differences may be the consequence of paradigm differences between methods as explained in DISCUSSION. Figure 6 presents the means and standard deviations of the estimates of the stiffness matrix coefficients, calculated across trials, inertial estimation methods, and perturbation types, at 0.135 s after the perturbation onset, expressed both in a Cartesian inertial reference frame and in joint space. Elbow stiffness KEE is higher than shoulder stiffness KSS for all subjects. The interjoint stiffness term KSE is very small compared with the other parameters. Figure 7 presents an analysis of the sensitivity of our technique to variations across inertial estimation methods. The mean stiffness is averaged across subjects and trials. Table 2 shows that different inertial estimation methods produce stiffness estimates that differ significantly. The influence of other parameters of the paradigm not explicitly accounted for is also shown as the result of a multivariate ANOVA, where subject is a random factor. The estimated value of the stiffness does not statistically depend on the perturbation direction. In principle, this could indicate that

the perturbation direction’s effect is small compared with the large total variance of a generally noisy measure. However, the total variance of our stiffness estimations (Fig. 5) is quite small, suggesting that the perturbation direction is indeed not critical. The variation of the stiffness components along the trajectory is statistically relevant: Figs. 5–7 indicate that stiffness decreases between the first (1/4d) and last (3/4d) points of measurement. As shown in Table 2, stiffness estimates did not significantly differ across sets of movements, suggesting that subjects performed the task without fatiguing. Symmetry The assumptions of symmetry and oscillatory behavior of the system are discussed in a separate theoretical and procedural analysis (Piovesan et al. 2012c). For vibrational symmetric systems, eigenvectors are orthogonal and the corresponding components of the vibration, measured at each DOF, are in either phase or counterphase. It can be demonstrated that the eigenvectors of asymmetric system matrices are not orthogonal, which in the time-frequency domain translates to a phase → shift between the components of vibration s j. For each DOF, we measured the time delay between the peaks of vibrational → components s j associated to the corresponding eigenvectors. Our analysis of variance found no significant delay between modes (P ⬍ 0.001), indicating that the two modes are indeed in phase and counterphase respectively, therefore supporting the assumption of a symmetric system. EMG analysis. Figure 4 shows the SEMG signals of all monitored muscles along the whole trajectory. The signals were normalized using the background activity at rest before



Fig. 6. Stiffness components represented in joint space (left) and Cartesian space (right), averaged across inertial models for individual subjects. Stiffness is calculated 0.135 s after perturbation onset. Since the latency between the SEMG signals and the generated force is ⬃0.050 s, the presented stiffness components reflect the SEMG at 0.085 s after the onset, which corresponds to the medium-latency (ML) interval. KEE, elbow stiffness; KSS, shoulder stiffness; KSE, interjoint stiffness. Color code: red, 1/4d; green, 1/2d; light blue, 3/4d. For each subject the total number of trials is 324 (3 perturbations ⫻ 12 movement sets ⫻ 9 inertial models).

2492


each movement. We ensured that each trial started with the minimum activation necessary to self-support the limb by monitoring the SEMG activity of each muscle in real time. Figure 4 shows the average across subjects of both SEMGs and kinematics. Starting and target contact times are marked with vertical dashed lines and correspond to 10% and 90%, respectively, of the total angular excursion. Instants at which a perturbation was delivered are marked with vertical dotted lines and correspond to 1/4d, 1/2d, and 3/4d in Fig. 1. Peak velocity occurs at 1/2d along the trajectory, where the shoulder and elbow angles are 45° and 90°, respectively, which is approximately the configuration in which the muscles have their maximum moment arm (Murray et al. 2002) and the optimal length to maximize joint torques. Muscle activation increases after velocity peaks even though the measured stiffness magnitude decreases in the same trajectory segment (see Fig. 4). At any activation level, muscle force is dependent on muscle length. Muscles at the end of the movement are either extended (biceps, deltoids) or contracted (pectoralis, triceps) away from their optimal length for force

generation. They thus exert a smaller force for their activation level than they would at their optimal level. This pattern suggests that muscle activation can be a unique predictor for estimating stiffness magnitude in isometric (postural) conditions, where muscle length is known, but not during movement, when the length of the muscle varies. Stiffness is actually higher before the velocity peak, even though the overall muscle activation is lower and the arm is mostly flexed (muscles away from their optimal length). This result can be understood by analyzing the pattern of reflexes. Figure 8 shows the reflexive components of the SEMG before and after the perturbations at each position and for each perturbation direction. A statistical analysis of the difference in SEMG signals between perturbed and nonperturbed trials shows that reflexes during perturbed trials are more pronounced before the velocity peak, as reported in Table 3. For all position and perturbation directions the background SEMG (BG) is not statistically different between perturbed and nonperturbed movements. For position 1/4d, most of the reflexive activity in perturbed trials is statistically significant for all muscles. LL reflexes are elicited only

Table 2. Statistical effects of factors influencing each stiffness term Source

d.f.

KXX

KXY

KYY

KSS

KSE

KEE

Subject Inertial method Perturbation epoch Set Perturbation direction

12 8 2 11 2

⬍0.001 ⬍0.001 ⬍0.001 0.824 0.354

0.001 ⬍0.001 ⬍0.001 0.492 0.530

0.005 ⬍0.001 ⬍0.001 0.313 0.371

0.003 ⬍0.001 0.027 0.492 0.300

0.003 ⬍0.001 0.068 0.375 0.972

⬍0.001 ⬍0.001 ⬍0.001 0.641 0.633

P values are shown. Table illustrates effect of inertial estimation methods on assessment of both Cartesian stiffness (KXX, KXY, KYY) and joint stiffness (KSS, KSE, KEE) components across subjects and trials. Subject is a random factor. The estimated value of the stiffness does not statistically depend on the perturbation direction. In principle, this could indicate that the perturbation direction’s effect is small compared to the large total variance of a generally noisy measure. However, the total variance of our stiffness estimations (Fig. 5) is quite small, suggesting that the perturbation direction is indeed not critical. J Neurophysiol • doi:10.1152/jn.01013.2012 • www.jn.org


Fig. 7. Stiffness representations averaged across subjects as a function of inertial models: Hanavan (1964) (HV), Dempster (1955) (DE), Chandler et al. (1975) (CH), Clauser et al. (1969) (CL), McConville et al. (1980) (MC), Zatsiorsky (2002) (Z2), Piovesan et al. 2011c (PI), Zatsiorsky and Seluyanov (1983) (Z1), de Leva (1996) (DL). Stiffness was estimated 0.135 s after perturbation onset. Color code: red, 1/4d; green, 1/2d; light blue, 3/4d.


2493

when the perturbation is delivered prior to the velocity peak, indicating that the higher value of stiffness measured in the first part of the movement is correlated with stronger reflexive activity. Our minimally invasive perturbations did not lead to the subjects correcting a trajectory voluntarily after the impulse. This can be inferred since voluntary level (VOL) in perturbed trials is not statistically different from the baseline. This finding is consistent with results presented by Popescu and Rymer (2000). DISCUSSION

We estimated joint stiffness during reaching movements with a new technique using single impulsive perturbations. The effects of the perturbations on the electromyographic activity of four arm muscles were analyzed. We found that the perturbations did not cause any measurable voluntary response but did elicit reflexive responses, especially when delivered before the movement velocity peak was reached. The stiffness estimates were repeatable across multiple measures at the same point along the movement trajectory, and we could estimate a time-varying stiffness profile by applying only one perturbation. To confirm the validity of our results, we measured the stiffness time profiles elicited by perturbations of different magnitudes and directions applied at three points along movement trajectories, compared the resulting estimates at these three points, and found them statistically comparable.

Muscle mechanics are well represented by third-order models (Piovesan et al. 2012b), and modeling arm mechanics as a second-order system is an approximation. A sufficient condition for a third-order muscle model to have a complex pole that originates an oscillatory behavior is that the stiffness of the tendon is less than eight times the stiffness of the muscle, independently of the damping (Piovesan et al. 2012b). When approximating a third-order muscle system with a second-order model, an oscillatory behavior is expected because the complex pole pair must be dominant (if the real pole were dominant, the third-order system would be approximated by a first-order model). Direct measure of the residual oscillation at the subject’s styloid process (Fig. 2) confirms that the system is underdamped. This is in line with the results obtained with stochastic nonparametric identification (Perreault et al. 2004; Schouten et al. 2001; Westwick and Perreault 2011). The estimation of stiffness at different points along the trajectory is of particular interest when studying motor control strategies within a single movement. We observed marked variations of stiffness magnitude along the movement trajectory, in accord with previously published results (Gomi and Kawato 1997). Our results are in agreement, although in Gomi and Kawato’s study movements were executed in the opposite direction (distal to proximal) compared with our protocol. In both studies, joint stiffness was higher in the first half of the movement prior to the movement velocity peak and then decreased as the movement proceeded. We also measured



Fig. 8. Reflexive activity compared across recorded muscles as a function of perturbation direction and position along the trajectory. Each panel represents the average across subjects. Black traces represent muscle activity in unperturbed trials, colored traces represent the average of trials with perturbations. For each subject and condition, unperturbed traces are obtained with the same number of trials as perturbed traces. The matching number of unperturbed trials for each subject was selected with a randomization algorithm.

2494


Table 3. P values for reflexes with respect to baseline BG

ML

LL

VOL

1/4d

1/2d

3/4d

1/4d

1/2d

3/4d

1/4d

1/2d

3/4d

1/4d

1/2d

3/4d

1/4d

1/2d

3/4d

0.51 0.29 0.92

0.53 0.13 0.07

0.11 0.96 0.09

0.08 0.45 1.00

0.09 0.38 0.01

0.59 0.54 0.19

0.01 0.20 0.95

0.19 0.20 0.07

0.98 0.57 0.40

0.01 0.96 0.13

0.67 0.50 0.05

0.31 0.12 0.85

0.06 0.59 0.35

0.95 0.26 0.13

0.78 0.71 0.11

0.56 0.18 0.07

0.97 0.89 0.51

0.48 0.23 0.23

0.17 0.79 0.02

0.29 0.25 0.40

0.55 0.18 0.12

0.04 0.05 0.02

0.18 0.04 0.02

0.29 0.43 0.07

0.04 0.66 0.00

0.19 0.38 1.00

0.38 0.10 0.07

0.05 0.80 0.05

0.39 0.45 0.52

0.68 0.06 0.69

0.11 0.74 0.98

0.86 0.50 0.35

0.49 0.58 0.91

0.00 0.77 0.36

0.61 0.53 0.22

0.03 0.05 0.93

0.71 0.16 0.43

0.05 0.94 0.30

1.00 0.05 0.06

0.80 0.21 0.48

0.51 0.23 0.52

0.56 0.97 0.68

0.66 0.76 0.61

0.54 0.75 0.77

0.92 0.71 0.16

0.24 0.20 0.23

0.35 0.71 0.07

0.41 0.18 0.61

0.06 0.04 0.13

0.99 0.69 0.08

0.71 0.03 0.21

0.03 0.01 0.01

0.69 0.15 0.01

0.97 0.75 0.04

0.05 0.12 0.03

0.76 0.29 0.28

0.28 0.98 0.14

0.43 0.38 0.32

0.84 0.25 0.10

0.67 0.15 0.74

Statistical analysis of the surface electromyographic (SEMG) signal difference between perturbed and nonperturbed trials. A P ⬍ 0.05 (in bold) indicates that the reflexive activity is statistically higher than the baseline level. Perturbations were applied along the trajectory at 3 distinct points located at 1/4, 1/2, and 3/4 of total distance d, between the starting and ending points measured along the y-axis. Force perturbations had 3 possible directions and 3 associated magnitudes (A ⫽ 3 N at 165°, B ⫽ 4 N at 45°, and C ⫽ 5 N at 285°, with 0° aligned to the x-axis). Defining the time at which each perturbation is applied to be 0 s, the average of the SEMG amplitude was evaluated in the following time intervals: background (BG; ⫺0.035 to 0 s), short-latency reflexes (SL; 0.015 to 0.050 s), medium-latency reflexes (ML; 0.050 to 0.085 s), long-latency reflexes (LL; 0.120 to 0.165 s), voluntary activity (VOL; 0.165 to 0.200 s). BLH, long head of biceps brachii; PCH, clavicular head of pectoralis; TLH, triceps brachii lateral head; DPH, deltoid posterior head.

higher muscle reflexive activity in the first half of the movement, showing a correlation between high stiffness and reflexive activity. Interestingly, Perreault and colleagues (2001, 2004), using a pseudorandom perturbation technique, also found higher stiffness than the values estimated with regression techniques (Darainy et al. 2004; Mussa-Ivaldi et al. 1985). This contrast could be due to segmental reflexes being elicited by the high-frequency components of pseudorandom perturbations, which consequently would increase the stiffness. Indeed, the ramp and hold perturbations used for some regressive techniques have lower power within high-frequency bands compared with pseudorandom perturbations (Schouten et al. 2008; van der Helm et al. 2002). We observed a marked reduction in stiffness about and after the movement velocity peak coinciding with the decrease in LL reflex activity. In vivo experiments have demonstrated that reflexal muscles react to a perturbation with a springlike behavior while stretched areflexal muscles behave as dampers because of the muscle yield property (Lin and Rymer 1998, 2001). These observations are consistent with our experimental results that show a marked increase in damping associated with the decrease in reflexive activity. Several factors are responsible for the attenuation of segmental reflexes. Reflexes are lower during voluntary movements than in postural tasks requiring self-support of the limb against gravity (Bawa and Sinkjaer 1999). This could be linked to the increasingly stabilizing effect of inertia as the velocity increases (Gallego et al. 2010). Visual feedback could also indirectly influence the modulation of reflexes and therefore joint stiffness. Franklin and colleagues (2007) reported a decrease in stiffness after the exclusion of visual feedback. The exclusion of visual feedback could induce subjects to decrease muscle activity, and consequently joint stiffness, to minimize the contact force at the fingertip as it impacts the table. This also suggests that a force control strategy could be employed in the absence of visual

feedback and would be consistent with results showing a decrease of reflexive activity when force control strategies are implemented (Forbes et al. 2011). Recent studies also suggest that a decrease in stiffness is used when interacting with an external force during rehabilitation robotic training when visual feedback is precluded (Piovesan et al. 2011a, 2013). The stiffness ellipses estimated with our technique were oriented with their major axes approximately aligned with the direction of the movement. This result is comparable to those obtained with regressive techniques employing force perturbations (Gomi and Kawato 1997). In three experiments in which stiffness was estimated through displacement perturbations, the major axes of stiffness ellipses were found to be oriented orthogonally to the movement (Darainy et al. 2007; Wong et al. 2009a, 2000b). Because stiffness was estimated in all these experiments during unfettered reaching movements, in approximately the same arm configuration, the orientation discrepancy might be due to differences in impedance at the points of contact for application of the force and displacement perturbations (de Vlugt et al. 2002). To minimize insertion errors, when force perturbations are used the impedance at the point of contact must be minimal (theoretically null), while for displacement perturbations it needs to be as high as possible in the direction of the perturbation, which requires the perturbation to be delivered by a highly rigid robotic device. Experimental evidence (Chib et al. 2006) suggests that the same force applied to the hand of a reaching subject can be interpreted as a force field to be crossed or as a surface to be circumvented depending on whether the impedance at the point of contact is low or high. Usually arm stiffness is measured after training, when the decision of either going through or going around the perturbing field is well established. If the perturbation is perceived as a force field to be overcome, the subject presumably increases stiffness in the direction of movement. By contrast, if the perturbation is perceived as a hard



PCH A B C DPH A B C BLH A B C TLH A B C

SL


ACKNOWLEDGMENTS The authors thank Dr. Ferdinando Mussa-Ivaldi and the members of the Robotics Lab at Northwestern University for their support. GRANTS The present research was supported by National Institute of Arthritis and Musculoskeletal and Skin Diseases Grant AR-048546. DISCLOSURES No conflicts of interest, financial or otherwise, are declared by the author(s). AUTHOR CONTRIBUTIONS Author contributions: D.P. conception and design of research; D.P. and A.P. performed experiments; D.P. and A.P. analyzed data; D.P., P.D., and J.R.L. interpreted results of experiments; D.P. and A.P. prepared figures; D.P. and A.P. drafted manuscript; D.P., A.P., P.D., and J.R.L. edited and revised manuscript; D.P., P.D., and J.R.L. approved final version of manuscript. REFERENCES Bawa P, Sinkjaer T. Reduced short and long latency reflexes during voluntary tracking movement of the human wrist joint. Acta Physiol Scand 167: 241–246, 1999. Bennett DJ, Hollerbach JM, Xu Y, Hunter IW. Time-varying stiffness of human elbow joint during cyclic voluntary movement. Exp Brain Res 88: 433– 442, 1992. Boashash B, O’Shea P. Use of the cross Wigner-Ville distribution for estimation of instantaneous frequency. IEEE Trans Signal Processing 41: 1439 –1445, 1993. Burdet E, Osu R, Franklin DW, Yoshioka T, Milner TE, Kawato M. A method for measuring endpoint stiffness during multi-joint arm movements. J Biomech 33: 1705–1709, 2000. Carter RR, Crago PE, Gorman PH. Nonlinear stretch reflex interaction during cocontraction. J Neurophysiol 69: 943–952, 1993.

Casadio M, Morasso P, Sanguineti V, Giannoni P. Minimally assistive robot training for proprioception enhancement. Exp Brain Res 194: 219 –231, 2009. Chandler RF, Clauser CE, McConville JT, Reynolds HM, Young JW. Investigation of Inertial Properties of the Human Body. Dayton, OH: Wright-Patterson Air Force Base, 1975. Chib VS, Patton JL, Lynch KM, Mussa-Ivaldi FA. Haptic identification of surfaces as fields of force. J Neurophysiol 95: 1068 –1077, 2006. Clarys JP, Marfell-Jones MJ. Anatomical segmentation in humans and the prediction of segmental masses from intra-segmental anthropometry. Hum Biol 58: 761–769, 1986a. Clarys JP, Marfell-Jones MJ. Anthropometric prediction of component tissue masses in the minor limb segments of the human body. Hum Biol 58: 771–782, 1986b. Clauser CE, McConville JT, Young JW. Weight, Volume, and Center of Mass of Segments of the Human Body. Dayton, OH: Wright-Patterson Air Force Base, 1969. Darainy M, Malfait N, Gribble PL, Towhidkhah F, Ostry DJ. Learning to control arm stiffness under static conditions. J Neurophysiol 92: 3344 –3350, 2004. Darainy M, Towhidkhah F, Ostry DJ. Control of hand impedance under static conditions and during reaching movement. J Neurophysiol 97: 2676 – 2685, 2007. de Leva P. Adjustments to Zatsiorsky-Seluyanov’s segment inertia parameters. J Biomech 29: 1223–1230, 1996. de Vlugt E, Schouten AC, van der Helm FC. Adaptation of reflexive feedback during arm posture to different environments. Biol Cybern 87: 10 –26, 2002. Dempster WT. Space Requirements of the Seated Operator; Geometrical, Kinematic, and Mechanical Aspects of the Body with Special Reference to the Limbs. Dayton, OH: Wright Air Development, 1955. Flash T. The control of hand equilibrium trajectories in multi-joint arm movements. Biol Cybern 57: 257–274, 1987. Forbes PA, Happee R, van der Helm FC, Schouten AC. EMG feedback tasks reduce reflexive stiffness during force and position perturbations. Exp Brain Res 213: 49 – 61, 2011. Franklin DW, So U, Burdet E, Kawato M. Visual feedback is not necessary for the learning of novel dynamics. PloS One 2: e1336, 2007. Frolov AA, Dufosse M, Rizek S, Kaladjian A. On the possibility of linear modelling the human arm neuromuscular apparatus. Biol Cybern 82: 499 – 515, 2000. Frolov AA, Prokopenko RA, Dufosse M, Ouezdou FB. Adjustment of the human arm viscoelastic properties to the direction of reaching. Biol Cybern 94: 97–109, 2006. Fulop SA, Fitz K. Algorithms for computing the time-corrected instantaneous frequency (reassigned) spectrogram, with applications. J Acoust Soc Am 119: 360 –371, 2006a. Fulop SA, Fitz K. A spectrogram for the twenty-first century. Acoustics Today 2: 26 –33, 2006b. Galleani L, Cohen L. Direct time-frequency characterization of linear systems governed by differential equations. IEEE Signal Process Lett 11: 721–724, 2004. Gallego JA, Rocon E, Roa JO, Moreno JC, Pons JL. Real-time estimation of pathological tremor parameters from gyroscope data. Sensors 10: 2129 – 2149, 2010. Gomi H, Kawato M. Human arm stiffness and equilibrium-point trajectory during multi-joint movement. Biol Cybern 76: 163–171, 1997. Hanavan EP. A Mathematical Model of the Human Body. Dayton, OH: Wright-Patterson Air Force Base, 1964. Hogan N. An organizing principle for a class of voluntary movements. J Neurosci 4: 2745–2754, 1984. Hondori HM, Tech AW. Smart mug to measure hand’s geometrical mechanical impedance. Conf Proc IEEE Eng Med Biol Soc 2011: 4053– 4056, 2011. Hondori HM, Khademi M, Lopes CV. Use of a portable device for measuring arm’s planar mechanical impedance during motion. In: IEEE EMBS Conference on Biomedical Engineering and Sciences (IECBES) 2012, p. 980 –983, 2012. Houk JC. Regulation of stiffness by skeletomotor reflexes. Annu Rev Physiol 41: 99 –114, 1979. Imran S, Jamil A, Ismail SS, Kashif FM. Techniques to obtain good resolution and concentrated time-frequency distributions: a review. EURASIP J Adv Signal Process 2009: 1– 43, 2009. Inman DJ. Vibration: With Control, Measurement, Stability. Englewood Cliffs, NJ: Prentice Hall, 1989.



surface, decreasing stiffness in the direction of the movement could represent a strategy to minimize contact force. These considerations suggest that mechanisms underlying stiffness modulation employ feedback information from tactile receptors, but this hypothesis requires further experimental evidence. Movement velocity might influence stiffness orientation as well, but further evidence is also needed to support this. The estimation of stiffness provided by our method is sensitive to the variability of inertial parameters. The set of anthropometric parameters we used was computed with the regression method of Zatsiorsky and Seluyanov (Zatsiorsky 2002), which we identified as the best available technique in our previous work (Piovesan et al. 2011c). Nevertheless, the results of the statistical analyses we presented are consistent across different inertial models, when the subject is a random factor. The method we have proposed to measure stiffness is particularly apt for studying movements in a noninertial environment. The deliverance of a small and brief force pulse does not interfere with contact impedance, which is null in unfettered movements. In addition, the capacity of the technique to estimate joint stiffness on a single-trial basis is suitable for the study of stiffness modulation during adaptation processes (Kolesnikov et al. 2011; Lackner and Dizio 1994; Shadmehr and Mussa-Ivaldi 1994) and in subjects with pathologies that reduce the repeatability of their movements (Casadio et al. 2009; Piovesan et al. 2011a, 2011b, 2012a).

2495

2496

ARM STIFFNESS MEASURE WITH TIME-FREQUENCY ANALYSIS Piovesan D, DiZio P, Lackner JR. A new time-frequency approach to estimate single joint upper limb impedance. Conf Proc IEEE Eng Med Biol Soc 2009: 1282–1285, 2009. Piovesan D, Morasso P, Giannoni P, Casadio M. Arm stiffness during assisted movement after stroke: the influence of visual feedback and training. IEEE Trans Neural Syst Rehabil Eng 21: 454 – 465, 2013. Piovesan D, Pierobon A, DiZio P, Lackner JR. Comparative analysis of methods for estimating arm segment parameters and joint torques from inverse dynamics. J Biomech Eng 133: 031003, 2011c. Piovesan D, Pierobon A, DiZio P, Lackner JR. Measuring multi-joint stiffness during single movements: numerical validation of a novel timefrequency approach. PLoS One 7: e33086, 2012c. Piovesan D, Casadio M, Mussa-Ivaldi FA, Morasso P. Comparing two computational mechanisms for explaining functional recovery in robottherapy of stroke survivors. Annu Int Conf IEEE Biomed Robot Biomechatronics 2012: 222, 2012a. Piovesan D, Pierobon A, Mussa-Ivaldi FA. Third-order muscle models: the role of oscillatory behavior in force control. Proc Int Mech Eng Cong Exposition ASME-IMECE, 2012b. Popescu FC, Rymer WZ. End points of planar reaching movements are disrupted by small force pulses: an evaluation of the hypothesis of equifinality. J Neurophysiol 84: 2670 –2679, 2000. Pruszynski JA, Kurtzer I, Scott SH. The long-latency reflex is composed of at least two functionally independent processes. J Neurophysiol 106: 449 – 459, 2011. Schouten AC, de Vlugt E, van Der Helm FC, Brouwn GG. Optimal posture control of a musculo-skeletal arm model. Biol Cybern 84: 143–152, 2001. Schouten AC, de Vlugt E, van Hilten JJ, van der Helm FC. Quantifying proprioceptive reflexes during position control of the human arm. IEEE Trans Biomed Eng 55: 311–321, 2008. Seki K, Perlmutter SI, Fetz EE. Sensory input to primate spinal cord is presynaptically inhibited during voluntary movement. Nat Neurosci 6: 1309 –1316, 2003. Shadmehr R, Mussa-Ivaldi F. Adaptive representation of dynamics during learning of a motor task. J Neurosci 14: 3208 –3224, 1994. Tatton WG, Lee RG. Evidence for abnormal long-loop reflexes in rigid Parkinsonian patients. Brain Res 100: 671– 676, 1975. Tsuji T, Morasso PG, Goto K, Ito K. Human hand impedance characteristics during maintained posture. Biol Cybern 72: 475– 485, 1995. van der Helm FC, Schouten AC, de Vlugt E, Brouwn GG. Identification of intrinsic and reflexive components of human arm dynamics during postural control. J Neurosci Methods 119: 1–14, 2002. Westwick DT, Perreault EJ. Closed-loop identification: application to the estimation of limb impedance in a compliant environment. IEEE Trans Biomed Eng 58: 521–530, 2011. Won J, Hogan N. Stability properties of human reaching movements. Exp Brain Res 107: 125–136, 1995. Wong J, Wilson ET, Malfait N, Gribble PL. The influence of visual perturbations on the neural control of limb stiffness. J Neurophysiol 101: 246 –257, 2009a. Wong J, Wilson ET, Malfait N, Gribble PL. Limb stiffness is modulated with spatial accuracy requirements during movement in the absence of destabilizing forces. J Neurophysiol 101: 1542–1549, 2009b. Zatsiorsky V, Seluyanov V. The mass and inertia characteristics of the main segments of the human body. In: International Congress of Biomechanics: Biomechanics VIII-B, edited by Matsui H, Kobayashi K. Champaign, IL: Human Kinetics, 1983, p. 1152–1159. Zatsiorsky VM. Best predictive regression equations for estimating inertial properties of body segments in males, Appendix A2.8. In: Kinetics of Human Motion, edited by Zatsiorsky VM. Champaign, IL: Human Kinetics, 2002, p. 600 – 601.



Johnson MT, Kipnis AN, Lee MC, Ebner TJ. Independent control of reflex and volitional EMG modulation during sinusoidal pursuit tracking in humans. Exp Brain Res 96: 347–362, 1993. Kearney RE, Hunter IW. Dynamics of human ankle stiffness: variation with displacement amplitude. J Biomech 15: 753–756, 1982. Kolesnikov M, Piovesan D, Lynch K, Mussa-Ivaldi F. On force regulation strategies in predictable environments. Conf Proc IEEE Eng Med Biol Soc 2011: 4076 – 4081, 2011. Lackner JR, Dizio P. Rapid adaptation to Coriolis force perturbations of arm trajectory. J Neurophysiol 72: 299 –313, 1994. Lin DC, Rymer WZ. Damping actions of the neuromuscular system with inertial loads: human flexor pollicis longus muscle. J Neurophysiol 85: 1059 –1066, 2001. Lin DC, Rymer WZ. Damping in reflexively active and areflexive lengthening muscle evaluated with inertial loads. J Neurophysiol 80: 3369 –3372, 1998. Loram ID, Lakie M, Di Giulio I, Maganaris CN. The consequences of short-range stiffness and fluctuating muscle activity for proprioception of postural joint rotations: the relevance to human standing. J Neurophysiol 102: 460 – 474, 2009. Loram ID, Maganaris CN, Lakie M. The passive, human calf muscles in relation to standing: the non-linear decrease from short range to long range stiffness. J Physiol 584: 661– 675, 2007a. Loram ID, Maganaris CN, Lakie M. The passive, human calf muscles in relation to standing: the short range stiffness lies in the contractile component. J Physiol 584: 677– 692, 2007b. Mah CD. Spatial and temporal modulation of joint stiffness during multijoint movement. Exp Brain Res 492–506, 2001. Masia L, Sandini G, Morasso PG. A novel mechatronic system for measuring end-point stiffness: mechanical design and preliminary tests. IEEE Conf Rehabil Robot 2011: 5975515, 2011. McConville JT, Churchill TD, Kaleps I, Clauser CE, Cuzzi J. Anthropometric Relationships of Body and Body Segment Moments of Inertia. Dayton, OH: Wright-Patterson Air Force Base, 1980. Murray WM, Buchanan TS, Delp SL. Scaling of peak moment arms of elbow muscles with upper extremity bone dimensions. J Biomech 35: 19 –26, 2002. Mussa-Ivaldi FA, Hogan N, Bizzi E. Neural, mechanical, and geometric factors subserving arm posture in humans. J Neurosci 5: 2732–2743, 1985. Mutha PK, Boulinguez P, Sainburg RL. Visual modulation of proprioceptive reflexes during movement. Brain Res 1246: 54 – 69, 2008. Nelson DJ. Cross-spectral methods for processing speech. J Acoust Soc Am 110: 2575–2592, 2001. Nichols TR, Houk JC. Improvement in linearity and regulation of stiffness that results from actions of stretch reflex. J Neurophysiol 39: 119 –142, 1976. Osu R, Gomi H. Multijoint muscle regulation mechanisms examined by measured human arm stiffness and EMG signals. J Neurophysiol 81: 1458 –1468, 1999. Perreault EJ, Kirsch RF, Crago PE. Effects of voluntary force generation on the elastic components of endpoint stiffness. Exp Brain Res 141: 312–323, 2001. Perreault EJ, Kirsch RF, Crago PE. Voluntary control of static endpoint stiffness during force regulation tasks. J Neurophysiol 87: 2808 –2816, 2002. Perreault EJ, Kirsch RF, Crago PE. Multijoint dynamics and postural stability of the human arm. Exp Brain Res 157: 507–517, 2004. Piovesan D, Casadio M, Morasso P, Giannoni P. Influence of visual feedback in the regulation of arm stiffness following stroke. Conf Proc IEEE Eng Med Biol Soc 2011: 8239 – 8242, 2011a. Piovesan D, Casadio M, Mussa-Ivaldi FA, Morasso PG. Multijoint arm stiffness during movements following stroke: implications for robot therapy. IEEE Int Conf Rehabil Robot 2011: 5975372, 2011b.

hand-related action verbs with and without emotional connotation.

The effect of arm weight support on upper limb muscle synergies during reaching movements.

Three-dimensional motion analysis of arm-reaching movements in healthy and hemispinalized common marmosets.

Submovements during reaching movements after stroke.

Children's development of reaching: temporal and spatial aspects of aimed whole-arm movements.

Stereotypical reaching movements of the octopus involve both bend propagation and arm elongation.

Object size modulates fronto-parietal activity during reaching movements.

Computing Arm Movements with a Monkey Brainet.

Repetitive Arm Movements During Sleep: A Polysomnographic Assessment.

Virtual trajectories, joint stiffness, and changes in the limb natural frequency during single-joint oscillatory movements.

A preferred pattern of joint coordination during arm movements with redundant degrees of freedom.

Cerebral cortical mechanisms of reaching movements.

Bimanual coupling effects during arm immobilization and passive movements.

A kinematic analysis of reaching and grasping movements in a patient recovering from optic ataxia.

An automatic experimental apparatus to study arm reaching in New World monkeys.

Neurons related to reaching-grasping arm movements in the rostral part of area 6 (area 6a beta).

Experimental Treatment of Ebola Virus Disease with TKM-130803: A Single-Arm Phase 2 Clinical Trial.

Arm and trunk movement kinematics during seated reaching within and beyond arm's length in people with stroke: a validity study.

Parcellation of sensorimotor transformations for arm movements.

Reorganization of muscle synergies during multidirectional reaching in the horizontal plane with experimental muscle pain.

Low-Dimensional Synergistic Representation of Bilateral Reaching Movements.

Minimizing endpoint variability through reinforcement learning during reaching movements involving shoulder, elbow and wrist.

Online gain update for manual following response accompanied by gaze shift during arm reaching.

Processes controlling arm movements in monkeys.