## Abstract

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.

- joint stiffness measure
- time-frequency analysis
- reflex modulation

the resistance to motion 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; 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 movement 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).

## 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 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 self-support 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 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

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).

### BSP Estimations

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

#### 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. (2012c). Classical methods for estimating stiffness assume the system to be time invariant (stationary). The typical model for a time-invariant 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: (1)
where
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. and are the first and second derivative of the position vector with 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): (2)
where matrix *S*_{m} is such that (3)
*P* is the eigenvector matrix, *I*_{n} is the identity matrix, and η_{j}^{2}(*t*) and Γ_{j}(*t*) are the natural frequency and modal damping ratio of the *j*th mode (Piovesan et al. 2012c). *Equation 1* can be recast as the following homogeneous diagonal system: (4)

Our technique estimates the solution of *Eq. 4* based on the analysis of the system's response to an impulsive perturbation in the time-frequency domain. At each instant *t*, the spectrogram of the impulse response of each DOF *i* presents as many amplitude peaks *A*_{j} = *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 *j*th 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) (5) (6)
where (7)
and *A*_{ji} is the instantaneous maximum of the spectrogram of the *j*th mode calculated for the impulse response of the *i*th DOF. The unknown variables of *Eq. 4*, η_{j}^{2}(*t*) and Γ_{j}(*t*), can be identified with *Eqs. 5–7*, by extracting the resonant frequency ω_{j} and the corresponding magnitude *A*_{j} 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 *A*_{j} 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
(Piovesan et al. 2012c), which can be written in the following form: (8)
where
is the eigenvector associated with the *j*th vibrational mode. The system is second order; therefore the frequency content of
and are the same, as can be seen from calculating the second derivative of *Eq. 8* with respect to time: (9)The oscillating terms in both *Eqs. 8* and *9* are in the form demonstrating that both
and oscillate at the same resonant 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: (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 frame *H* by means of a time-varying rotation matrix *R*_{H′}^{H}(*t*), namely: (11)The position of the hand along the trajectory is needed to calculate *R*_{H′}^{H}(*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*(θ): (12)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 analysis follows a free response, the term .

#### 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 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.

## RESULTS

### Time-Frequency Analysis

The magnitude of hand acceleration measured by the accelerometer ensemble during typical perturbed trials is shown in 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.

### 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 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 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 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 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 *K*_{EE} is higher than shoulder stiffness *K*_{SS} for all subjects. The interjoint stiffness term *K*_{SE} 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
. For each DOF, we measured the time delay between the peaks of vibrational components 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 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 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 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 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).

## 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.

## ACKNOWLEDGMENTS

The authors thank Dr. Ferdinando Mussa-Ivaldi and the members of the Robotics Lab at Northwestern University for their support.

- Copyright © 2013 the American Physiological Society