## Abstract

Computational models of motor-unit populations are the objective implementations of the hypothesized mechanisms by which neural and muscle properties give rise to electromyograms (EMGs) and force. However, the variability/uncertainty of the parameters used in these models—and how they affect predictions—confounds assessing these hypothesized mechanisms. We perform a large-scale computational sensitivity analysis on the state-of-the-art computational model of surface EMG, force, and force variability by combining a comprehensive review of published experimental data with Monte Carlo simulations. To exhaustively explore model performance and robustness, we ran numerous iterative simulations each using a random set of values for nine commonly measured motor neuron and muscle parameters. Parameter values were sampled across their reported experimental ranges. Convergence after 439 simulations found that only 3 simulations met our two fitness criteria: approximating the well-established experimental relations for the scaling of EMG amplitude and force variability with mean force. An additional 424 simulations preferentially sampling the neighborhood of those 3 valid simulations converged to reveal 65 additional sets of parameter values for which the model predictions approximate the experimentally known relations. We find the model is not sensitive to muscle properties but very sensitive to several motor neuron properties—especially peak discharge rates and recruitment ranges. Therefore to advance our understanding of EMG and muscle force, it is critical to evaluate the hypothesized neural mechanisms as implemented in today's state-of-the-art models of motor unit function. We discuss experimental and analytical avenues to do so as well as new features that may be added in future implementations of motor-unit models to improve their experimental validity.

## INTRODUCTION

Computational models of motor-unit populations are, by definition, the objective implementations of the mechanisms by which neural and muscle properties are hypothesized to give rise to electromyograms (EMGs) and force. Unlike regression-based inferences, these so-called forward or predictive models (Baker and Lemon 1998; Fuglevand et al. 1993) are causal in that they use mathematical equations that simulate the flow of events from neural command to electromagnetic signals recorded at the electrodes or forces produced at the tendons. The degree to which a model is able to replicate well-established experimental relations is a test of how well those equations reflect the underlying mechanisms and a test of our hypotheses of neuromuscular function. Thus the rigorous study of how to make these models as experimentally valid as possible is critical to improve both our understanding of neuromuscular function, and our ability to use EMGs to confidently infer command signals and forces in the neuromuscular system.

Briefly, the EMG signal is the sum of motor-unit potentials detected by electrodes and is the most common experimental means to estimate the neural drive to a muscle (Lacquaniti and Soechting 1982; Schieber 1995; Thoroughman and Shadmehr 1999; Valero-Cuevas et al. 1998) and muscle force (Inman et al. 1952; Milner-Brown and Stein 1975). Because many factors influence the EMG signal (Basmajian and De Luca 1985), most researchers agree that it has important limitations (Day and Hulliger 2001; Farina et al. 2004a). Due to these limitations and the current absence of other approaches to examine interactions of entire motor-unit populations, computational models are fast becoming the dominant means to establish the mechanisms relating EMG to muscle force and neural drive (Baker and Lemon 1998; Fuglevand et al. 1993; Stegeman et al. 2000). Structure-based EMG models (reviewed in McGill 2004; Stegeman et al. 2000) use motor units as their fundamental building blocks and promise to indicate how to best use EMG to estimate muscle force and the neural drive to muscle in a behaving animal.

Critically assessing these computational models of motor units is relevant and urgent because in the absence of other computational means, these models are routinely used to study and debate the effects of motor neuron and muscle parameter values on muscle function. Some current debates center on peak discharge rates of motor units (Connelly et al. 1999; Enoka and Fuglevand 2001; Johns and Fuglevand 2004), recruitment ranges (De Luca et al. 1982; Moritz et al. 2005), number of motor units (Hamilton et al. 2004; Heckman and Enoka 2004), and range in innervation numbers (Enoka and Fuglevand 2001). However, the variability and uncertainty of the parameters used in these models confounds assessing how realistic they are. Given the excessive computational cost needed to model the behavior of an entire motor-unit pool, typically requiring millions of calculations, simulations have previously been restricted to examining the influence of independent changes in single parameters on single outputs (Jones et al. 2002; Keenan et al. 2005; Zhou and Rymer 2004). Validation efforts to date have not addressed multivariate interactions, the potential confound of an inappropriate model structure (Valero-Cuevas 2005), intra- and inter-subject variability in motor-unit properties (Feinstein et al. 1955), nor the fact that many properties across an entire population of motor units are largely unknown (Enoka and Fuglevand 2001).

We use multivariate Monte Carlo simulations to exhaustively and simultaneously explore the sensitivity of model predictions to uncertainty and variability in nine of the most commonly investigated motor-unit parameters. Motor-unit models should, at the very least replicate known relationships among force, force variability, and EMG amplitude and be robust to parameter variability across experimentally identified ranges for motor neuron and muscle properties. It is well documented that both EMG amplitude and force variability scale with mean force level (Bigland-Ritchie 1981; Enoka et al. 2003; Jones et al. 2002; Milner-Brown and Stein 1975). This study examined the ability of the most popular motor-unit model to replicate experimental EMG/force and force/force-variability relations in the presence of variability and uncertainty in parameter values. A preliminary account of part of these findings has been published in abstract form (Keenan and Valero-Cuevas 2006).

## METHODS

We used Monte Carlo simulations to test the ability of the most popular motor-unit model available (Fuglevand et al. 1993) to robustly predict well-established relations among EMG, force, and force variability when changing the values for motor neuron and muscle parameter values across experimentally identified ranges. Our computer simulations are based on the forward model of recruitment of a population of motor units and production of muscle force (Fuglevand et al. 1993) with the addition of an EMG model (Farina and Merletti 2001) that simulated the surface EMG with a planar volume conductor with muscle, fat, and skin tissues. Motor-unit potentials and the EMG signal were computed at 4,096 sample/s and muscle force at 500 sample/s. Monte Carlo methods, as their name implies, use exhaustive repeated simulations where the specific value of each free model parameter is assigned by random chance from a preassigned range, and the output is computed (Santos and Valero-Cuevas 2006; Valero-Cuevas et al. 2003). Each iteration consists of a simulation that uses one given set of values for the free parameters (see details in the following text). On convergence the results are in the form of the range and variability in the output in response to changes in the parameter values—which is a rigorous evaluation of the robustness and realism of the structure of the model in relation to the variability/uncertainty in parameter values. “Convergence” in Monte Carlo simulations has the specific meaning that the number of simulations performed has sufficed to produce stable estimates of the statistical parameters of the output distributions. Importantly, convergence certifies that further iterations are no longer useful because they would not alter the results. Thus converged Monte Carlo simulations are a computationally rigorous means to thoroughly explore large parameter spaces (Santos and Valero-Cuevas 2006; Valero-Cuevas et al. 2003). The model and Monte Carlo iterations are implemented in Matlab version 7.0 (The Mathworks, Natick, MA).

As described in the extensive literature using this model (Jones et al. 2002; Keenan et al. 2006a; Yao et al. 2000) each iteration of our Monte Carlo simulations comprised six main sequential steps: *1*) determining the recruitment and discharge times of a population of motor neurons in response to the level of excitatory drive (we used 6 levels for each Monte Carlo simulation); *2*) generating motor-unit potentials from a given number, location, and conduction velocities of the muscle fibers for each motor unit; *3*) simulating the surface EMG by summing the trains of motor-unit potentials; *4*) generating twitch forces from the given number of muscle fibers for each motor unit; *5*) implementing the sigmoidal relation between motor-unit force and discharge rate frequency as a function of the instantaneous discharge rate and contraction time of the unit; and *6*) simulating muscle force by summing motor-unit forces. On convergence, the cumulative results from the Monte Carlo simulations provide statistical distributions of average full-wave rectified EMG amplitudes, and the mean and SD of force that result from the reported variability and measurement uncertainty in motor-unit parameters (Table 1).

Because the original equations (Farina et al. 2004b; Fuglevand et al. 1993) and our implementation (Keenan et al. 2005, 2006a,b) of the model have previously been published, the following section will only review key components of the model and focus on our novel Monte Carlo implementation.

### Motor neuron pool model

The distributions of properties across the motor-unit pool were based on the size principle (Henneman 1957), and these included innervation number, twitch force, recruitment threshold, motor-unit territory, and conduction velocity of motor-unit potentials. The total number of motor neurons in the pool varied from 150–500 and the number of muscle fibers from 50,000–250,000 (Table 1).

Activation of the motor-unit pool is modeled as a ramp-and-hold function, with a 1-s ramp increase in excitation to a constant level for 10 s. Excitation ranges from 0 to 1, with maximal excitation defined as the level of activation necessary to bring the last recruited motor neuron to its assigned peak discharge rate. The distribution of recruitment thresholds for the motor neurons is represented as an exponential function with many low-threshold neurons and progressively fewer high-threshold neurons. The range over which recruitment occurs varies from ≤30 to 80% of maximal excitation (Table 1), consistent with the wide range found experimentally (De Luca et al. 1982; Moritz et al. 2005). Each motor unit began discharging at 8 pulse/s (pps) once excitation exceeded the assigned recruitment threshold for that unit, and discharge rates increase linearly with increases in excitation. The interspike interval for motor-unit discharge is modeled as a random process with a Gaussian probability distribution function, with the coefficient of variation [CV = 100 × (SD/mean) %] of the interspike intervals fixed for all motor units at 20%. The influence of the correlated discharge among motor units (synchronization) was not investigated as it has minimal influence on normalized EMG amplitudes (Keenan et al. 2005; Zhou and Rymer 2004).

High-threshold motor neurons have been reported to obtain either higher (Gydikov and Kosarov 1974; Moritz et al. 2005), lower (De Luca et al. 1982), or similar (Erim et al. 1996) peak discharge rates with respect to low-threshold motor neurons. Given that evidence for only one discharge rate strategy has not been clearly established in vivo, peak discharge rates are assigned for the first- and last-recruited motor neuron independently and are drawn from a uniform distribution that varied from 25 to 50 pps (Table 1).

### Muscle force generation model

Motor unit twitch force is modeled as the impulse response of a critically damped, second-order system (Milner-Brown et al. 1973a). A large number of units produce small forces, whereas relatively few units generated large forces. Twitch force is modeled as directly proportional to the number of fibers assigned to each motor neuron (innervation number), and the range in innervation numbers, calculated by dividing the innervation number of the last-recruited motor unit by the first-recruited motor unit, varies from 1 to 100× (Table 1). Variation in the range in innervation numbers, and thus twitch forces in the model, is simulated because experimental reports using glycogen depletion techniques find ranges <2×, whereas estimates of ranges in innervation numbers using twitch and tetanic motor-unit forces, which are highly correlated with innervation number (Bodine et al. 1987; Kanda and Hashizume 1992), can extend ≤100× (reviewed in Enoka and Fuglevand 2001). Innervation numbers increase exponentially from the smallest to the largest motor unit. There is an inverse relation between twitch force and twitch contraction time. The first- and last-recruited unit are assigned a twitch contraction time to peak of 90 and 30 ms, respectively. All motor units follow the established sigmoidal relation between motor-unit force and motor-unit discharge rate (Bigland and Lippold 1954; Fuglevand et al. 1993; Kernell et al. 1983). Total force of the muscle is calculated as the linear summation of all the individual motor-unit forces.

### Surface EMG simulations

The planar volume conductor consists of an isotropic fat (3-mm thick) and skin layer (1-mm thick), and an anisotropic muscle layer. Each tissue layer is homogeneous and the conductivity properties of each layer are similar to those reported by Farina et al. (2004b). The cross-section of the muscle is simulated as circular with a radius proportional to the assigned number of muscle fibers (Table 1). Average fiber length varies from 4 to 16 cm (Table 1), and the center of the innervation zones is located at 50% of the length of the fibers and in the same transverse location within the volume conductor. The end-plate and insertion of each fiber into the tendons varies randomly (uniform distribution) over a range of 5 mm. The electrodes are simulated as circular (4-mm diam) and arranged in bipolar pairs (10-mm interelectrode distance) along the direction of the muscle fibers. The electrodes are located halfway between the innervation zone and the distal tendon and are placed directly over the center of the muscle in the radial direction. The simulation of the intracellular action potential is based on the analytical description of Rosenfalk (1969). Motor units have mean muscle fiber conduction velocities that vary from 3 to 4 m/s and have a spread (SD) that varies from 0.0 to 0.5 m/s (Farina et al. 2000) (Table 1), with the smallest motor units assigned the slowest conduction velocities (Andreassen and Arendt-Nielsen 1987). The motor-unit potential comprises the sum of the action potentials of the muscle fibers belonging to the motor unit and the surface EMG is simulated by summing the trains of motor-unit potentials.

The distribution of the motor units within the muscle (Johnson et al. 1973; Milner-Brown and Stein 1975) is determined by randomly selecting the *x-y* coordinates corresponding to the center of each motor-unit territory within the circular cross-section of the muscle. The fibers of a motor unit are randomly scattered in the motor-unit territory (Stålberg and Antoni 1980) with a density of 20 fibers/mm^{2} (Armstrong et al. 1988) and interdigitated with fibers belonging to many other units. The territories of the largest motor units, therefore are greater than those of the smallest motor units (Bodine et al. 1988). To restrict the distribution of the fibers of a motor unit (Bodine et al. 1988; Stålberg and Antoni 1980), motor-unit territories are modeled as circular (Buchthal et al. 1959). The radius of the motor-unit territory is based on the density of 20 fibers/mm^{2} and the assigned innervation number of the motor unit. However, when a portion of the motor-unit territory is constrained by the muscle boundary, the radius of the territory of the unit is augmented to accommodate the required number of motor-unit fibers within the new territory while maintaining fiber density.

### Monte Carlo simulations

We investigated how changes in five motor neuron and four muscle parameter values (Table 1) affect the predicted relations among EMG, force, and force variability. These nine parameters are specifically chosen because they are intensively investigated in animal and human experiments and also because the influence of these properties on EMG and force production is the focus of the majority of studies using motor-unit models (e.g., Fuglevand et al. 1993; Jones et al. 2002; Keenan et al. 2005). Three muscle properties (mean and variability in motor-unit conduction velocities, and fiber length) were varied to specifically test their influence on the surface EMG. The Monte Carlo approach is demonstrated in Fig. 1 for the generation of the surface EMG with a population of 150 motor units. For each Monte Carlo iteration, we defined a set of parameter values by sampling at random from uniform distributions for each of the nine parameters bounded by the experimental range of values commonly reported for each (Table 1). For each iteration, an activation pattern (motor-unit recruitment and discharge rates) was simulated at six levels of excitation (arrows in Fig. 1) based on the motor neuron property values drawn. For the anatomical factors previously described in methods (Fig. 1), and the motor neuron and muscle properties drawn for each iteration, motor-unit potentials are generated for each motor unit, and an EMG signal is simulated at each excitation level by summing all trains of motor-unit potentials. The same Monte Carlo approach is used to simulate muscle force.

### Data analysis

Force and full-wave rectified EMG are averaged over a 10-s window during constant excitation, and force variability is quantified as the SD of the detrended force. Convergence was declared when there was <2% change in the running mean and coefficient of variation of EMG amplitude and force for the last 20+% of the simulations (Santos and Valero-Cuevas 2006; Valero-Cuevas et al. 2003). Note that each Monte Carlo iteration itself required convergence at every one of the six levels of excitation.

Because EMG is influenced by the distribution of motor-unit territories within the muscle, published simulations generally involve an arbitrary number (e.g., 20) of random motor-unit locations (Farina et al. 2002; Lowery et al. 2003), while holding all other parameters fixed (e.g., discharge patterns, conduction velocities, and number of fibers for each motor unit). To ensure that this variability did not affect our results, EMGs within each Monte Carlo iteration were also computed iteratively with randomly located motor-unit territories until EMG amplitudes also met the 2% convergence criteria. On average, 107 such iterations sufficed.

Two fitness criteria were chosen to indicate a reasonable match between simulated and experimental EMG/force and force/force-variability relations. These relations were determined by regression analysis done in MATLAB. First, force variability is known to scale with mean force (Enoka et al. 1999; Moritz et al. 2005; Slifkin and Newell 2000), and Jones and Wolpert (2002) reported that the mean slope of the log-log regression line between SD and mean force was 1.05. We chose a slope >0.75 and <1.25 to indicate a log-log linear (slope = 1) scaling consistent with signal-dependent noise (Harris and Wolpert 1998). And second, EMG amplitude is known to scale with force (Bigland-Ritchie 1981; Inman et al. 1952; Milner-Brown and Stein 1975), with the EMG/force relation described as either linear or nonlinear—with EMG increasing less rapidly than force at moderate contraction levels. We chose a slope in the regression line between EMG and force of <1.05 as a liberal estimate of a match with experimental observations so as to not over-constrain the model output a priori. To calculate EMG/force relations, EMG and force were normalized with respect to their values at maximal excitation for each Monte Carlo iteration. Convergence was assessed using output distributions of EMG amplitude (Fig. 1) and force, instead of using output distributions of the EMG/force and force/force-variability slopes. This resulted in a more rigorous sampling of the parameter space because more than twice as many Monte Carlo iterations were necessary to satisfy the convergence criteria on the former than on the later output distributions.

On convergence, we identified the iterations that met the experimental fitness criteria and analyzed two aspects of their associated sets of parameter values. First, we established and compared the ranges of values of each model parameter for which the model produces experimentally valid simulations. This explicitly identifies the parameters to which the model is most sensitive. Second, we rank ordered these sets of valid parameter values based on how well they met the experimental fitness criteria. This provides specific sets of parameter values for others to run experimentally valid simulations with different fitness emphases.

## RESULTS

The initial Monte Carlo simulations converged after 439 iterations, but only three simulations approximate the experimental relations for force/force-variability and EMG/force. Figure 2 shows the 439 Monte Carlo iterations plotted against the required experimental relations in force/force-variability and EMG/force. Horizontal and vertical dashed lines indicate boundaries for the fitness criteria representing an approximate match with experimental data (Fig. 2*C*), with quadrant IV containing those iterations that match both relations. In general, the model showed a predominant nonlinear scaling of force variability (Fig. 2*A*, green) and EMG (Fig. 2*B*, red) with mean force across 6 excitation levels, in contradiction to what is observed experimentally. Only 3 of the 439 Monte Carlo iterations resulted in an approximate match with experimental relations (Fig. 2, black). Histograms for the nine motor neuron and muscles properties are depicted in Fig. 3 for the parameter values that lead to an approximate match with experimental findings.

We found a tendency toward a tradeoff between fitness criteria where satisfying one fitness criterion was negatively associated with satisfying the other. Specifically, the slopes of the EMG/force relations were analyzed as a function of the slopes of the force/force-variability relations (Fig. 2*C*) with a first-order regression model, and the two were slightly positively correlated (*r* = 0.38, *P* < 0.05). This suggests that using a set of parameters that replicate experimental EMG-force relations tended not to replicate the force/force-variability relations. In fact, close inspection of the valid parameter sets suggests that satisfying each fitness criterion involved preferentially sampling from different regions of the parameter space. For example, parameter values for peak discharge rate of the last recruited motor neuron that resulted in valid force/force-variability relations were always <35 pps, while fewer values <35 pps resulted in valid EMG/force relations (Fig. 4).

We performed a second set of Monte Carlo simulations, now preferentially sampling from the parameter ranges identified in Fig. 3, to better characterize the subset of the parameter space that might lead to simulations compatible with experimental data. Specifically for this second set of simulations, the probability of selection from the truncated range bounded by observations in Fig. 3 was 0.8, whereas the probability of sampling from the remaining portion of the range was 0.2 (Fig. 3, gray region). Convergence was met after 424 iterations for the second set of simulations (Fig. 5), and 65 iterations met our fitness criteria. Figure 6 depicts force/force-variability (*A*) and EMG/force relations (*B*) for those 65 valid iterations.

Examination of the histograms for the parameter sets (Fig. 7) from the second set of simulations (seen in Figs. 5 and 6) show that the model is not robust to parameter variability and is most sensitive to motor neuron properties. This lack of robustness and greater sensitivity is made evident by the fact that three of the nine histograms cover <50% of the experimentally-identified ranges, and all three are motor neuron properties (Fig. 7). Specifically, the parameter sets that met with our experimental fitness criteria required peak discharge rates for the first- and last-recruited motor neuron to be <33 pulses/s, and recruitment ranges to be >60% maximal excitation. In contrast, each muscle property could vary across its full range and still produce valid EMG/force and force/force-variability relations. In addition, a final Monte Carlo simulation injecting small uniformly distributed perturbations (±5%) to the best parameter set (i.e., that most closely approximated the linear scaling of EMG and force variability with force) led to valid solutions for 132 iterations, after converging with 232 iterations (only 56.9% success rate). That best parameter set (Fig. 5, blue asterisk) was selected from the 863 simulations shown in Figs. 2 and 5. The *inset* in Fig. 5 depicts the fitness for those 232 iterations.

Table 2 presents the 50 parameter sets from all 1,095 (i.e., 439 + 424 + 232) Monte Carlo iterations that best approximate experimental findings, ranked in order of their proximity to the nominal fitness point (1,1), which is the reproduction of the commonly reported linear scaling (slope = 1) of force variability and EMG amplitude with mean force. As the slopes of simulated EMG/force relations were always >1 while satisfying both fitness criteria (Figs. 2*C* and 5), it was not possible to assess experimentally observed nonlinear (slope <1) EMG/force relations.

## DISCUSSION

Most importantly, our work systematically evaluates the mechanisms by which neural and muscle properties are hypothesized to give rise to EMG and force. The equations in the model are, in their essence, the objective implementation of these hypotheses. Finding that the model is most sensitive to neural properties indicates that to advance our understanding of EMG and muscle force, it is critical to evaluate and potentially revise the equations representing the neural mechanisms in this popular model. As discussed in the following text in detail, it is likely critical to re-evaluate how recruitment ranges and peak discharge rates affect EMG and muscle force. In addition, we find that the model, as implemented in today's literature, tends to trade-off between the two well-established experimental relations we used as validation criteria: the EMG/force and force/force-variability relations. Nevertheless, we found parameter sets for which the model approximates both these experimental relations. We begin by discussing limitations and finish by suggesting experimental and computational means to improve our understanding and representation of neural mechanisms in the model, as well as additional features that may be added to the model, to make it as experimentally valid as possible.

Our results should be interpreted with respect to the limitations of the modeling approach, which we nevertheless believe do not affect the validity of our conclusions. All Monte Carlo simulations are limited by assuming parameter independence, resulting in parameter sets and output distributions that may be broader than are likely to exist in reality. However, when there is a lack of experimental evidence describing whether and how the free parameters co-vary (Table 1), it is necessary to assume parameter independence to obtain an unbiased sense of all possible model outputs. Additionally, this first large-scale computational study did not vary all possible parameters in the model because the number of simulations to convergence grows exponentially with the number of free parameters. Rather we focused on those parameters that are most commonly investigated and segregate into intuitively neural versus muscle properties. The other possible parameters were retained as fixed to values most commonly used in the literature (Fuglevand et al. 1993; Keenan et al. 2005). Our results now enable and guide the future exploration of the hypotheses implemented in the model by the systematic expansion of the number of free parameters to determine if a better match with experimental data are possible. As stated in the preceding text, the choice of our fitness criteria was motivated by a desire to approximate the well-established experimental relations for the scaling of EMG amplitude and force variability with mean force. A limitation of our fitness criteria could be that the slope of the force/force-variability relationship on a log-log scale tends to mitigate the impact of force variability at low activation levels (Laidlaw et al. 2000). Future studies could certainly refine or change fitness criteria. Last, alternative motor-unit model structures have been reported (e.g., simulating the synaptic inputs onto the motor neuron pool) (Baker and Lemon 1998; Lowery and Erim 2005). However, those models likely share similar hypothesized mechanisms and may therefore be susceptible to similar limitations as those presented here. Nonetheless future work can perform similarly thorough exploration of their parameter spaces.

The dominant sensitivity of the model to motor neuron parameter values indicates the equations representing neural function may require refinement. This conclusion is reached using two lines of reasoning. First, when the model is most sensitive to a few parameters (i.e., defined here as spanning a narrow range of values for a group of valid simulations, Fig. 7), it means that the equations using those parameters are most critical to producing those valid solutions.^{1} In this case, we find that the equations using parameters associated with neural properties of the motor-unit populations are most influential to model fitness. Second, the model, as implemented in the literature, needs refinement because it does not produce simulations whose fitness is comfortably within the experimentally established ranges for *any* combination of parameter values [i.e., the distribution of valid simulations (Fig. 5) does not include the nominal fitness point (1,1)]. The reason this can be concluded is that a converged Monte Carlo simulation indicates that further sampling of the parameter space will not produce qualitatively different results (Santos and Valero-Cuevas 2006; Valero-Cuevas 2005); thus the nominal fitness point (1,1) will not be the mode of the distribution of valid simulations even if the model is run for an infinite number of iterations. Furthermore, while satisfying both fitness criteria, EMG amplitude was never less than force at intermediate levels of activation (Fig. 6*B*) as is sometimes observed experimentally (Bigland-Ritchie 1981; Lawrence and De Luca 1983). Therefore the model warrants improvement, and this improvement may be brought about at this point by refining the equations in the model representing the neural mechanisms of motor-unit function.

The narrow ranges of valid motor neuron parameter values (Fig. 7), which resulted in part by the trade-off between EMG/force and force/force-variability relations (Fig. 2*C*), are difficult to reconcile with current debates in the literature over those parameters. Specifically, peak discharge rates for the first and last recruited motor unit had to be <33 pps, and the range over which recruitment occurred had to be >65% maximal excitation to reproduce experimentally valid solutions. In contrast, experimental studies report peak discharge rates ≤50 pps (Bellemare et al. 1983; Kamen et al. 1995) and recruitment complete by <50% maximal force (De Luca et al. 1982; Milner-Brown et al. 1973b). Also, the need to preferentially sample different parameter values (e.g., peak discharge rates—Fig. 4) to independently satisfy each fitness criterion must necessarily be an inherent artifact of the model's structure. This disparity between parameter values necessary for the model to approximate realistic simulations versus parameter values known to exist experimentally likely limits the utility this model structure to study muscle function when assuming peak discharge rates ≥35 pps (Jones et al. 2002; Zhou and Rymer 2004) and recruitment ranges <50% maximal excitation (Fuglevand et al. 1993; Keenan et al. 2005). However, having identified these narrow parameter ranges is particularly useful because they identify which aspects of the model structure (i.e., equations hypothesized to represent neural mechanisms of motor-unit function) would benefit most from systematic evaluation. For example, peak discharge rates define maximum excitation in the model (see methods), determine motor-unit recruitment and discharge rate at each level of activation, and establish when discharge rates saturate with increases in activation of the motor-unit pool.

Importantly, our exhaustive search of the nine-dimensional parameter space using 1,095 Monte Carlo iterations did identify sets of parameter values that resulted in an acceptable match with the experimental observations describing EMG/force and force/force-variability relations (Figs. 2*C* and 5, quadrant IV, Table 2). Although experiments and motor-unit models have consistently characterized relations between both force and force variability (Jones et al. 2002; Laidlaw et al. 2000; Moritz et al. 2005) and EMG (Bigland-Ritchie 1981; Inman et al. 1952; Milner-Brown and Stein 1975), to our knowledge no previous study has evaluated motor-unit models by comparing their predictions against both of these well-established relations. For example, previous univariate approaches addressed which factors influence force/force-variability relations, including: discharge rate variability (Moritz et al. 2005), range in twitch amplitudes (Jones et al. 2002), and number of motor units (Hamilton et al. 2004). In addition, previous univariate approaches addressed which factors influence EMG/force relations, including motor-unit synchronization (Yao et al. 2000), electrical and mechanical properties at the individual motor-unit level (Zhou and Rymer 2004), amplitude cancellation (Day and Hulliger 2001; Keenan et al. 2005), and recruitment range and peak discharge rates (Fuglevand et al. 1993; Milner-Brown and Stein 1975). By performing an exhaustive multivariate analysis that compared modeled and experimental EMG/force and force/force-variability relations, we obtain a true indication of the capabilities of this model structure as implemented in today's literature. A list of 50 valid parameter sets is provided in Table 2, rank ordered by their proximity to the nominal fitness point (1,1). Given the current model structure, these values can be used as starting points to reproduce functionally realistic experimental data. As with many complex models, however, caution is warranted because the fitness of the model lacks robustness to experimentally realistic parameter variability. Simulations in the vicinity of these parameter sets are not guaranteed to be as experimentally valid (e.g., Fig. 5, *inset*).

This study suggests specific computational means to improve our understanding and representation of the mechanisms of motor-unit function. The results suggest reevaluating the gain functions describing the scaling of motor-unit discharge rate with excitation and the scaling of motor-unit force with discharge rate because both gain functions are influenced by peak discharge rates and recruitment ranges. The first gain function in the model assumes a linear scaling of motor-unit discharge rate with increased excitation (Kernell 1965; Schwindt and Calvin 1972), although important nonlinearities are known to exist (Hultborn et al. 2003; Kernell 1965). The second gain function assumes a sigmoid relation between motor-unit force and discharge rate (Rack and Westbury 1969), normalized to the contractile properties of the motor unit (Kernell et al. 1983), though motor-unit force and contraction time are not always correlated (Bigland-Ritchie et al. 1998). Moreover there is likely substantial context-dependent variability in the scaling described by both gain functions across motor units, which has not been implemented in the model. Also additional features not previously implemented in the model could be included to test whether its experimental validity improves. First, excitatory drive is not constant (as commonly modeled) but may vary systematically with increasing contraction intensity, for example, due to signal-dependent noise (Harris and Wolpert 1998). Thus an increase in variability of excitatory drive with increasing contraction intensity may attenuate the unrealistic drop-off found in variability of force at high force levels (Fig. 2*A*). Second, there may be a preferential distribution of low-threshold motor units deep in some muscles, inferred from the finding that slow, oxidative muscle fibers are more densely represented distant from the skin in some human muscles (Johnson et al. 1973; Lexell et al. 1983; cf. Edgerton et al. 1975). Thus at low levels of excitation, more active muscle fibers may be situated deeper in the muscle and thereby contribute less to the surface-detected EMG signal, potentially attenuating the unrealistic increase in EMG at low-to-moderate levels of force (Fig. 2*B*).

This work suggests iterations between modeling and experimental work will improve our understanding of the mechanisms of muscle function. For example, the feasibility of isolating a small number of axons from the motor neuron pool to a muscle in situ has been demonstrated (Day and Hulliger 2001; Perreault et al. 2003). This holds the promise of allowing the delivery of known activation patterns (e.g., specific recruitment and rate coding schemes) simulating voluntary contractions to validate the modeled gain functions for the generation of EMG and muscle force. In addition, alternative “muscle outputs” can be used as fitness criteria depending on the features of muscle function being investigated, such as power spectral densities. The present study is therefore an example of how computational modeling of neuromuscular systems could direct future experimental work to better understand the mechanisms that give rise to muscle function and EMG.

## GRANTS

This work was supported by National Science Foundation Grant 0237258 and National Institutes of Health Grants R21-HD-048566, R01-AR-050520, and R01-AR-052345.

## Acknowledgments

We thank Drs. Madhusudhan Venkadesan and Veronica J. Santos for helpful comments on the manuscript.

The contents of this work are solely the responsibility of the authors and do not necessarily represent the official views of the National Institute of Arthritis and Musculoskeletal and Skin Diseases (NIAMS), the National Institute of Childhood and Human Development (NICHD), the National Institutes of Health, or the NSF.

## Footnotes

↵1 The converse is also true: when a larger range of values is used to produce valid solutions, it means that the equations using those parameters are less critical to producing valid solutions.

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 © 2007 by the American Physiological Society