## Abstract

The perturbation method has been used to measure stiffness of the human arm with a manipulator. Results are averages of stiffness during short perturbation intervals (<0.4 s) and also vary with muscle activation. We therefore propose a novel method for estimating static arm stiffness from muscle activation without the use of perturbation. We developed a mathematical muscle model based on anatomical and physiological data to estimate joint torque solely from EMG. This model expresses muscle tension using a quadratic function of the muscle activation and parameters representing muscle properties. The parameters are acquired from the relation between EMG and measured torque. Using this model, we were able to reconstruct joint torque from EMG signals with or without co-contraction. Joint stiffness is directly obtained by differentiation of this model analytically. We confirmed that the proposed method can be used to estimate joint torque, joint stiffness, and stiffness ellipses simultaneously for various postures with the same parameters and produces results consistent with the conventional perturbation method.

## INTRODUCTION

The human musculoskeletal system has a spring-like property called stiffness, essential to our interaction with our surroundings and controlled by motor command and neural feedback systems. Two degree of freedom manipulators have been used to measure stiffness during posture (Darainy et al. 2004; Dolan et al. 1993; Flash and Mussa-Ivaldi 1990; Gomi and Osu 1998; Mussa-Ivaldi et al. 1985; Osu and Gomi 1999; Perreault et al. 2001; Tsuji et al. 1995) and movement (Burdet et al. 2001; Franklin et al. 2003; Gomi and Kawato 1996, 1997). Mussa-Ivaldi et al. (1985) developed a perturbation method to measure hand stiffness by using a manipulator to displace the subject's hand during maintenance of a given posture. Dolan et al. (1993) and Tsuji et al. (1995) extended the perturbation method to include measurement of dynamic components such as viscosity and inertia, as well as stiffness. Gomi and Kawato (1996, 1997) developed the parallel link drive air-magnet floating manipulandum (PFM) to measure dynamic hand stiffness during two joint arm movements. In the same laboratory, the most complete study was performed by Gomi and Osu (1998), who found that stiffness was dependent on endpoint force properties during force regulation tasks. Although these studies established the concept and basic experimental approach to measuring stiffness, the perturbation method is tedious and relies on trial-to-trial repeatability. A multiple-input and multiple-output recursive least-squares system identification algorithm was implemented to provide real-time estimates of stiffness (Bennett et al. 1992; Perreault et al. 2002, 2004). The remarkable study was performed by Perreault et al. (2002), who provided their subjects with real-time visual feedback of endpoint stiffness every 500 ms. They used a robotic manipulator for the characterization of stiffness with a real-time parametric identification algorithm on measured force and position data from stochastic displacement perturbation. Although this perturbation method could be updated at each sample instant (10 ms), a large manipulator was still needed for perturbation. This creates difficulty in estimating arm stiffness on the sagittal plane because the weight of the manipulator itself exerts influence on the task. It is also difficult to apply to daily tasks such as lifting an object.

One of our primary goals was to estimate stiffness using only EMG measurement. Most of the previous studies have reported that muscle activation is closely related to stiffness (Darainy et al. 2004; Flash and Mussa-Ivaldi 1990; Franklin et al. 2003; Gribble et al. 2003; Kearney and Hunter 1990; Mussa-Ivaldi et al. 1985; Osu and Gomi 1999; Osu et al. 2002; Perreault et al. 2001; Tsuji et al. 1995;). Osu and Gomi (1999) attempted to reconstruct joint stiffness from measured EMG levels and determine conversion factors to match the EMG to the stiffness measured by the perturbation method. However, the conversion factors were different among the positions and the data groups. Based on the previous study, Osu et al. (2002) proposed an index of muscle co-contraction around the joint (IMCJ) that was defined as the summation of absolute values of muscle quasi-tension calculated from EMG signals. They estimated a linear relationship between the IMCJ and joint stiffness measured by the perturbation method, and so they converted the IMCJ into stiffness units (rIMCJ) to agree with the PFM measured stiffness. However, they did not consider cross-joint stiffness (Rcj), so a Jacobian matrix could not be used to transform the rIMCJ into endpoint stiffness. In addition, it was previously assumed (Osu and Gomi 1999; Osu et al. 2002) that moment arm was constant, which is an oversimplification.

This study proposes a mathematical myokinetic model (Mykin model) that is able to estimate joint torques from EMG levels. In this model, the human arm on the horizontal plane was modeled as a two-link manipulator with six monoarticular muscles and two biarticular muscles. Each muscle tension was expressed by a quadratic function of the muscle activation represented by rectified, filtered, and normalized EMG signals. For each of the experimental subjects, we calculated the parameters of the eight muscles from the muscle activations and torque relationships. Joint stiffness was directly obtained from muscle activations by a stiffness equation, which differentiates the torque equation. To confirm the effectiveness of this method, we studied the relation between stiffness estimated by the proposed Mykin model and stiffness measured using the PFM. During maintained postures, we used the proposed method to estimate joint torque and stiffness simultaneously in time series and also to draw hand stiffness ellipses directly from EMG signals. This is the first report of a method for estimating stiffness and stiffness ellipses from only EMG signals and without the use of hand perturbation.

## METHODS

Nine subjects, ranging in age from 22 to 41 yr (8 males, right-handed; 1 female, ambidextrous) participated in this study. Each subject first performed *experiment 1* to calibrate the parameters of the Mykin model and the parameters of the reference total co-contraction levels (TCL). The trials in *experiment 1* took 2 h, and 30-min rest periods were allowed during the Mykin parameter calibration. After this calibration, we used the Mykin model to estimate stiffness and joint torques simultaneously. *Experiment 2* was performed to evaluate the accuracy of the proposed method by comparing it. *Subjects A–E* accomplished the perturbation experiment while in a relaxed state, and *subjects F–I* executed it during the co-contraction state. Each participant performed *experiments 1* and *2* on the same day. *Subject D* was asked to repeat *experiment 2* 3 days after performing the first experiment (Table 1).

### Experimental apparatus

We used two servo systems to measure stiffness. The PFM was used to evaluate arm stiffness while the subject was in a relaxed state (*subjects A–E*). Joint angle of the subject was calculated from hand position, which was measured with position sensors (409,600 pulse/rev) of the PFM. The force exerted by the hand was measured using a force sensor (resolution, 0.059 N). Position and force were sampled at 500 Hz, and EMG was sampled at 2 kHz (Fig. 1*A*). A linear motor system was also used to measure hand stiffness during subject regulated co-contraction (*subjects F–I*). The arm posture of the subject was measured with an OPTOTRAK system (resolution: 0.02 mm), the marker position of which was sampled at 200 Hz. Hand displacement was also measured with the linear motor system (resolution: 0.0005 mm). The force vector between the hand and the handle was measured with a force sensor (resolution: 0.05 N for both the *x*- and *y*-axes) attached to the motor table (Fig. 1*B*).

Each subject sat on a chair in front of the manipulandum with his shoulders fixed to the back of the chair by means of a shoulder harness (Fig. 1, *A* and *B*). The right forearm was securely coupled to the handle using a rigid custom molded thermoplastic cuff and supported against gravity by a horizontal beam. The cuff immobilized the wrist joint, permitting the exertion of shoulder and elbow joint torque only in the horizontal plane. Nine hand positions were examined in *experiments 1* and *2*. These positions formed a 3 × 3 grid, with 10 cm in each direction between hand positions. These hand positions created shoulder postures with joint angles from 12 to 70° and elbow postures with joint angles from 48 to 110°. EMG, position, and force of the hand were recorded during all trials of experiments.

##### EXPERIMENT 1.

Each muscle of the Mykin model has parameters representing muscle properties. In *experiment 1*, the subject performed two kinds of isometric contraction tasks to calibrate the Mykin model parameters from the EMG–torque relationships. Generally, humans can produce two kinds of regulation in the isometric condition.

*Task 1* was a force-regulation task performed using as little co-contraction as possible. A monitor in front of the subject displayed the current force vector and a moving circle indicated the target force and force direction (Fig. 1*C*). The subjects were asked to exert 20 N of force in 16 different radial target directions in the plane of the table at each of the nine hand positions. Each radial target direction was separated by 22.5°. One trial took 24 s, during which the target circle moved 1.5 rev. Each subject performed 4 trials in each position, resulting in a total of 36 trials.

*Task 2* was a posture maintenance and force regulation task, with five different levels of co-contraction. A monitor displayed the current force vector (Fig. 1*D*). TCL was displayed in three bar graphs for shoulder, elbow, and total TCL during the experiment. TCL is defined as the summation of the absolute values of muscle tension (1) (2) After *task 1*, we instantly calculated β_{ji} for each position. Reference lines were marked on the current TCL bar graph. The reference lines were determined by taking the average TCL while regulating a 5-Nm joint torque to 1. The subjects were asked to generate self-paced co-contractions in accordance with TCL while not producing any external force (keeping the value of the force vector at 0). Each subject was asked to perform 2 trials in each position resulting in a total of 18 trials.

During the tasks, the subject's right hand was kept at a specified position by the manipulandum. After performing *experiment 1*, we calculated the parameters of the muscles from the muscle activations and torque relationships.

##### EXPERIMENT 2.

The accuracy of the proposed method was tested by comparing the estimated stiffness with stiffness measured using the perturbation method. In *experiment 2*, the subject was asked to position his arm in one of nine hand positions and maintain TCL using the TCL bar graphs on the monitor. To avoid voluntary responses to perturbation by the subject, the handle and the subject's arm were hidden by a cover, and the subject was asked not to push the handle in any direction (which could be determined from the monitor) or voluntarily react to the perturbations.

While in a relaxed state, an external disturbance was applied to the subject's hand with the PFM. The temporal changes of hand displacement and hand force were measured with the PFM during the perturbation period, and stiffness was calculated by the perturbation method. At the same time, EMG was recorded, and stiffness was estimated using the Mykin method. Stiffness ellipses from the perturbation method were computed for subjects on the basis of ∼40 servoed displacements with eight different directions according to the conventional analysis, during which the subjects were asked to maintain posture in each position. Stiffness ellipses for the proposed method were also computed from the average of 40 trials during a period 0.4 s before perturbation. Each subject was asked to perform 360 trials in each position. The PFM and setup are described in detail elsewhere (Gomi and Kawato 1996, 1997). Note that we used two different perturbation systems. Because these experiments required a considerable amount of time and caused subjects to tire, we were unable to perform both experiments on the same day.

For the co-contraction state, stiffness was measured by applying small force perturbations. The linear motor system pushed the subject's hand 10 mm in the *x*+-direction and pulled it back within 0.4 s (Fig. 1*E*). The displacement amplitude caused by the force perturbation was nearly the same in all trials. The subject was asked to position his hand in each position, generate co-contraction, and maintain it using the TCL bar graphs. External force disturbance was applied to the subject's hand only along the *x*+-axis by the linear motor. Each subject was asked to perform 25 trials in each position, resulting in a total of 225 trials.

### Mykin muscle model

Earlier arm models have been proposed using two kinematic degrees of freedom in the horizontal plane (Feldman et al. 1990; Flanagan et al. 1993; Gribble et al. 1998; Katayama and Kawato 1993; Osu and Gomi 1999; Prilutsky 2000). Based on these studies, we suggested the Mykin model, in which the human arm in the horizontal plane is modeled as a two-link manipulator with six monoarticular muscles and two biarticular muscles. In Fig. 2, the mechanical action of the bicep is shown as a rack-and-pinion gear and a spring connected in series (Ghez 2000). Here, we used the Kelvin-Voigt model, consisting of an elastic element for static isometric contraction (Özkaya and Nordin 1991). Muscle tension *T* is determined from muscle stiffness *k*(*u*) and the stretch length of a muscle [*l*_{r}(*u*) − *l*(θ)] as follows (3) Here, the parameters *k*_{0} and *k*_{1} of the muscle stiffness *k*(*u*) are the intrinsic elasticity and the elasticity, respectively. The muscle stiffness *k*(*u*) is represented by linear functions of muscle activation *u* (Katayama and Kawato 1993; Osu and Gomi 1999). For muscles, two different changes of muscle length are distinguished. First, *l*_{r}(*u*) denotes the deviation of equilibrium length. We assumed that *l*_{r}(*u*) can be contracted by only the muscle activation *u. L*_{bias} represents a bias term and Δ*L*_{eq} is the deviation of equilibrium length between the equilibrium length *L*_{eq}(0) and the current equilibrium length *L*_{eq}(*u*). *L**_{bias} denotes a bias term that is set so that [*l*_{r}(*u*) − *l*(θ)] is not negative. The parameter *l*_{0} is the intrinsic rest length when *u* is zero, and *l*_{0} is a constant.

Second, *l*(θ) denotes the current muscle length with the current joint angle θ. *L*(0) denotes the muscle length when the joint angle is 0; *L*(θ) denotes the current muscle length at the current joint angle. It can be simplified as arc (*A*^{T}θ). *Equation 4* expresses each muscle's tension *T*_{i} as a quadratic function of the muscle activation *u*_{i} (4) where *a*_{ji} > 0(*i*: 1, 2, …,8; *j* = *s, e*) is the length of the moment arm. Each moment arm **A**_{n} for different postures is denoted as shown in *Eq. 5* (5) Joint torque τ_{j} is determined using muscle tension *T*_{i} acting on joint *j* (6) Joint stiffness is defined by the following differential operator (7) Assuming that muscle activation *u*_{i} is proportional to the EMG level, we estimate stiffness directly using *Eq. 7*. From the estimated joint stiffness *R*_{ji}, the hand stiffness *K* can be obtained using the following equation (8) Here, *J* denotes the Jacobian matrix of kinematics transformation. From this hand stiffness matrix *K*, a stiffness ellipse can be drawn. To summarize the hand stiffness ellipse, we use the size *A*, shape *R* (ratio of the major and minor axis length), and orientation φ (its relative angle to the shoulder–hand direction) of the stiffness ellipse as represented in *Eqs. 9*–*11* (9) (10) (11) where , and φ_{h} is the angle between shoulder–hand direction and *x*-axis. To compare the proposed method to the perturbation method, we used the relative error δ*x*, which is defined as (12) where *S*_{Mykin} is a characteristic value estimated from the proposed method, and *S*_{PFM} is a characteristic value obtained from the perturbation method.

### Muscle activation u

Several studies have investigated the relationship between surface EMG signals and muscle force (Basmajian and De Luca 1985; Clancy and Hogan 1997; Clancy et al. 2006; Gottlieb and Agarwal 1971; Inman et al. 1952; Koike and Kawato 1995; Maton et al. 1987). Surface EMG signals are spatiotemporally convoluted action potentials of the muscle membranes and involve not only descending central motor commands, but also reflex motor commands generated from sensory feedback signals. The muscle activation therefore is presumed to contribute to an increase in joint torque and stiffness. The rectified, filtered, and normalized EMG signals were used as the muscle activation *u* of *Eqs. 3*, *4*, *6*, and *7*. We recorded EMG using pairs of Ag-AgCl surface electrodes in a bipolar configuration. EMG was sampled at 2 kHz with a 12-bit resolution. The signal was digitally rectified, integrated for 0.5 ms, sampled at 200 Hz, and filtered with a second-order, low-pass filter (cut-off frequency: 2.2 Hz) as follows (13) where *h*_{j} is the FIR filter, *j* is discrete time, and *EMG* represents the rectified and integrated EMG signal. The second-order frequency response of the filter *H*(*s*) is defined as (14) where ω_{n} and ξ denote natural frequency and damping coefficient, respectively. The impulse response of *H*(*s*) is (15) The coefficients *h*_{j} in *Eq. 13* were acquired by digitizing *h*(*t*) (Koike and Kawato 1995). The filtered EMG (*fEMG*_{i}) signals were normalized.

A simple normalization technique was used to operate the Mykin model because of intrasubject variability in maximum vocuntary Contraction (MVC), the need for preliminary training, and the difficulty in performing a maximum exertion. In our model, the EMG signals of each muscle were normalized with the different minimum value according to the contribution to torque. As shown in Fig. 3, we replotted the maximum values, the mean values, and the minimum values of the filtered EMG signals for each range of torque from the force regulation data of position 1. Because EMG values over the minimum for each range of torque were suggested to generate surplus torque stemming from co-contraction, we normalized and found the best fit line for the minimum values using a least-squares method. We set the minimum EMG level of each muscle while regulating an 80-Nm joint torque to 1. We assumed that the 80 Nm would be adequate for the maximum torque considering the stiffness parameter. Even if we used a larger value such as 100 Nm, it made no difference because the stiffness parameter *k*_{1} became smaller.

The normalized EMG levels were used as the muscle activation *u* for estimating the Mykin model parameters. During the experiment, surface EMG was recorded from eight muscles as shown in Fig. 4*A*. For flexion/extension of the shoulder joint, we measured EMG from the deltoid clavicular part (DELC), pectralis major clavicular head (PMJC), and deltoid-scapular part (DELS). For shoulder–elbow double-joint muscles, we measured EMG from the biceps long head (BILH) and triceps long head (TRIO). For flexion/extension of the elbow joint, we measured EMG from the brachioradialis (BRAD), triceps lateral head (TRIA), and triceps medial head (TRIM).

### Calibration of Mykin model parameters

To keep the Mykin model parameters from converging on unrealistic values, we restrained their ranges as follows Murray et al. (2000) measured moment arm and muscle architecture in the same cadaveric upper extremity specimens and calculated fascicle lengths. Based on this study, we adjusted the range of moment arm from the minimum (the bony radius elbow = 0.015 m, shoulder = 0.023 m determined from fMRI photographs as shown in Fig. 4*B*) to the maximum (the peak moment arms of the study of Murray et al.). The ranges of the elastic parameters were chosen by considering the static stiffness measured in previous studies (Katayama and Kawato 1993; Tsuji et al. 1995). The bias term *l*_{0} was also constrained so that [*l*_{r}(*u*) − *l*] would not be negative. The range of the length parameter *l*_{1} was assumed by considering the stretch of muscle tendon at the maximal isometric contraction, measured in vivo using ultrasonography (Asakawa et al. 2002; Muramatsu et al. 2001; Reeves and Narici 2003).

After performing *experiment 1*, the 27 trials of *task 1* (3 trials per position) and the 9 trials of *task 2* (1 trial per position) were used to optimize the Mykin parameters for each subject. A total of 172,800 data points (864 s) of normalized EMG and torque data were used for parameter optimization of one subject. A constrained linear regression method was conducted to estimate the parameters ( and *A*_{n}) for each muscle to satisfy *Eq. 6* using MATLAB. The typical parameters for *subject C* are listed in Table 2. The TRIM moment arm of *subject C* seemed to be rather constant, whereas the BRAD moment arm changed. This tendency agrees well with the results reported by Murray et al. (1995). In addition, the moment arms of TRIM and BILH estimated from the optimization were similar to the moment arms (*d1, d2*) measured from MRI images (Fig. 4*B*).

Nine sets of the remaining data were used to test the accuracy of the proposed method. We reproduced time-varying joint torque from the EMG signals and compared them with the calculated joint torque from the force sensor.

## RESULTS

### Estimated joint torque

Numerous studies have investigated the relationship between surface EMG and joint torque. This relation would provide a capable tool for musculoskeletal assessment in various applications, including clinical biomechanics, prosthetics control, and ergonomics. These studies estimated joint torque using a neural network or low-pass filtering. However, for our method, we used a mathematical musculoskeletal tension model. To test the proposed model, we reproduced joint torque from the nine sets of the remaining data. Typical plots of the estimated joint torque during the two types of isometric tasks are shown in Fig. 5. Figure 5*A* shows the results of *task 1* obtained when the hand of *subject C* was in *position 1*. The arrows in the circles above the plots indicate the target directions of the force vectors. The subject generated a 20-N force in each of 16 different target directions according to the direction of the arrow. It was found that flexor muscles (DELC and BRAD) were easy to co-contract with extensor muscles. The shoulder flexor muscle DELC usually showed cross-talk. This co-contraction and S/N ratio of the EMG signals caused some deterioration of torque estimation. Despite that, the estimated torque (solid line) agrees well with the measured torque (dotted line) in all directions because mover agonists (DELS, PMJC, TRIO, BILH, and TRIA) related well with torque. The correlation coefficients based on the entire time series between the measured and estimated torques for each subject are summarized in Table 3. The average and the SD of the correlation coefficients about each position were calculated from 24-s test trials (Fig. 5*A*). The correlation coefficients of *subjects A–E* were lower than those of *subjects F–I*. This tendency was caused by the weak support of the force sensor by the arm of PFM. The correlation coefficients for *position 1* of all subjects were higher than those of other positions since we normalized on the basis of data at *position 1*. The correlation coefficient of *subject D* (another day) was calculated with his previous Mykin model. This result implies that our normalization method based on torque is effective independent of recording instance.

The raw EMG signals of Fig. 5*B* were higher than those of Fig. 5*A*. Despite that, Fig. 5*B* shows that the proposed method is always able to produce zero torque, even if co-contraction levels of muscles are different. As shown in Fig. 6, the maximum difference between the torque zero and the estimated torque was below ±2 Nm. In particular, the results of *subjects B, D*, and *E* were below ±1 Nm while regulating in the range of TCL 4 to TCL 5. The difference also increased with increased TCL, concurrent with reports of nonlinear increases in EMG with muscle force in high muscle activation (Basmajian and De Luca 1985).

The estimated joint torque evidently agrees well with the measured torque in the two force-regulation tasks. Note that, with our simple model, we were able to estimate time-varying torque even for various positions, with the same parameters. The close estimation of joint torque for various positions indicates that the Mykin model can also predict joint stiffness solely with EMG measurement.

### Comparison with the perturbation method

As mentioned in the introductory remarks, Mussa-Ivaldi et al. (1985) developed an experimental perturbation method for measuring stiffness. We used this conventional perturbation method to evaluate the effectiveness of our proposed method. Joint stiffness was calculated from the results of hand stiffness using a Jacobian matrix. For the relaxed state, the average and SD of joint stiffness measured by the PFM for six subjects (*A–E*) were 5.78 ± 2.29 (Rss), 2.51 ± 1.44 (Rse), 1.16 ± 1.9 (Res), and 5.80 ± 2.01 (Ree). In the case of the Mykin model, those were 8.16 ± 1.90 (Rss), 1.71 ± 0.72 (Rse), 1.71 ± 0.72 (Res), and 3.62 ± 0.86 (Ree). The shoulder stiffness (Rss) values were similar to the elbow stiffness (Ree) values from the result of the PFM. In the case of Mykin model, the shoulder stiffness (Rss) values were higher than the elbow stiffness (Ree) values. Two possible explanations could account for these differences. The first is the effect of the physical measurement error of the length and the angle of the arms on the Jacobian computation. The second is that there may have some other contributions by muscles besides those measured in the present experiments. Consequently, these differences produced an effect on endpoint stiffness, as we see from Table 4 that the other elements except the *Kyy* elements are different from those of Mykin model. Even if the elements of stiffness estimated from the Mykin method did not exactly agree with those from the perturbation method, the stiffness ellipses between the two methods do show some strong similarities as shown in Fig. 7. At any given position, the shape and orientation of the ellipses do not change substantially from one subject to another. However, the stiffness magnitudes of *subject E* (female) are smaller than those of the other subjects, which have been well reflected in the Mykin method's estimated stiffness ellipses for *subject E* (Fig. 7). The results obtained from *subject D* on the different days are similar (Fig. 7). The size, shape, and orientation obtained at given positions are also shown in Table 5. The averages, with SD, of the Mykin/PFM ratio of size, shape, and orientation were 1.000 ± 0.268, 0.968 ± 0.351, and 0.995 ± 0.239, respectively (Table 5). The shape results of the proposed method tended to be more circular than those of the perturbation method in *positions 2–4* and more elongated in *positions 7–9*. Also, orientation using the proposed method was more counterclockwise rotated. These differences may have resulted from the accuracy of the torque estimation dependant on the fixation of a force sensor to the arm of PFM, the S/N ratio of EMG signals, and a large range of joint angle. The average and SD of the relative errors of size, shape, and orientation between the perturbation method and the Mykin method were 0.224 ± 0.154, 0.342 ± 0.204, and 0.216 ± 0.138, respectively. Size and orientation characteristics of stiffness ellipses obtained by the proposed method can be summarized as follows: *1*) stiffness ellipses became more anisotropic (elongated) than those of the perturbation method as the hand location approached the distal position and more isotropic (circular) as the hand location moved to the proximal position; *2*) the major axes of stiffness ellipses tended to be roughly oriented toward the shoulder of the subject; in proximal posture, those tended to be more directed toward the elbow than shoulder; and *3*) in a displacement from right to left at constant hand to shoulder distance, stiffness ellipses underwent a counterclockwise rotation.

As shown in Fig. 8, hand stiffness obtained by the two methods increased according to the co-contraction task. The correlation coefficient of *subjects F* and *H* were higher than those of *subjects G* and *I*. The accuracy of stiffness estimation depends on the accuracy of torque estimation. Because the external disturbance was applied only along the *x*+-axis, the stiffness estimation is influenced by the shoulder torque estimation.

## DISCUSSION

Based on the Kelvin-Voigt model, we proposed a mathematical muscle model using a quadratic function for muscle activation. With our simple model, we were able to obtain a close estimation of joint torque for force-regulation tasks and produce zero torque with or without co-contraction. In addition, the proposed method, which uses an equation of differentiated torque, is effective in estimating stiffness. We confirmed that our results correlated with stiffness measured by the perturbation method in various positions. With this study, we wish to establish a basis for the direct estimation of joint torque and stiffness using solely EMG. This is the first report of a method for estimating stiffness ellipses without using perturbation measurements.

### Mykin model based on physiological characteristics of muscles

A model consisting of a two-joint planar arm has been used to simulate shoulder and elbow rotation in the horizontal plane (Gribble and Ostry 2000; Gribble et al. 1998; Katayama and Kawato 1993). Katayama and Kawato used a six-muscle model to simulate virtual trajectories and stiffness ellipse and determined the values of the elastic parameter to match the range of stiffness measured by Bennett et al. (1992). Gribble et al. (1998) proposed a muscle model based on the λ version of the equilibrium-point hypothesis for virtual trajectory simulations. Their model parameters related to stiffness and damping were set to match the measured impedance in static posture (Bennett et al. 1992; Gomi and Kawato 1996; Tsuji et al. 1995). In contrast, our parameters of the presented tension model were optimized to agree with the real torque measured with a force sensor. Consequently, our parameters *k*_{1} in the muscle stiffness *k*(*u*) are 10 times higher than *k*_{1} set by Katayama and Kawato (1993). These parameter differences brought about the size difference in stiffness ellipses. The hand stiffness values of the Mykin model were ∼6–10 times larger than those of the simulation of Katayama and Kawato. However, our length parameters *l*_{1} are 1/4–1/10 times smaller than those of Katayama and Kawato. We set the range of the length parameter *l*_{1} to match the results reported by Reeves and Narici (2003), such that muscle tendon elongation increased to 4.7 ± 1.1 mm at maximal tendon force, determined using ultrasonography. We assumed that the change of muscle length at the maximal isometric contraction resulted from only muscle activation and brought about the stretch of the tendon. Therefore the increased tendon length at the maximal contraction could be considered as the maximal length *l*_{1} parameter.

### Mykin model as an estimator for joint torque

EMG-driven biomechanical models used to examine force have been developed for the ankle (Hof et al. 1987), knee (Lloyd and Besier 2003; Lloyd and Buchanan 1996, 2001; Olney and Winter 1985), shoulder (Laursen et al. 1998), and elbow (Buchanan et al. 1998; Manal et al. 2002). Because direct measures of muscle force in vivo are difficult, these models were typically validated to external joint moments measured using an inverse dynamics approach. These models, applying Hill's model, need expert subjects or preliminary training because they require the subject to generate nearly MVC for normalization. Without training, the MVC obtained by a beginner could be as much as 20–40% less than that of an expert. In addition, it is difficult to generate the same maximum torque each time.

In the Mykin model, it was not important to determine the maximum EMG activation because we only used the force range of the calibration. The proposed normalization technique standardized the minimum EMG values for the generated torque. It was easy for beginner subjects because they were only asked to regulate 20 or 30 N of force for calibration. We attempted to implement a 40-N force calibration task but found it difficult for beginners to generate the proper force. In addition, the accuracy of torque estimation may decrease for forces >50 N because the nonlinearity between EMG and force increases. The accuracy of the proposed method also depends on the S/N ratio of EMG signals, co-contraction, and joint angle range.

Among the related EMG-force based studies, a few have reported on the use of EMG signals as inputs of mathematical muscle models for estimating torque and stiffness in various positions or movements. The previous models that used EMG to match and estimate joint torque and stiffness used conversion parameters (Franklin et al. 2003; Koike and Kawato 1995; Osu and Gomi 1999; Osu et al. 2002). Although these parameters change as the joint angle changes, they assumed them to be constant in different positions or movements. Because these parameters mixed both moment arm and a conversion factor, this assumption of constant moment arm is oversimplified. Kuechle et al. (1997) reported that the moment arm of PMJC might have changed 20% and that of DELS might have changed 30% for the movements. Therefore such assumptions may have led to errors in torque estimation for various positions. The proposed method, by contrast, sets the moment arm at optimized values along the nine positions, which concur with the pattern of the moment arm values measured from MRI images. Consequently, our model was successful at reproducing joint torque and stiffness in the given range of joint angles.

However, our model still has some limitations. The *l*(0) term does not represent the nonlinearity of the force–length relationship and could affect the optimization of the moment arm. Besides, the proposed model does not use the continuous moment arm but only a constant moment arm, which is optimized by posture. Therefore it could cause a reduction in accuracy of estimated torque and stiffness, especially in the range beyond the optimized position or in arm movements. Additionally, recent studies report EMG-driven musculoskeletal models that provide continuously differentiable moment arm curves (Lloyd and Besier 2003; Manal et al. 2002). We expect to incorporate a function of the moment arm curves into a future version of the Mykin model. By including continuously differentiable moment arms in Mykin model, we will be able to model the component of joint stiffness that depends on change in moment arm with change in joint angle, which is an implicit element of joint stiffness defined by *Eq. 7*.

### Comparison with previous stiffness studies

Our stiffness ellipses were similar to stiffness ellipses independently measured with mechanical perturbations. We compared our results with previous observations that reported elbow stiffness during maintained postures. Mussa-Ivaldi et al. (1985) characterized the field of elastic forces as a stiffness ellipse and reported that the elbow stiffness values ranged from 10 to 40 Nm/rad MacKay et al. (1986) found that the elbow stiffness values measured in different subjects ranged from 2 to 12 Nm/rad. Bennett et al. (1992) measured the elbow joint mechanical stiffness during a single joint cyclic movement and reported that elbow stiffness ranged from 2 to 9 Nm/rad. Tsuji et al. (1995) extended the perturbation method to include impedances such as stiffness, viscosity, and inertia, and reported that elbow stiffness ranged from 4 to 13 Nm/rad. They showed that an increase in the handle gripping force led to co-contractions, which increased stiffness. Our results ranged from 2.8 to 5.2 Nm/rad and were in the range of results from previous studies. Discrepancies of perturbation results may have been caused by difference in individuals, experimental setup, and instructions.

### Advantages of the Mykin model

There are two controversial points with regard to the perturbation method. The first is that the perturbation increases net torque. This torque has a large dynamic component caused by the force–velocity relation of the extensor and flexor and has a smaller static component caused by the force–length relationship (Loeb and Ghez 2000). This external force from perturbation itself exerts influence on the maintenance of posture and obstructs the measurement of natural stiffness. The second is that results vary because stiffness changes with muscle activation, and the variance of the results increases as muscle activation increases. Kearney and Hunter (1983) reported that muscle spindle sensitivity seems to change with perturbation amplitude or velocity. Useful results therefore cannot be obtained without performing several experimental trials. These repeated tasks also bring about difficulties for the subjects. The tedious nature of the system causes the subjects to tire and lose concentration. Therefore Osu and Gomi (1999) tried to reconstruct joint stiffness from EMG levels using conversion factors to match the EMG to the PFM measured stiffness. Based on previous studies reporting that muscle activation is closely related to stiffness, Osu et al. (2002) consequently proposed IMCJ to evaluate joint elasticity during movement. These earlier models used parameters or “gains” with no physiological basis, compromising construct validity. We developed a joint and stiffness estimation technique based on physiological characteristics of muscles.

The primary advantage of the proposed method is that it does not require a manipulator for actual stiffness measurements. This creates several benefits. First, it is possible to estimate joint torque and stiffness simultaneously. The added flexibility would be useful in various applications including clinical biomechanics, prosthetics control, and ergonomics assessment. Second, it reduces preparation time for experiments. It also extricates both the operator and the subject from repetitive and tedious tasks. Third, it can be used to estimate arm stiffness on the sagittal plane because its only requirements are a torque sensor and EMG measurements for the parameters’ calibration and to predict an equilibrium point in static posture via a general function consisting of stiffness, determined from EMG, and torque, measured by the torque sensor. However, in this study, we were only able to estimate stiffness of the arm during maintained posture. Future research will expand the Mykin model by a viscosity term and estimate torque and stiffness in movement or co-contraction.

## GRANTS

This work was supported by grants from the Ministry of Education, Science, Sports and Culture, Grant-in-Aid for Scientific Research (B), 16360169, and the JST CREST Program to Y. Sakurai.

## Acknowledgments

We thank S. Masaki (Advanced Telecommunications Research Institute International) and H. Fukuyama (Kyoto University) for the MRI image measurements and R. Osu and M. Kawato (ATR) for the PFM experiments.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2009 the American Physiological Society

## REFERENCES

- Asakawa et al. 2002.↵
- Basmajian and De Luca 1985.↵
- Bennett et al. 1992.↵
- Buchanan et al. 1998.↵
- Burdet et al. 2001.↵
- Clancy et al. 2006.↵
- Clancy and Hogan 1997.↵
- Darainy et al. 2004.↵
- Dolan et al. 1993.↵
- Feldman et al. 1990.↵
- Flanagan et al. 1993.↵
- Flash and Mussa-Ivaldi 1990.↵
- Franklin et al. 2003.↵
- Franklin et al. 2003.
- Ghez 2000.↵
- Gomi and Kawato 1996.↵
- Gomi and Kawato 1997.↵
- Gomi and Osu 1998.↵
- Gottlieb and Agarwal 1971.↵
- Gribble et al. 2003.↵
- Gribble and Ostry 2000.↵
- Gribble et al. 1998.↵
- Hof et al. 1987.↵
- Inman et al. 1952.↵
- Katayama and Kawato 1993.↵
- Kearney and Hunter 1983.↵
- Kearney and Hunter 1990.↵
- Koike and Kawato 1995.↵
- Kuechle et al. 1997.↵
- Laursen et al. 1998.↵
- Lloyd and Besier 2003.↵
- Lloyd and Buchanan 1996.↵
- Lloyd and Buchanan 2001.↵
- Loeb and Ghez 2000.↵
- MacKay et al. 1986.↵
- Manal et al. 2002.↵
- Maton et al. 1987.↵
- Muramatsu et al. 2001.↵
- Murray et al. 2000.↵
- Murray et al. 1995.↵
- Mussa-Ivaldi et al. 1985.↵
- Olney and Winter 1985.↵
- Osu et al. 2002.↵
- Osu and Gomi 1999.↵
- Özkaya and Nordin 1991.↵
- Perreault et al. 2001.↵
- Perreault et al. 2002.↵
- Perreault et al. 2004.↵
- Prilutsky 2000.↵
- Reeves and Narici 2003.↵
- Tsuji et al. 1995.↵