INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Recently, a new optimal control model has been proposed that captures several invariant features of human movements by making a single physiological assumption. The Task Optimization in the Presence of Signal-dependent noise model (TOPS; Hamilton and Wolpert 2002; Harris and Wolpert 1998) assumes that there is noise in the motor command and that the amount of noise scales with the motor command's magnitude. In the presence of such noise, the same sequence of intended motor commands, if repeated many times, will lead to a probability distribution of position at the end of the movements. Modifying the sequence of motor commands can control aspects of this probability distribution. In the TOPS model, the task specifies how aspects of the distribution are penalized, and this forms the cost. For example, in a simple aiming movement, the task is to minimize the final error, as measured by the variance. Provided the signal-dependent noise in the motor command has a constant coefficient of variation (CV), i.e., scales linearly with the magnitude of the motor command, this model accurately predicts the invariant kinematics of arm movement trajectories, Fitts' law, and the two-thirds power law (Fitts 1954; Lacquaniti et al. 1983; Morasso 1981). The essential feature of TOPS is that by default the motor command to the muscles includes physiological noise and that the SD of this noise is proportional to the magnitude of the motor command.

Isometric contractions of the hand muscles exhibit variability in force production that is proportional to the mean force exerted (Enoka et al. 1999; Galganski et al. 1993; Laidlaw et al. 2000; Schmidt et al. 1979; Slifkin and Newell 1999). Where is this signal-dependent noise in force output generated; is it in the planning or execution of motor commands? The variability in continuous isometric force production is thought to arise from the statistical variability and synchrony in the discharge of motoneurons supplying the muscle (Laidlaw et al. 2000; Semmler and Nordstrom 1998; Semmler et al. 2001; Yao et al. 2000). Any process contributing signal-dependent noise (SDN) at any stage in the evolution of a motor command will affect the outflow of action potentials to the muscle. However it is not known whether the contractile mechanism of muscle itself contributes to the increase in variability with the strength of the motor command. It could be that SDN is part of the peripheral machinery of the skeletal muscles that contributes mechanical noise proportional to the amount of activation. If some component of SDN results from physiological characteristics of skeletal muscle, then this will be added on to any noise in the motor command.

The purpose of this study was to localize the sources of the SDN in the neuromuscular system. To do this, we first compared the variability in force production during voluntary isometric contractions and contractions elicited by neuromuscular electrical stimulation (NMES). NMES is a methodology that uses a series of electrical pulses to generate a muscle contraction (Baker et al. 2000). The contraction is elicited indirectly through stimulation of the motor axons with the electrodes usually located over the muscle motor point. The contractions studied were extensions of the distal phalanx of the thumb, a movement that is produced by the action of only one muscle, the extensor pollicis longus (EPL). By using NMES, we could experimentally estimate the noise generated by the peripheral machinery and if that noise had signal-dependent features. Since many of the variables of interest were not directly available for testing and manipulation in human experiments, we used a model of the motor-unit pool to determine how variables such as of order of recruitment, rate coding, and the statistics of motoneuronal firing may contribute to SDN. The main hypothesis was to test whether the SD of force output tended to scale either isometrically [SD(force)  Mean(force)^1.0] or allometrically according to the square root of average force output [SD(force)  Mean(force)^0.5]. Isometric, or linear, scaling implies that the relationship is characterized by a constant CV (CV = SD/mean). A relationship between SD and mean force output that is characterized by a constant Fano factor (Fano = variance/mean) will scale allometrically as a square root (e.g., a Poisson process). We found that the variability of force production tends to scale linearly according to mean force output, and that this scaling is conferred by the physiological properties of the motor-unit pool.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Human experiments

Five healthy young adults participated in the study (3 men and 2 woman, age range 26-37 years). An additional two subjects were tested with similar results, but because they were not tested in all conditions, we have excluded them from this report. All participants indicated that they were right-hand dominant and gave informed consent to the experimental procedures. The experiments were approved by the Committee for Ethics in Human Experimentation at University College London.

EXPERIMENTAL SET-UP. The participants sat with their right arm resting on a table. The forearm was positioned midway between pronation and supination and was stabilized with a vacuum splint. The ulnar surface of the hand rested on the table in a relaxed posture and the hand was fixed in place. The thumb rested on the specially shaped palmer support with the distal phalanx extending beyond the edge of the support and the interphalangeal joint was slightly flexed by 10-20°. The metacarpal and proximal phalanx of the thumb were secured to a foam-backed aluminum splint that was individually fit to each subject and extended onto the forearm along the radius. The immobilization of the hand and forearm allowed isolated movement of the distal phalanx of the thumb by the actions of the flexor pollicis longus (FPL) and EPL muscles. Although we did not directly test for co-contraction by measuring electromyograph (EMG) in FPL, stabilizing the forearm, hand, and thumb minimized this potential problem.

ELECTRICAL STIMULATION. An Odstock 4-channel neuromuscular stimulator (NMES, Department of Medical Physics and Biomedical Engineering, Salisbury District Hospital, Salisbury, U.K.) was modified for PC control. The modification consisted of replacing the analog potentiometers controlling pulse amplitude with digital ones (DS1267-50; Dallas Semiconductor) and interfacing this to the PC via the standard parallel port. The connections from the computer were optically coupled (74OL6000; Fairchild Semiconductors) to maintain subject isolation. Modulation of force output with NMES can be achieved by changes in pulse duration and/or amplitude, which controls recruitment, and pulse frequency that controls the firing rate of motor units. We used a fixed pulse width of 300 µs, a fixed pulse frequency of 25-30 pulses per s (pps), and varied the stimulus amplitude to control the strength of the contraction. At these stimulation rates, the NMES should produce near tetanic contractions of motor units. To circumvent the well-known problem of fatigue during NMES, we used a 1:3 duty cycle, i.e., 7 s on followed by 21 s off.

PROCEDURE. Prior to the start of the experiment proper, the subjects were asked to slowly increase their effort to a maximum voluntary contraction against the force transducer and hold this maximum for 3 s while receiving verbal encouragement. The average peak force during the last 2 s of the contraction was calculated over the three trials. This value was loosely considered the maximal voluntary contraction (MVC), and all forces are reported with respect to this value.

Model of the motor-unit pool

The modeling component of this study relied heavily on a model of the motor-unit pool by Fuglevand and colleagues that has been described in detail (Fuglevand et al. 1993; Yao et al. 2000). The previous studies using this model considered the output of both force and EMG, whereas we have used only those portions of the model concerned with force production. Where our implementation of the model closely follows the previously published version, we will briefly describe the main assumptions and parameters used. In those areas where our implementation departs from the original model of Fuglevand et al. (1993), a more thorough description is given. The model of the motor-unit pool was implemented in the MATLAB environment. The duration of the simulations was 5 s with a time step of 0.5 ms.

RECRUITMENT AND RATE-CODING. The model consisted of a pool of 120 motoneurons that was excited by an excitatory drive distributed uniformly across the motoneuron pool [E, measured in arbitrary excitation units (eu)]. The recruitment threshold of each motoneuron (RTE) was defined by an excitation value that was compared with E to determine whether a given motoneuron was active. The distribution of recruitment thresholds across the pool was modeled by an exponential relationship resulting in a pool with a relatively greater number of low- than high-threshold motoneurons. In some simulations, Gaussian distributed noise with several different SDs was added to the RTE for each motoneuron of the pool. This was done to simulate the finding in human motor unit studies that, while units are generally recruited from small-to-large twitch units, there is some noise in the overall orderly recruitment pattern. Also recent experimental data in the cat and subsequent theoretical analysis has emphasized that fluctuation of spike threshold likely contributes to experimental observations of variability in motoneuron spike trains and synchrony between motoneurons (Binder and Powers 2001; Powers and Binder 2000).

ISI VARIABILITY. In recent years the role of membrane noise in triggering action potentials in human and cat motoneurons has been actively investigated (Kudina 1999; Matthews 1996; Piotrkiewicz 1999; Powers and Binder 2000). A key finding of this research has been the demonstration that the coefficient of variation of the distribution of interspike intervals (ISIs) increases with the mean of the ISI distribution, whereas it was previously held that the coefficient of variation remained constant with changes in mean ISI (Clamann 1969; Fuglevand et al. 1993). It is stressed that these new findings on the statistical distribution of ISIs are applicable within a special region of motoneuron firing called the sub-primary range (Kudina 1999; Matthews 1996; Person and Kudina 1972; Piotrkiewicz 1999). In the sub-primary range, the ISI histogram has a long tail, or positive skew (Matthews 1996), that is better fit to a Rayleigh as opposed to a Gaussian distribution. Outside of the sub-primary range, i.e., at rates >10 pps, the distribution of ISIs becomes more Gaussian in nature.

Data analysis

In most subjects the removal of visual feedback led to drift of the force output even though they were asked to maintain a constant effort (worst case shown in Fig. 1). To eliminate this slow drift component from the force data, trend removal was done on each 4 s data segment using a 2nd order polynomial (Bendat and Piersol 1986). By doing the trend removal, we exclude the waning of the motor command as a source of noise that could contribute to SDN. The final 4 s period of force data during each experimental or simulation trial was used for further analysis. The force was digitally filtered using a 5th order Butterworth filter with a low-pass cutoff of 25 Hz. The cutoff setting was empirically determined by fast Fourier transform (FFT) analysis that showed that >99% of the power in the signal fell between DC 25 Hz for voluntary contractions; thus the bandwidth of the noise signal lies in this frequency range. This setting for the low-pass filter had the added benefit of removing the periodicity due to the slightly unfused contraction at the stimulus rate of 25-30 pps in the NMES condition. The mean force was calculated from the raw data, regardless of any nonstationary trends. The SD was calculated for each trial and averaged across trials of the same target force.

View larger version (27K): [in this window] [in a new window]   Fig. 1. Force output during voluntary and neuromuscular electrical stimulation (NMES)-induced contractions. Isometric force output for subject 1 is shown during voluntary contractions to 3 different target levels (left) and during NMES contractions to similar levels (right). The 4 s period without visual feedback that was analyzed is indicated by the bar over the top trace. The 2nd and 4th columns illustrate this 4 s period of data following filtering and correction for nonstationarity of the force output. SD was calculated from these data and compared with the mean force over the same 4 s period in the data shown in columns 1 and 3. In the voluntary condition, the variability is greater for larger forces, but this relationship is not evident in the NMES condition.

Two features of the relationship between mean force and force variability were analyzed: scaling and magnitude. The main hypothesis was to test whether the SD of force output tended to scale either isometrically (<force>^1.0) or allometrically according to the square root of average force output (<force>^0.5). Isometric, or linear, scaling implies that the relationship is characterized by a constant coefficient of variation (CV = SD/mean). A relationship between SD and mean force output that is characterized by a constant Fano factor (Fano = variance/mean) will scale allometrically as a square root (e.g., a Poisson process). The scaling factor was determined by regression analysis. All regression analysis, both simple and complex models, was done in the MATLAB computing environment. The regression lines were fit by the least squares method. The model of SDN that was fit was
which was logarithmically transformed prior to the regression analysis to
where the slope b is the scaling factor between the mean and SD of force output. The magnitude of the variability was calculated using the regression coefficients and solving for sd with mean = 100% MVC.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

The experimental data from the five subjects will be initially presented to demonstrate the presence of SDN during voluntary isometric contractions and its absence during the NMES condition. All force measurements are normalized to MVC that ranged between 9.7 and 15.6 N, with a mean (±SD) of 12.7 ± 2.2 N. The scaling of the relationship between mean force and force variability (i.e., SD) will be compared with the theoretical relationship proposed by the TOPS model, that is linear scaling. This is equivalent to a log-log relationship between the SD and mean with a slope of 1.0. Following this, simulation results will be presented that highlight the organizational features of the motor-unit pool that are necessary for the linear scaling of SDN.

Experimental data

VOLUNTARY AND NMES-INDUCED CONTRACTIONS. Voluntary isometric contraction of the EPL muscle at increasing levels of mean force resulted in proportional increases in the variability of the force during the steady-state period. This is illustrated for a representative subject in Fig. 1. In the first column, the raw data are illustrated during the voluntary condition at three contraction levels. The 4 s period without visual feedback, indicated by the bar, was cut out of the raw data for further analysis. The second column illustrates these 4 s of data following filtering and removal of nonstationary trends (see METHODS). The data illustrate that the force variability increased as mean levels of force increased.

View larger version (14K): [in this window] [in a new window]   Fig. 2. The scaling of force variability with respect to mean force output. A: data from subject 1 is shown. The SD of the force output is positively correlated to the mean force only in the voluntary condition. B: same data as above, plotted on a log-log scale with regression lines to determine whether the signal-dependent noise (SDN) exhibits linear scaling, i.e., slope = 1.0. Slope in the voluntary condition is 0.99 and in the NMES condition is 0.10, but the latter was not significant.

MIXED CONTRACTIONS. On the whole, the NMES condition generated a fixed level of force variability across a wide range of contraction levels. Thus combined voluntary and stimulation-induced contractions should result in a regression line with the same slope as the voluntary condition offset by the constant variability associated with the NMES condition. Asking subjects to produce a constant target force on top of different levels of stimulation-induced contraction tested this hypothesis. The results averaged across all subjects are presented in Fig. 3. In Fig. 3A, the mean total force at each of six different stimulus levels is shown to be relatively constant (). The voluntary contribution to the total force was estimated by subtracting the mean of the NMES-induced force at each stimulus level (Fig. 3A, ). The predicted voluntary force decreased as the stimulus levels increased and this pattern was mirrored in the force variability as the stimulus levels increased (Fig. 3B).

View larger version (19K): [in this window] [in a new window]   Fig. 3. Group means for the experimental data. A: target force in the mixed contraction condition was approximately 70% maximal voluntary contraction (MVC) and remained relatively constant across the 6 stimulus levels. Voluntary contribution to the target force decreased as the stimulus levels increased. B: SD of the force output in the mixed condition decreased with increasing stimulus levels. C: average data from the 5 subjects (±SD) for the 3 experimental conditions. Data for the mixed condition plots the predicted voluntary force at each stimulus level from A against the SD of force at each stimulus level from B. The slopes (±confidence interval with  = 0.05) followed by r2 and P values for the log-log regressions were: voluntary, 1.07 ± 0.55, 0.88, 0.006; mix, 0.57 ± 0.67, 0.59, 0.08; NMES, 0.29 ± 0.23, 0.76, 0.02.

Model of motor-unit physiology

To determine which aspects of voluntary force production were responsible for generating the linear scaling of SDN, a model of a motor-unit pool was used (Fuglevand et al. 1993). With the model, we tested how the range of motor-unit twitch amplitudes, motoneuron recruitment thresholds, and motoneuron firing variability contributed to the scaling of SDN in the simulated force output.

STRUCTURE OF THE MOTOR-UNIT POOL. Initially it was necessary to determine if the default model produced SDN in the simulated force output. This is illustrated in Fig. 4, A-F. Figure 4, A, C, and E, illustrates the results of the simulations of a single motor unit receiving a stochastic spike train input. The results of the simulations of the pool of 120 motor-units are illustrated in Fig. 4, B, D, and F. These data illustrate the difference between allometric scaling according to the square root and linear scaling. The top traces show the initial 3 s of simulated force output at different levels of excitatory drive. Note that, due to the stochastic firing of the motoneuron, the force output of the single motor unit never reaches tetanus. For one of the single unit simulations, the minimal rhythmic firing rate was set to 5 pps to illustrate the high variability resulting from unfused twitches (Fig. 4A). As illustrated in Fig. 4C, this results in an initial period of decreasing force variability as the twitches become partially fused. However, this is over a limited range following which the variability starts to increase in proportion to the mean force output. The scaling of variability with mean force for the single unit is marginally better fit by a square root function (dashed line) compared with a linear fit (Fig. 4C, solid line). This demonstrates the utility of using the log-log regression to differentiate between linear and square root scaling. In the log-log plot the regression line has a slope of 0.47, which is close to the theoretical value of 0.5 for a square root function (Fig. 4E). Square root scaling between the mean and SD is a hallmark of a process with a constant Fano factor, e.g., shot noise process (Poisson process played through a linear twitch filter; Cox and Miller 1965). It turns out that the scaling of any single unit in the pool follows a square root function. However, the scaling of the output of the whole motor-unit pool was a better fit to a linear model (Fig. 4D, solid line) compared with the square root function (dashed line). The regression analysis on the logarithmically transformed data resulted in a slope of 0.88 (Fig. 4F), which was within the confidence intervals for the mean slope in our experimental data (Fig. 3; Table 1). These results demonstrate that the force output of the whole motor-unit pool displays linearly scaled SDN comparable to the experimental data and thus may be used to explore which aspects of motor systems physiology give rise to SDN.

View larger version (33K): [in this window] [in a new window]   Fig. 4. Force output variability from the model. A: simulations are illustrated for 7 different levels of excitatory drive in a high-threshold unit with a large, fast twitch. Only the first 3 s of the 5 s simulation are shown to allow visual resolution of the twitches. B: simulated force output from the default pool of 120 motor units at 8 different levels of excitatory drive. C: mean force and SD of force for a single simulated motor unit. Both linear (solid line) and square root functions (dashed line) were fit to the data. Root mean square (RMS) error was 50.6842 for the linear fit and 50.6573 for the square root fit. D: mean force and variability for the simulated motor unit pool with default parameters. RMS error was 60.4364 for the linear fit and  for the square root fit. E: log-log regression analysis for the single unit data clearly reveals that force variability does not scale linearly, because the slope, 0.47, is nearer the value of 0.5 for square root scaling. F: log-log regression analysis for the whole motor unit pool exhibits linear scaling with a slope of 0.88. This is within the confidence intervals of the experimental data and close to the theoretical value of 1.0 for linearly scaled SDN. G: decreases in the number of motor units had no effect on the slope of the log-log regression for the whole motor unit pool. H: range of twitch forces in the motor unit pool was a significant factor affecting the slope of the log-log regression. If the range was too small, the slope fell sharply to values close to 0.5. I: range of recruitment thresholds across the motor unit pool was a significant factor affecting the slope of the log-log regression. If all the motor units were recruited in a narrow range of thresholds, the slope fell to values close to 0.5.

RECRUITMENT ORDER. Simulations were done to explore the effect of recruitment order on the scaling of SDN during voluntary contractions. The three recruitment schemes tested are illustrated in Fig. 5A: orderly (), reversed (), and random recruitment (). The distribution of threshold and twitch amplitudes was nonhomogenous in each of the three conditions with the majority of motor-units having smaller twitch amplitudes. The order of recruitment had a notable effect on the relationship between excitatory drive and the mean force output, as illustrated in Fig. 5B. In the random condition the force output rapidly reached maximal output within a narrow range of the excitatory drive. The excitation/force relationship was less steep in the case of reversed recruitment but it was still steeper than in the orderly recruitment scheme.

View larger version (39K): [in this window] [in a new window]   Fig. 5. Recruitment effects on the pattern of force output variability. A: three recruitment schemes (orderly, ; reversed, ; random, ) are illustrated by plotting the peak twitch force for each motor unit against its recruitment threshold. In each case, there are more small twitch motor units than large. B: input/output relationship for the motor unit pool with the 3 recruitment schemes is illustrated. The type of recruitment scheme has a strong effect on the shape of this relationship. Several lines are seen from repeated simulations with random recruitment as this scheme gave slightly different results C: log-log regression analysis is shown in which the slope of the orderly recruitment is the only one that approximates the empirical data, i.e., a slope of 0.88. D: force output variability across the range of force output is shown as the coefficient of variation. This facilitates comparison to previous empirical data and highlights the nonlinear changes in coefficient of variation (CV) of force at low force levels.

NMES SIMULATIONS. To confirm the experimental NMES results, we simulated the NMES condition with the default model. In these simulations, the 120 motor units of the pool were recruited randomly with respect to twitch size and fired with a fixed rate of 30 pps. The simulated force output was low-pass filtered in the same manner as the experimental data to exclude the periodicities in the force output due to partially unfused contraction at the stimulus rate, yet retaining the bandwidth of the noise signal determined from the voluntary condition. The magnitude of the force SD at 100% MVC was 0.02 (%MVC), compared with 0.40 ± 0.17 for the experimental results (Table 1). Thus the NMES simulations result in a near complete lack of noise.

ISI VARIABILITY. The previous simulations have demonstrated that many of the functional features of the organization of motor-unit pools bias the force output to produce linearly scaled SDN. However, none of these simulations have examined the role of noise in the motor command to the motoneurons. The presence of noise in the excitatory drive to the motoneurons will affect the distribution of the ISIs (Matthews 1996). The statistical distribution of the motoneuron spike trains may be an important feature in determining the scaling and magnitude of SDN in the motor output.

View larger version (21K): [in this window] [in a new window]   Fig. 6. The effect of interspike interval (ISI) variability on the pattern of force output variability. A: different types of ISI variability simulated are shown. In the default condition, the variability of the ISI distributions had a constant CV (CV = 0.2). Simulations were also performed with ISI variability matching that of the human sub-primary range data, and with 4 different levels of constant SDs, as indicated on the figure. B: log-log regression analysis revealed little effect of ISI variability on the scaling of signal-dependent noise, i.e., all have similar slopes. However, the absolute magnitude of the force output variability is closely correlated with ISI variability. Data for the simulations incorporating the sub-primary range have not been added to the graph, because they could not be distinguished from the constant CV condition. The slope of the line for the sub-primary range simulations was 0.85. All regression lines had r2  0.90 and P < 0.001.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

The main findings of the study were that voluntary isometric contractions are characterized by the presence of linearly scaled SDN over a wide range of force output. This noise is not a result of peripheral neuromuscular noise because it was not present in the NMES condition. Instead the linear scaling of the SDN is a characteristic of the process of graded voluntary contractions, whether those occur in isolation or on top of an NMES induced contraction. This feature of SDN in isometric contractions was captured in a model of the motor-unit pool where it was found that the scaling of SDN depended on the range of motor-unit forces, the distribution of recruitment thresholds and orderly recruitment. The magnitude of the noise was closely correlated with the variability in the discharge of motoneurons, and by inference fluctuations in premotoneuronal drive, but this variability in and of itself did not affect the linear scaling relationship of SDN.

Voluntary contractions and force variability

The relationship between motor-output variability and the magnitude of the force output has been previously investigated for both discrete and continuous force production protocols (Enoka et al. 1999; Laidlaw et al. 2000; Schmidt et al. 1979; Slifkin and Newell 1999). All of these studies found that the variability in force output increased monotonically with an increase in mean force output, similar to what we have reported here. None of these previous studies were particularly concerned with the scaling relationship between variability and mean force as it relates to the optimization of motor planning.

It has been previously shown that the increased variability in force that accompanies aging is likely due to an increased variability in the discharge rates of motoneurons (Laidlaw et al. 2000). These authors also showed that elderly subjects had a higher frequency of double discharges that would also increase the force variability. While we did not examine the effects of double discharges, it was clear from our modeling results that the variability of motoneuron ISIs had a strong effect on the magnitude of the variability at a given force level without affecting the scaling of the SDN (Fig. 6B). Thus our results support their conclusion that an increase in the variability of motoneuron ISIs will give rise to an increased magnitude of force output variability.

These authors also pointed out that isometric contractions of the first dorsal interosseous are characterized by a CV that is higher for lower forces before reaching a plateau of about 2-5% at a contraction level of 10% MVC (Enoka et al. 1999; Laidlaw et al. 2000). There was some evidence for a similar monotonically decreasing CV in our experimental data; however, we did not sample the forces below a target value of 20% MVC and thus cannot comment further. However, the modeling results presented in Fig. 5D can be directly compared with the previous experimental evidence for an initially high CV that decreases to a plateau. This feature of the CV was particularly susceptible to changes in recruitment order. So while in the control condition with orderly recruitment there was some evidence for a decreasing CV, this became particularly evident in the random and reversed recruitment schemes. Thus it could be that noise in recruitment thresholds accompanying aging may in part give rise to the enhancement of CV at low force levels in old subjects (Laidlaw et al. 2000).

Another factor affecting the overall magnitude of force variability at a given mean force level is the synchronization of motoneuron discharges (Datta and Stephens 1990; Semmler and Nordstrom 1998; Semmler et al. 2001; Yao et al. 2000). Our simulations revealed the presence of linearly scaled SDN in the motor output in the absence of synchronization, and therefore it seems that while synchronization will affect the magnitude of the variability (Semmler et al. 2001; Yao et al. 2000), it is likely to do so without changing the scaling of the SDN. It could be postulated that if synchronization of motoneuron discharges within a pool changed with the overall level of excitatory drive, then this too could be an important factor contributing to the scaling of SDN. However it remains to be empirically determined whether synchronization varies as a function of excitatory drive.

NMES condition

The main purpose of this condition was to determine if peripheral neuromuscular noise contributed to the scaling of SDN. The classical descriptions of muscle force output in response to electrical stimulation were concerned primarily with the mean force output in response to a particular stimulus paradigm (e.g., Rack and Westbury 1969). While it was clear from this earlier work that staggered stimulation of groups of motor units could produce a smoother force output compared with synchronous stimulation, force variability at different levels of mean force was not examined. We found that the variability of the force output in all but one case was not related to the mean force output generated by NMES. In the single case that did show a significant relationship, the variability decreased with an increase in mean force output, opposite the direction of the empirically measured SDN. Thus we conclude that the mechanics of muscular contraction, while they may contribute some fixed amount of noise, do not contribute to SDN. Simulations of the NMES condition supported this conclusion.

To appreciate why neither the experimental data nor the model exhibit SDN in the NMES condition, we must highlight the differences in the way the muscle is generating force in the two conditions. We should initially point out that with the stimulus parameters used in our study there was no evidence for reflex evoked involuntary muscle contractions as has been recently described (Collins et al. 2001). Therefore we are confident that the force output in our NMES experimental paradigm is primarily due to direct stimulation of the alpha motor axons. However there are notable differences between the NMES condition and the voluntary condition including the following: 1) all active motor units fire synchronously in response to the stimuli (25-30 pps); 2) increases in mean force output are achieved by recruitment alone, i.e., there is no rate coding; and 3) recruitment order during NMES is more random compared with the voluntary condition. Most importantly, at the stimulation rates used, it is likely that most motor units will be contracting near tetanically (Macefield et al. 1996; Nathan and Tavi 1990; Thomas et al. 1991). Thus our NMES paradigm is sensitive only to mechanical noise generated during tetanic contractions summed across active motor units. This is in contrast to the voluntary condition where even at high levels of excitation the motoneurons fire in a stochastic manner giving rise to variability in the force output (e.g., Fig. 4A). The synchronous discharge of the motor units in the NMES condition would tend to increase the force variability, but because the contractions are tetanic, the effect of synchrony is minimized. In addition, the effects of synchrony are summation of twitches at the stimulus rate which is >25 pps and therefore outside the bandwidth of the SDN signal determined from the voluntary condition. The random recruitment order arises due to the relationship between percutaneous stimulation, the distance of the axons from the current source, and the range of alpha motor axon diameters at the motor point (see Singh et al. 2000). While the recruitment of axons to electrical stimulation is biased so that larger axons will be excited at lower stimulus amplitudes, this is true primarily in the condition where the nerve is in intimate contact with the electrodes. Percutaneous stimulation, as in this study, will recruit the alpha motor axons in a more random order because the most important factor determining excitability will be distance from the electrodes. Furthermore, even if distance from the electrodes was not a factor, it is not clear that that recruitment due to stimulation would proceed from large twitch to small twitch motor units. This is because the positive correlation between axon diameters, as measured by conduction velocity, and twitch force output in the cat hind limb is not evident in humans (Bigland-Ritchie et al. 1998). Taking all these factors into consideration, the result is that the NMES condition is likely to recruit motor units randomly with respect to their twitch tension.

Twitches, threshold, and recruitment order

The simulation results showed that the pattern of increased force variability with increases in mean force that characterizes SDN depended on the large range of twitch forces, the distribution of recruitment thresholds, and the orderly recruitment of motor-units in a pool, i.e., the underlying physiology of the motor-unit pool. A motor-unit pool composed of units with the same twitch amplitude did not generate SDN, nor did a pool in which all motor-units had the same recruitment threshold (Fig. 4, G and H). Even given a motor-unit pool with a broad distribution of twitch amplitudes and recruitment thresholds, it was necessary to recruit these in an orderly fashion to reproduce the pattern of SDN seen in the experimental data.

In the default model, the twitch range was 1-100 and the recruitment threshold range was 30, meaning that the last motoneuron was recruited at about 60% MVC. Can these values be justified based on experimentally determined ranges in human studies? A recent review by Chan et al. (2001) has tabulated the contractile properties of human motor units from a number of upper and lower limb muscles. Important considerations in interpreting the human data are the technique used and the number of motor units sampled. The first dorsal interosseous (FDI) muscle has been widely studied using the techniques of spike-triggered averaging (STA) and intramuscular stimulation (IMS). Although each of these techniques has some drawbacks in terms of sampling, there appears to be no systematic difference in the range of peak twitch force reported with the two techniques in studies with >150 motor units. The reported range of peak twitch forces is from <1 to >100 mN, thus our default range is appropriate for the FDI muscle. It is less clear that a similar range of peak twitch forces exists in other human muscles. In the thenar muscles, the range appears more compressed, but this may simply be a function of low numbers of motor units sampled. For the extensor carpi radialis muscle, there is a large discrepancy with a range of 40 reported in one study and a range of >100 reported in another. The model predicts that if the range of twitch tensions is <50, then the scaling of SDN will depart from linear (Fig. 5H). However, this is restricted to the case of a muscle acting in isolation, which is not the case across most articulations. It remains to be empirically determined what the scaling of SDN is across a joint with multiple synergists. The default range of recruitment threshold is less problematic. For hand muscles, recruitment occurs over the first 50% MVC, but this range is extended to 85% MVC for limb muscles (reviewed by Enoka and Fuglevand 2001). The model results suggest that as long as recruitment occurs over at least the first 30% MVC, the scaling of SDN will tend to be linear (Fig. 5I).

While recruitment studies on human motor units have emphasized the importance of orderly recruitment during isometric contractions, it is also clear that there is some noise in the overall orderly recruitment pattern (Desmedt and Godaux 1977; Milner-Brown et al. 1973). This noise is mainly apparent in the lower threshold units where the differences between recruitment thresholds are small. Additionally, recent experimental data in the cat and subsequent theoretical analysis has emphasized that fluctuation of spike threshold likely contributes to experimental observations of variability in motoneuron spike trains and low levels of synchrony between motoneurons commonly reported in human motor unit studies (Binder and Powers 2001; Powers and Binder 2000). While this latter work has been primarily concerned with variation in threshold between spikes, it indicates that recruitment threshold for a motoneuron is not static but likely exhibits some variability.

One of the key factors determining the force output of a motor unit, and therefore the range of twitch forces in a motor-unit pool, is the innervation number. It has been empirically determined and theoretically estimated that the innervation numbers and resulting motor unit force outputs in human muscle are not homogenously distributed (Enoka and Fuglevand 2001; Garnett et al. 1979; Thomas et al. 1990). Instead the distribution of motor units according to force is skewed with a majority of motor units producing small forces. Theoretical studies have concluded that the optimal distribution of motor-unit forces depends on the probability distribution function (pdf) of the forces generated by the muscle. If this pdf is monotonically decreasing, then an optimal distribution of twitch/tetanic force outputs for a fixed number of motor-units, N, and a constant value of MVC will also be monotonically decreasing. That is, if the muscle produces small force outputs more frequently than large force outputs, then to produce the finest resolution of force with respect to usage, it will be optimal to have a greater number of small than large twitch motor-units (Senn et al. 1997; Tax and Denier van der Gon 1991). In addition, such a distribution of motor unit forces is optimal from an information theoretical point of view in a motor-unit pool that regulates force by pure recruitment modulation (Senn et al. 1997).

Optimization of motor output

There have been many post hoc explanations for the benefit of orderly recruitment as well as much experimental work designed to determine the physiological explanation for orderly recruitment first enunciated by Henneman (Binder and Mendell 1990; Henneman et al. 1965). As argued by Senn et al. (1997), orderly recruitment according to the size of motor unit force output minimizes the error between the input, modeled as required force, and the force output. We have shown that the sequela of this pattern of recruitment is SDN in the force output.

We have shown a particular scaling of SDN, linear scaling over a wide range of force output. There are two reasons for emphasizing this particular scaling of SDN. First, the experimental data support this pattern of SDN; the mean slope from the log-log regression analysis was 1.05 ± 0.48 (Table 1). Second, the value of the log-log slope has significant effects on the type of control strategy used in the TOPS model (Harris and Wolpert 1998). If the slope were 0.5, i.e., noise with a square root scaling, then the optimal control strategy for minimizing endpoint errors would be bang-bang control with the motor commands taking either values of zero or its maximum value (unpublished observations). Conversely, with a slope of 1.0, i.e., SDN with a constant CV, the optimal strategy is one generating a continuous, i.e., smoothly varying motor command covering the entire range of possible values, and the output trajectories match experimental observations (Harris and Wolpert 1998) The effect of nonlinear decreases in the CV in the low force ranges, shown in Fig. 5D and highlighted in some experimental studies (Enoka et al. 1999; Laidlaw et al. 2000), on the TOPS model have yet to be evaluated.

Thus it would appear that the cost of optimization of the force output at the level of a single muscle is linearly scaled SDN in the force output. The CNS accounts for the SDN in planning movements so that the optimal trajectory is constrained by this feature of the motor system.