2178 IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 58, NO. 8, AUGUST 2011 A Model of the Surface Electromyogram in Pathological Tremor Jakob Lund Dideriksen, Member, IEEE, Roger M. Enoka, and Dario Farina*, Senior Member, IEEE Abstract—The study developed a novel multiscale model for simulating the surface electromyogram (EMG) of an antagonistic pair of muscles during pathological tremor. By combining and expanding mathematical descriptions from motor units to limb kinematics, the model constitutes the first attempt to simulate the surface EMG and the individual motor unit activity under the influence of descending voluntary command, oscillatory noise in the descending signal, and afferent feedback when controlling a freely moving limb to achieve a predefined angular trajectory. The oscillatory noise was adjusted to simulate various types of pathological tremor. The simulations replicated previously reported experimental results for the power spectral density of the surface EMG, the angular velocity of the limb, and single motor unit activity. The model provides a powerful tool for extracting information about how the surface EMG can be used to describe tremor in various conditions, including different tremor frequencies and intensities, that cannot be achieved solely with experimental approaches. Index Terms—Closed-loop systems, electromyography, neuromuscular model, tremor. I. INTRODUCTION ATHOLOGICAL tremor is an involuntary rhythmic oscillation of a body part that emerges with the development of some neurological conditions, such as Parkinson’s disease and cerebellar dysfunction [1]. Tremor arises from oscillations in the descending drive to the muscles from the motor cortex [2], and is maintained and enhanced by stretch reflexes [3] and the biomechanical properties of the limbs [4]. Current methods for the treatment of tremor include medication, deep brain stimulation [5], functional electrical stimulation (FES) [6], [7], and orthotic exoskeletons [8]. These methods, however, are either not effective for all tremor patients, or too impractical for daily use. In the absence of an effective treatment, there is a need to enhance our knowledge on tremor, so that new approaches can be developed to manage the disorder. P Manuscript received July 24, 2010; revised January 5, 2011; accepted February 1, 2011. Date of publication February 22, 2011; date of current version July 20, 2011. This work was supported in part by the European Union (EU) project TREMOR. Asterisk indicates corresponding author. J. L. Dideriksen is with the Center for Sensory-Motor Interaction, Department of Health Science and Technology, Aalborg University, Aalborg, Denmark 9220 (e-mail: jldi@hst.aau.dk). R. M. Enoka is with the Department of Integrative Physiology, University of Colorado, Boulder, CO 80309 USA (e-mail: roger.enoka@colorado.edu). *D. Farina is with the Department of Neurorehabilitation Engineering, Bernstein Center for Computational Neuroscience, Georg-August University of Göttingen, Göttingen, Germany 37075 (e-mail: dario.farina@bccn. uni-goettingen.de). Digital Object Identifier 10.1109/TBME.2011.2118756 Computational models are valuable tools to analyze interdependences among involved elements in physiological systems, including the genesis of pathological tremor. For example, a model of tremor based on a Hill-type muscle model has been developed to examine the contribution of afferent feedback to centrally generated tremor [3]. Similarly, other modeling studies have addressed the influence of proprioceptive reflexes and the gain in these reflex pathways on motor stability [9]–[13]. Although most studies of pathological tremor measure the activity with inertial sensors, surface EMG recordings can provide a more detailed description of tremor. For example, it is difficult to identify the muscle from which the tremor arises by using inertial sensors as the oscillations can be transferred across joints and coactivation of noninvolved muscles may mask the source of the tremor. The aim of this study was to develop a multiscale model of pathological tremor capable of simulating both limb dynamics, and surface EMG signals generated during tremor. The approach was to expand and combine mathematical descriptions of motor unit activity through to limb kinematics into one integrative model. The neuromuscular model of isometric force production proposed by Fuglevand et al. [14] was expanded to describe dynamic contractions, and combined with models of afferent feedback and a model that generates surface motor unit action potentials EMG. Tremor was imposed as an oscillatory noise embedded in the simulated descending drive to a pair of antagonistic muscles that controlled the angular displacement of a limb about a joint. Representative simulations are presented to demonstrate the influence of tremor on angular displacement and surface EMG. II. METHODS Fig. 1 depicts the structure of the model for one of the two antagonistic muscles (first dorsal interosseus muscle (FDI) and second palmar interosseus) that were modeled for the control of index finger abduction/adduction [15]. Although the model and the simulations are presented for one muscle pair, the model can be adjusted for any other muscle group and joint system. The muscle parameters used in this study are reported in Table I, and other model parameters are listed in Table II. Muscle parameters for the FDI (Table I) were obtained from anatomical and experimental studies [16], [17], except for the relations between muscle and tendon length and between musculotendon length and moment arm (distance from joint to the tendon attachment on the index finger), which were estimated from the anatomy of the index finger. The other model parameters (Table II) were either adopted or derived from the original description as explained in the following section. The two antagonistic muscles 0018-9294/$26.00 © 2011 IEEE DIDERIKSEN et al.: MODEL OF THE SURFACE ELECTROMYOGRAM IN PATHOLOGICAL TREMOR 2179 TABLE II COEFFICIENT VALUES Fig. 1. In response to net synaptic input, the motor neuron pool generates trains of action potentials for activated motor units (MU). The MU action potentials serve as input to a model that simulates a surface EMG signal and another model that activates the muscle to produce a net muscle force to control joint angle. The muscle force signal activates Golgi tendon organs that contribute a feedback signal (Ib afferent) to the net synaptic input. Similarly, angular displacements and velocities activate muscle spindles that contribute (Ia afferent) to the net synaptic input. The PID algorithm receives input that comprises the difference between an intended and actual angular trajectory and produces an output that is denoted as the descending drive. Asterisks indicate submodels that run concurrently in both antagonistic muscles, whereas the other submodels run separately for each muscle. PID: proportional, integrative, derivative (see description in the text). TABLE I PHYSIOLOGICAL MODEL PARAMETERS motor neurons, along with a term describing the synaptic noise NSI(t) = DD(t) + KTO · TO(t) + KGTO · GTO(t − dGTO ) + KM SP · MSP(t − dM SP ) + FSN(t) where DD denotes the descending drive from supraspinal centers, GTO and MSP indicate the afferent input from Golgi tendon organs and muscle spindles, respectively, KTO , KGTO , and KM SP are the gain factors, and TO denotes the tremor oscillations. Unless stated otherwise, the gains were set to 10, so that the afferent feedback provided ∼30% of the net synaptic input, as has been observed experimentally during isometric contractions [18]. Both types of afferent input were implemented with a delay (dGTO and dM SP ) to account for the conduction time from the muscle to the spinal cord. DD was estimated using a proportional–integral derivative (PID) control algorithm, described by the following equation: DD(t) = Kp · AE(t) + Ki · AE(t) + Kd · AE(t) were assumed to be similar, and the tremor was imposed on both muscles with similar amplitudes in opposite phases. The positive and negative output of the control algorithm was used as input to the agonist and antagonist muscles, respectively. The joint kinematics were estimated from the net force exerted by the two muscles. A. Net Synaptic Input to Motor Neurons The input to the motor neuron pool model comprised the sum of inputs from four sources, and represents the spatial and temporal summation of the net synaptic inputs received by the (1) (2) where Kp , Ki , and Kd are the gain factors (see Table II for values). The values of these gains were set to produce experimentally observed variability during constant force contractions in normal conditions. AE indicates the joint angular error and corresponded to the differences between the target joint angle and the angle simulated by the model. AE was updated at 3 Hz and averaged with a moving-average filter that had a length equal to one period of the imposed tremor frequency. This specification resulted in the PID algorithm compensating for the average angular error and not for the tremor oscillations. The PID control algorithm has been used previously to simulate the control of muscle force by the central nervous system [19]. An oscillatory noise (TO), representing the central origin of the tremor, was added to the descending drive. The TO input was implemented as band-pass filtered (pass-band width: 1 Hz around a predefined tremor frequency) white Gaussian noise with an amplitude between −0.5 and 0.5. The intensity of the tremor was controlled by KTO , whose value thereby represented the peak-to-peak amplitude of the descending oscillations. The frequency of pathological tremor is usually in the range of 3– 13 Hz [20]. 2180 IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 58, NO. 8, AUGUST 2011 GTO acted as a sensor of the force applied to the tendons, and was modeled by the following transfer function, as earlier described [21]: (1 + (s/(0.15)))(1 + (s/(0.15)))(1 + (1/16)) . (1 + (s/(0.2)))(1 + (s/2))(1 + (s/37)) (3) MSP denoted the change and rate of change in muscle length. The proposed model implements the simplified version of a previously proposed model [22], according to which the average intrafusal muscle fiber experiences similar dynamics as the whole muscle HGTO (s) = MSP(t) = lm · v̄l0.3 (4) where lm denotes normalized muscle length and v̄l represents normalized muscle velocity. The spindle gain was represented by KM SP (1). Although the different types of synaptic input correspond to the summed excitatory and inhibitory postsynaptic potentials, they were represented in the model by a constant value. As the summation of the postsynaptic potentials results in a fluctuating membrane potential, a random value from a Gaussian distribution with zero mean and standard deviation proportional to the magnitude of the total synaptic input was added to the net synaptic input [19]. B. Motor Neuron Population and Muscle The models of the motor neuron population and the muscle were adopted from those developed by Fuglevand et al. [14]. The inputs to both models comprised the motor unit action potentials that were characterized with equations describing the motor unit recruitment, rate coding, and saturation properties. The motor unit activity enabled a model to simulate the total active force (Fact ) produced by the muscle. The twitch forces simulated for isometric contractions with the Fuglevand model were adjusted according to the force-length and force-velocity relations described by Zajac [23]. C. Joint Kinematics The passive viscoelastic joint force was represented with terms for the passive properties of muscle, tendons, and other tissues, as well as joint surface friction [24]. The expression comprised one term for the elastic force (Fe ) and one term for the hysteresis (Fh ) Fve (t) = a · (Fe (t) + Fh (t)) (5) where a denotes the gain that constrained the maximum achievable angle to the defined range of motion (see Tables I and II for values). Fe and Fh were defined as Fe (t) = b0 · (eb 1 ·∆ θ (t) − eb 2 ∆ θ (t) ) (6) ¯ ¯ ¯ Fh (t) = sign(θ̇(t)) · (c1 + c2 · θ̇2 ) · |θ̇(t)|n , (7) these terms (Table II) were adopted from the original implementation [24]. Limb displacement was produced by a shortening contraction of the muscle, which was assumed to be proportional to the total force applied to the muscle and was expressed as a percentage of the maximal voluntary contraction (MVC) force. Fload refers to the force required to support a predefined added inertial load, as described by the equation Ftot (t) = Fact (t) + Fve (t) + Fload . (8) The relation between acceleration and the normalized force was described by the equation accm (t) = accm ax · Ftot (t) (9) where accm ax indicates the maximal acceleration the muscle can achieve during a shortening contraction, and F denotes the force expressed as a fraction of the maximal voluntary force. The maximal acceleration was derived from experimental data for ballistic contractions with FDI [17]. Changes in muscle velocity and length were obtained by single and double integration of the acceleration. The relation between the change in musculotendon length and the change in limb angle was derived from the law of cosines by the equation   2 2 2 lmt − lm a − (1 − lm a ) (10) θ(t) = 180 − arccos −2 · lm a · (1 − lm a ) where lmt denotes the musculotendon length normalized to its slack length, and lm a indicates the moment arm expressed as a fraction of the slack length (see Table I for values). Tendon length was assumed to be constant, and was expressed as a fraction of muscle slack length. D. EMG Generation The surface EMG signal was simulated with a model that comprised a multilayer cylindrical volume conductor [25], using the same parameter settings as in Keenan et al. [26]. The shape of each simulated intracellular action potential depended on the instantaneous motor unit conduction velocity, which corresponded to the velocity at which the action potential propagated along the muscle fibers. The conduction velocity of each motor unit was implemented as a function of the instantaneous discharge rate to represent the velocity-recovery function of the muscle fibers [27]. The surface EMG signals were sometimes simulated ten times for each condition by randomly assigning the locations of the muscle fibers within the muscles to replicate the variability observed for measurements on different subjects. E. Experimental Recording where ∆θ denotes the normalized difference in instantaneous ¯ angle and slack angle, θ̇ represents the normalized angular velocity, and b0 , b1 , b2 , c1 , c2 , n are coefficients. The values for The simulated data were compared with experimental measurements of surface EMG and limb movement from one woman (66 years) with essential tremor. EMG was recorded from the wrist extensors and movement was detected by using an DIDERIKSEN et al.: MODEL OF THE SURFACE ELECTROMYOGRAM IN PATHOLOGICAL TREMOR 2181 Fig. 3. Force (a) and joint angle (b) during the first 5 s of an 8◦ abduction movement with (thin line) and without (bold line) tremor. The dashed lines in (a) represents the net muscle force. Force for the antagonist muscle is represented in negative values. Force for the antagonist muscle was equal to zero during the simulation without tremor. Fig. 2. Metacarpophalangeal joint angle, agonist force, agonist EMG, and the discharge times of the smallest motor unit as the muscle exerted an abduction force to achieve a static target angle of 5◦ and was subjected to 5-Hz (a) and 10-Hz (b) tremor, respectively. accelerometer placed at the wrist, while the patient held her arms outstretched against gravity holding a mass of 500 g. III. RESULTS The simulation results reported in this paper were obtained from voluntary contraction forces in the range 0 to 20% of the maximal voluntary contraction force to match the forces that occur most often during activities of daily living. The maximal level of tremor was set at a peak-to-peak angle of approximately 9◦ (obtained by setting KTO = 20), which corresponds to substantial tremor [28]. In this case, the oscillation amplitude was equivalent to ∼45% of the excitation level required to produce maximal force. Simulated tremor frequencies ranged from 5–10 Hz, which covers most of the range of typical tremor frequencies [1]. Fig. 2 depicts the simulated angle, force and EMG for the agonist muscle, and discharge times of the smallest motor unit during an abduction with a static target angle of 5◦ that included moderate tremor (KTO = 15) at 5 and 10 Hz. Although the amplitude of the imposed oscillations was similar for both frequencies, the fluctuations in angle and force were significantly more pronounced for the lower-frequency tremor. The greater tremor amplitude for the lower frequency can be explained by more action potentials, and thus, more force with longer tremor intervals. The longer tremor intervals also included an increased tendency for motor units to discharge multiple action potentials in each tremor period. Fig. 3 shows the simulated force for both the agonist and antagonist muscles and the metacarpophalangeal joint angle in conditions with and without tremor during an 8◦ abduction. The tremor was set at 8 Hz with KTO = 15. When tremor was present in both muscles, agonist muscle force had to be greater than in the no-tremor condition to achieve the desired angle. The greater Fig. 4. The smoothed power spectral density for angular velocity (a), (c), (e) and agonistic EMG (b), (d), (f) during a contraction with a static target angle of 0◦ with (8 Hz) or without tremor and inertial loads of 0% MVC (a), (b), 10% MVC (c), (d), and 20% MVC force (e), (f). The solid line indicates the condition without tremor, whereas the dashed and dotted lines indicate the conditions with tremor (8 Hz) with amplitudes (KT O ) of 10 and 20, respectively. agonist muscle force involved greater modulation of motor unit activity by the descending oscillations, thereby producing force fluctuations of a higher magnitude than the antagonist muscle. Fig. 4 indicates the power spectral density of the angular velocity and agonist EMG during simulations with a static target angle of 0◦ , inertial loads of 0%, 10%, and 20% of MVC force, and values for KTO of 10 and 20 in the presence (8 Hz) and absence of tremor. The simulated power spectral density for the EMG is similar to that observed experimentally [29], [30]. The magnitude of the spectral peak at the tremor frequency increased with load, but tended to remove the harmonic peaks for the EMG. The appearance of a low-frequency spectral peak for angular velocity (Fig. 4e) has been observed experimentally [4]. Fig. 5 compares the power spectrum of the EMG and histograms of the distribution of intervals between limb 2182 IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 58, NO. 8, AUGUST 2011 Fig. 6. The relation between EMG burst amplitude above baseline and peakto-peak angle during contractions with static target angles of 0◦ , inertial loads of 0% MVC (circles), 10% MVC (squares), and 20% MVC (triangles), and with tremor amplitudes of 10 au (no fill), 15 au (gray), and 20 au (black). Error bars indicate the variability in the EMG amplitude induced by different locations of the activated fibers in the muscle. Fig. 5. The smoothed power spectrum of the EMG (a), (b) and the interburst interval histograms (c), (d) in experimental (a), (c) and simulated (b), (d) conditions. Both simulated and experimental data were from a 4-s interval. In the experiment the subject held her arms outstretched against gravity and a added mass of 500 g, whereas the simulation was performed with an inertial load of 20% MVC and a tremor amplitude of 4 au. Tremor frequency of both conditions was approximately 6 Hz. acceleration for the simulated and experimental conditions. The simulation was performed with a static target angle of 0◦ , 20% MVC inertial load, a tremor frequency of 6 Hz, and tremor amplitude (KTO ) of 4 au. The power at the spectral peak around the tremor frequency (±1 Hz) expressed as a percentage of the total power was similar in the two examples (experimental: 1.1%, simulated: 1.2%). Furthermore, the distributions of the intervals between acceleration peaks were similar as expressed by the kurtosis (experimental: 8.4 au, simulated 7.6 au). Fig. 6 plots the relation between the average rectified agonist EMG amplitude during tremor-bursts and the peak-to-peak angle of the tremor (8 Hz) with tremor amplitudes (KTO ) of 10, 15, and 20, and with inertial loads of 0%, 10%, and 20% MVC force for a static target angle of 0◦ . The angle was bandpass filtered (band: 7–9 Hz, second-order Butterworth filter) to isolate the angular displacements due to tremor, and exclude those due to variability in the descending drive. In general, the EMG burst amplitude and peak-to-peak angle were correlated. However, the EMG amplitude was close to zero for the heaviest load (20% MVC) and tremor amplitudes (KTO ) of 15 and 20 despite the presence of significant 8-Hz oscillations. This result indicates that the amplitude of the agonist EMG above baseline is not always a reliable index of the magnitude of the tremor fluctuations. The cancellation of EMG amplitude due to summation of single motor unit potentials [26] did not vary across conditions (∆ cancellation <4%). Fig. 7 comprises the first- and second-order interspike interval (ISI) histograms for all active motor units during contractions with static target angles of 0◦ , 8-Hz tremor, KTO set to 15, and inertial loads of 0%, 10%, and 20% MVC force. Fig. 7. First- (a), (c), (e) and second-order (b), (d), (f) ISI histograms for all motor units during a static contraction with a target angle of 0◦ with inertial loads of 0% MVC (a), (b), 10% MVC (c), (d), and 20% MVC (e), (f). The dashed vertical lines in (a), (c), (e) indicates the half and full tremor period (125 ms), whereas the dashed vertical lines in (b), (d), (f) indicates the full and double tremor period. The first-order histogram indicates the average temporal distance between successive motor unit discharges, whereas the second-order histogram denotes the average temporal distance between every second motor unit discharge. The histograms indicate an increased tendency for double discharges in each tremor period with an increase in load. The close similarity between the first- and second-order histogram suggests that each motor unit exhibits similar trends throughout the contraction and does not tend to shift between single and double discharges DIDERIKSEN et al.: MODEL OF THE SURFACE ELECTROMYOGRAM IN PATHOLOGICAL TREMOR Fig. 8. The smoothed power spectral density for angular velocity (a), (c) and agonist EMG (b), (d) during a static contraction at a target angle of 0◦ with inertial load of 10% MVC, either 5-Hz (a), (b) or 10-Hz tremor (c), (d), and afferent gains of 0 (solid line), 20 (dashed line), and 40 (dotted line). in each tremor period. These observations are consistent with experimental findings [30]. Fig. 8 compares the influence of different afferent gains on the power spectral density of the angular velocity and the agonist EMG during contractions with a static target angle of 0◦ with 10% MVC added load in conditions with imposed 5- and 10-Hz tremor, both with KTO set to 15. The afferent gains [KGTO , KM SP , and (1)] were set to 0, 20, and 40, which were equivalent to ratios for the contributions by descending drive and total afferent input of approximately 1:0, 1:0.66, and 1:1.33 in the simulated conditions. The power spectral densities show that the power of the tremor frequency was enhanced by more afferent input, whereas the baseline did not change with afferent input. Increases in afferent feedback also seemed to increase the magnitude of the second- and third-harmonic relative to the peak in the tremor frequency. Remarkably, there was no spectral peak at the tremor frequency for the EMG but only harmonics peaks when the afferent feedback was greatest during the 5-Hz tremor. IV. DISCUSSION The work described in the current report was based on the combination and expansion of seven previously published computational models that involved anatomical, biomechanical, neuromuscular, and physiological descriptions of muscle activity. Such an approach provided a significant advance in the development of a detailed multiscale model of pathological tremor to examine the tremor-associated changes in motor unit activity [31] and the EMG amplitude [29], [30]. Furthermore, the model can be easily applied to study human movement in nonpathological conditions by removing the central oscillations (setting amplitude of TO = 0). The representative simulation results were consistent with a number of experimental observations on tremor, including the power spectral densities for joint angle and surface EMG, single motor unit discharge, and the characteristics of an in- 2183 dividual with essential tremor (Fig. 5). Furthermore, the results (Fig. 8) confirm and extend the findings that afferent input contributes significantly to tremor amplitude for both angular velocity and surface EMG [3]. Although the representative simulations mainly provided a qualitative validation of the model, the multiscale structure of the model provides the foundation for more extensive examinations of different model parameters in future studies. As described previously, tremor accompanies a number of pathological conditions and exhibits variable characteristics, including intensities, frequencies, and conditions in which it occurs. Accordingly, tremor is usually classified as resting, postural, or action tremor [1]. The proposed model is able to simulate both tremor of varying intensities and frequencies, and these different classifications of tremor. For example, action tremor, which occurs during voluntary contractions, can be simulated by implementing the tremor amplitude (1) as a function that is inversely related to the amplitude of the descending drive, and the converse approach can be used to study resting tremor. Due to the complexity and multiscale structure of the model, there are a number of limitations and simplifications that need to be acknowledged when interpreting the results. One limitation involved the motor neuron model, which was derived from experimental observations on slow isometric contractions. The model, therefore, was unable to simulate the experimental observations of three single motor unit discharges in each tremor period with ISIs <40 ms [31], even though the main trend of the ISI histograms (Fig. 7) was consistent with the experimental results. In the current version of the motor neuron model, the discharge time # n + 1 was estimated from the excitation level at the time of discharge # n and did not accommodate rapid changes in excitation level, as occurs during tremor. Furthermore, the model does not account for changes in the location of the surface EMG electrode relative to the muscle fibers, as occurs during movements. Although such changes in the recording system influence the shape of each single motor unit action potential as they appear on the surface of the skin [32], the changes were assumed to be negligible and unlikely to change the compound surface EMG signal significantly as the simulated changes in angle were relatively small (<4◦ ). Finally, tendon length was assumed to be constant in the model, even though tendons have a nonlinear stress-strain curve [22]. Despite these limitations, the model was able to reproduce the features of the surface EMG signals recorded during tremor. Furthermore, the relative significance of these various limitations can be determined, for example with a sensitivity analysis. Earlier, tremor-detection algorithms based on inertial data have been proposed for a tremor-suppression system that uses FES [33]. The muscular activity underlying tremor, however, is described more accurately with surface EMG signals than with recordings obtained from inertial sensors. Therefore, an EMGbased tremor detection system might provide a more effective FES-suppression of tremor, thereby minimizing the patient discomfort and muscle fatigue that are inevitably associated with FES. The proposed model could be extended to develop such algorithms and test them in various tremor conditions, with a more 2184 IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 58, NO. 8, AUGUST 2011 systematic approach than can be achieved by relying solely on experimental recordings. In conclusion, this study presents a novel multiscale model of pathological tremor that is able to simulate surface EMG signals consistent with experimental observations of pathological tremor. REFERENCES [1] G. Deuschl, P. Bain, M. Brin, Y. Agid, L. Benabid, R. Benecke et al., “Consensus statement of the movement disorder society on tremor,” Movement Disorders, vol. 13, no. Suppl. 3, pp. 2–23, 1998. [2] J. Raethjen, R. B. Govindan, F. Kopper, M. Muthuraman, and G. Deuschl, “Cortical involvement in the generation of essential tremor,” J. Neurophysiol., vol. 97, no. 5, pp. 3219–3228, 2007. [3] D. Zhang, P. Poignet, A. P. Bo, and W. T. Ang, “Exploring peripheral mechanism of tremor on neuromusculoskeletal model: A general simulation study,” IEEE Trans. Biomed. Eng., vol. 56, no. 10, pp. 2359–2369, Oct. 2009. [4] M. E. Héroux, G. Pari, and K. E. Norman, “The effect of inertial loading on wrist postural tremor in essential tremor,” Clin. Neurophysiol., vol. 120, no. 5, pp. 1020–1029, 2009. [5] D. E. Vaillancourt, M. M. Sturman, L. V. Metman, R. A. E. Bakay, and D. M. Corcos, “Deep brain stimulation of the VIM thalamic nucleus modifies several features of essential tremor,” Neurology, vol. 61, no. 7, pp. 919– 925, 2003. [6] A. Prochazka, J. Elek, and M. Javidan, “Attenuation of pathological tremors by functional electrical stimulation I: Method,” Ann. Biomed. Eng., vol. 20, no. 2, pp. 205–224, 1992. [7] M. Javidan, J. Elek, and A. Prochazka, “Attenuation of pathological tremors by functional electrical stimulation II: Clinical evaluation,” Ann. Biomed. Eng., vol. 20, no. 2, pp. 225–236, 1992. [8] M. Manto, E. Rocon, J. Pons, J. M. Belda, and S. Camut, “Evaluation of a wearable orthosis and an associated algorithm for tremor suppression,” Physiol. Meas., vol. 28, no. 4, pp. 415–425, 2007. [9] I. Fukumoto, “Computer simulation of Parkinsonian tremor,” J. Biomed. Eng., vol. 8, no. 1, pp. 49–55, 1986. [10] M. N. Oguztoreli and R. B. Stein, “The effects of multiple reflex pathways on the oscillations in neuro muscular systems,” J. Math. Biol., vol. 3, no. 1, pp. 87–101, 1976. [11] A. Prochazka, D. Gillard, and D. J. Bennett, “Implications of positive feedback in the control of movement,” J. Neurophysiol., vol. 77, no. 6, pp. 3237–3251, 1997. [12] M. Santillán, R. Hernández-Pérez, and R. Delgado-Lezama, “A numeric study of the noise-induced tremor in a mathematical model of the stretch reflex,” J. Theor. Biol., vol. 222, no. 1, pp. 99–115, 2003. [13] R. B. Stein and M. N. Oguztoreli, “Tremor and other oscillations in neuromuscular systems,” Biol. Cybern., vol. 22, no. 3, pp. 147–157, 1976. [14] A. J. Fuglevand, D. A. Winter, and A. E. Patla, “Models of recruitment and rate coding organization in motor-unit pools,” J. Neurophysiol., vol. 70, no. 6, pp. 2470–2488, 1993. [15] Z. Li, H. J. Pfaeffle, D. G. Sotereanos, R. J. Goitz, and S. L. Woo, “Multi-directional strength and force envelope of the index finger,” Clin. Biomech., vol. 18, no. 10, pp. 908–915, 2003. [16] M. D. Jacobsobn, R. Raab, B. M. Fazelli, R. A. Abrams, M. J. Botte, and R. L. Lieber, “Architectural design of the human intrinsic hand muscles,” J. Hand Surg. [Am]., vol. 17, pp. 804–809, 1992. [17] S. Baudry and J. Duchateau, “Postactivation potentiation in a human muscle: Effect on the rate of torque development of tetanic and voluntary isometric contractions,” J. Appl. Physiol., vol. 102, no. 4, pp. 1394–1401, 2007. [18] V. G. Macefield, S. C. Gandevia, B. Bigland-Ritchie, R. B. Gorman, and D. Burke, “The firing rates of human motoneurones voluntarily activated in the absence of muscle afferent feedback,” J. Physiol. (Lond.), vol. 471, pp. 429–443, 1993. [19] J. L. Dideriksen, D. Farina, M. Bækgaard, and R. M. Enoka, “An integrative model of motor unit activity during sustained submaximal contractions,” J. Appl. Physiol., vol. 108, pp. 1550–1562, 2010. [20] M. Hallett, “Overview of human tremor physiology,” Movement Disorders, vol. 13, no. Suppl. 3, pp. 43–48, 1998. [21] A. Prochazka and M. Gorassini, “Ensemble firing of muscle afferents recorded during normal locomotion in cats,” J. Physiol. (Lond.), vol. 507, no. 1, pp. 293–304, 1998. [22] J. C. Houk, W. Z. Rymer, and P. E. Crago, “Dependence of dynamic response of spindle receptors on muscle length and velocity,” J. Neurophysiol., vol. 46, no. 1, pp. 143–166, 1981. [23] F. E. Zajac, “Muscle and tendon: Properties, models, scaling, and application to biomechanics and motor control,” Crit. Rev. Biomed. Eng., vol. 17, no. 4, pp. 359–411, 1989. [24] A. Esteki and J. M. Mansour, “An experimentally based nonlinear viscoelastic model of joint passive moment,” J. Biomech., vol. 29, no. 4, pp. 443–450, 1996. [25] D. Farina, L. Mesin, S. Martina, and R. Merletti, “A surface EMG generation model with multilayer cylindrical description of the volume conductor,” IEEE Trans. Biomed. Eng., vol. 51, no. 3, pp. 415–426, Mar. 2004. [26] K. G. Keenan, D. Farina, R. Merletti, and R. M. Enoka, “Amplitude cancellation reduces the size of motor unit potentials averaged from the surface EMG,” J. Appl. Physiol., vol. 100, no. 6, pp. 1928–1937, 2006. [27] J. L. Dideriksen, D. Farina, and R. M. Enoka, “Influence of fatigue on the simulated relation between the amplitude of the surface electromyogram and muscle force,” Philos. Trans. A Math. Phys. Eng. Sci., vol. 368, no. 1920, pp. 2765–2781, Jun. 2010. [28] D. Duval, A. F. Sadikot, and M. Panniset, “The detection of tremor during slow alternating movements performed by patients with early Parkinson’s disease,” Exp. Brain. Res., vol. 154, pp. 395–398, 2004. [29] J. Marusiak, A. Jaskólska, K. Kisiel-Sajewicz, G. H. Yue, and A. Jaskólski, “EMG and MMG activities of agonist and antagonist muscles in Parkinson’s disease patients during absolute submaximal load holding,” J. Electromyogr. Kinesiol., vol. 19, no. 5, pp. 903–914, 2009. [30] J. Raethjen, M. Lindemann, H. Schmaljohann, R. Wenzelburger, G. Pfister, and G. Deuschl, “Multiple oscillators are causing Parkinsonian and essential tremor,” Movement Disorders, vol. 15, no. 1, pp. 84–94, 2000. [31] C. N. Christakos, S. Erimaki, E. Anagnostou, and D. Anastasopoulos, “Tremor-related motor unit firing in Parkinson’s disease: Implications for tremor genesis,” J. Physiol. (Lond.), vol. 587, no. 20, pp. 4811–4827, 2009. [32] D. Farina, R. Merletti, M. Nazzaro, and I. Caruso, “Effect of joint angle on EMG variables in leg and thigh muscles,” IEEE Eng. Med. Biol. Mag., vol. 20, no. 6, pp. 62–71, Nov./Dec. 2001. [33] L. Z. Popović, T. B. Šekara, and M. B. Popović, “Adaptive band-pass filter (ABPF) for tremor extraction from inertial sensor data,” Comput. Methods Programs Biomed., vol. 99, no. 3, pp. 298–305, Sep. 2010. Jakob Lund Dideriksen (M’10) received the M.Sc. degree in biomedical engineering, in 2009, from the Center for Sensory Motor Interaction, Department of Health Science and Technology, and is currently working toward the Ph.D. degree on the characteristics of neural signals in pathological tremor and the use of these signals as the basis for automatic tremor suppression systems at The International Doctoral School in Biomedical Science and Engineering, Aalborg University, Aalborg, Denmark. His research interests include motor control, muscle fatigue, functional electrical stimulation for rehabilitation purposes, and computational modeling of biomedical signals related to these areas. DIDERIKSEN et al.: MODEL OF THE SURFACE ELECTROMYOGRAM IN PATHOLOGICAL TREMOR Roger M. Enoka received the undergraduate training in physical education at the University of Otago, Dunedin, New Zealand, in 1970 the M.S. degree in biomechanics, and Ph.D. degree in kinesiology from the University of Washington in Seattle, in 1976 and 1981, respectively. He has held faculty positions in the Department of Exercise and Sport Sciences and the Department of Physiology at the University of Arizona, Tucson, and in the Department of Biomedical Engineering at the Cleveland Clinic Foundation. He is currently a Professor in the Department of Integrative Physiology at the University of Colorado, Boulder. His research interests include neuromuscular mechanisms that mediate the acute adjustments and chronic adaptations experienced by the neuromuscular system of humans. 2185 Dario Farina (M’01–SM’09) received the M.Sc. degree in electronics engineering from Politecnico di Torino, Torino, Italy, in 1998, and the Ph.D. degree in automatic control and computer science and in electronics and communications engineering from the Ecole Centrale de Nantes, France, and Politecnico di Torino, respectively, in 2002. During 2002–2004, he was a Research Assistant Professor at Politecnico di Torino, and during 2004– 2008, an Associate Professor of biomedical engineering at Aalborg University, Aalborg, Denmark. From 2008 to 2010, he was Full Professor of motor control and biomedical signal processing and Head of the Research Group on Neural Engineering and Neurophysiology of Movement at Aalborg University. In 2010, he was appointed Full Professor and Founding Chair of the Department of Neurorehabilitation Engineering at the University Medical Center Göttingen, Georg-August University, Germany, within the Bernstein Center for Computational Neuroscience. His research focuses on biomedical signal processing and modeling, neurorehabilitation, and neural control of movement. Within these areas, he has authored or coauthored approximately 200 papers in peer-reviewed journals and about 300 among conference papers/abstracts, book chapters and encyclopedia contributions. He is an Associate Editor of the journal Medical & Biological Engineering & Computing, and a member of the editorial boards of the Journal of Electromyography and Kinesiology and of the Journal of Neuroscience Methods. He is also an Associate Editor of the IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. He is the recipient of the 2010 IEEE Engineering in Medicine and Biology Society Early Career Achievement Award for his contributions to biomedical signal processing and to electrophysiology. Prof. Farina is the Vice-President of the International Society of Electrophysiology and Kinesiology (ISEK).