## Abstract

This study focuses on neuromuscular mechanisms behind ankle torque and EMG variability during a maintained isometric plantar flexion contraction. Experimentally obtained torque standard deviation (SD) and soleus, medial gastrocnemius, and lateral gastrocnemius EMG envelope mean and SD increased with mean torque for a wide range of torque levels. Computer simulations were performed on a biophysically-based neuromuscular model of the triceps surae consisting of premotoneuronal spike trains (the global input, GI) driving the motoneuron pools of the soleus, medial gastrocnemius, and lateral gastrocnemius muscles, which activate their respective muscle units. Two types of point processes were adopted to represent the statistics of the GI: Poisson and Gamma. Simulations showed a better agreement with experimental results when the GI was modeled by Gamma point processes having lower orders (higher variability) for higher target torques. At the same time, the simulations reproduced well the experimental data of EMG envelope mean and SD as a function of mean plantar flexion torque, for the three muscles. These results suggest that the experimentally found relations between torque-EMG variability as a function of mean plantar flexion torque level depend not only on the intrinsic properties of the motoneuron pools and the muscle units innervated, but also on the increasing variability of the premotoneuronal GI spike trains when their mean rates increase to command a higher plantar flexion torque level. The simulations also provided information on spike train statistics of several hundred motoneurons that compose the triceps surae, providing a wide picture of the associated mechanisms behind torque and EMG variability.

- torque variability
- plantar flexion torque
- force variability
- motor control
- neuromuscular model
- motoneuron pool model

motor tasks such as rising from a chair, standing, walking, reaching, pointing, grasping or moving an object always involve some degree of within-trial and between-trial variability. Torque variability during isometric contractions has been shown to be positively correlated with motor performance variability in healthy young adults (Magalhães and Kohn 2012; Mello et al. 2013). A similar association has been described in the elderly with postural instabilities or in pathologies of the neuromuscular system (Chu and Sanger 2009; Enoka 2008). Force (or torque) fluctuations generated during a muscle contraction are due to either random or deterministic mechanisms. The former include those related to the randomness of motor unit (MU) discharge times (Moritz et al. 2005), while a deterministic source of torque variability is associated with the phasic mechanical response of each muscle unit when it generates a twitch. However, when the whole set of MUs comprising a single muscle is taken into account, additional mechanisms should be considered, such as recruitment and correlated activity (synchronism) of the MUs (Jones et al. 2002; Yao et al. 2000).

Torque variability has been mostly studied during the performance of a steady isometric contraction, and the usual finding is that torque standard deviation (SD) increases with mean torque level. For some muscles torque SD was found to be proportional to the mean torque (Jones et al. 2002; Tracy 2007a), which is sometimes referred to as signal-dependent noise (Harris and Wolpert 1998). Nonetheless, other muscles exhibit a nonlinear relation between torque SD and the mean torque (Christou et al. 2002; Chu and Sanger 2009; Taylor et al. 2003a). The causes for such different relations can be studied from either a top-down approach, for example, using optimal control theories (Scott 2004), or from a bottom-up approach, in which the goal is to perform computer simulations of mathematical models of the neuromuscular system so as to investigate the underlying mechanisms that produce the observed scaling of torque variability with contraction intensity.

Jones et al. (2002a) and Jesunathadas et al. (2012) performed computer simulations of models of the neuromuscular system based on the phenomenological model developed by Fuglevand et al. (1993) to investigate the motor output variability produced by the first dorsal interosseous (FDI) and tibialis anterior muscles, respectively. They argued that the MU organization (the recruitment range and the orderly recruitment of MUs) was the key feature to obtain the relations between torque variability and mean torque that match the experimental results. Moritz et al. (2005) used a neuromuscular model of the FDI muscle also based on that developed by Fuglevand et al. (1993) to provide evidence that the MU discharge rate variability may also influence the scaling between isometric torque variability and muscle (FDI) contraction level. These and other studies in the literature attribute the general findings of scaling between isometric torque variability and muscle contraction level primarily to physiological features of the spinal cord and the muscle units.

To test even further the phenomenological neuromuscular model mentioned above (Fuglevand et al. 1993), Keenan and Valero-Cuevas (2007) added to the original Fuglevand model a mathematical description of the generation of electrical activity by the muscle fibers. They included in their analyses not only the torque variability dependence on mean torque level, but also the mean absolute value (MAV) of the electromyogram (EMG) as a function of mean torque level. They performed thousands of simulations (Monte-Carlo method) and showed that only few cases could reproduce experimental data for both torque variability vs. mean torque and EMG MAV vs. mean torque relations. In these few cases, the model's parameters were out of the physiological range, leading the authors to argue that the usage of a more refined model of neuronal mechanisms could produce more satisfactory results.

In the above-cited studies, the neuromuscular models employed for the analyses do not have a biophysical representation of spinal motoneurons (MNs), so that the spike trains are generated according to specific rules. As there is no provision to include the influence of neurons presynaptic to the MN pool, this phenomenological model can only provide information from the spinal cord downward to the muscle fibers. Therefore, the first hypothesis of this study was that a modeling approach based on biophysical MN models and with presynaptic inputs acting on the MNs through time-varying conductances would increase the biological realism of the neuromuscular model and improve the reproducibility of available torque and EMG data. Closely associated with this first hypothesis is the second hypothesis that the statistics of the premotoneuronal commands to the MN pool could contribute relevantly to the relations mentioned above involving torque SD and EMG MAV as a function of mean torque level.

Some considerations are needed with respect to the second hypothesis mentioned above. Frequently in computational neuroscience, stochastic point processes have been used as input commands to evaluate the output spike train statistics of neuron models (Lindner 2006). Poisson point processes have a constant coefficient of variation (CV) (equal to 1) irrespective of the mean intensity, and this may be interpreted as a condition in which the central nervous system has a lower control over the premotoneuronal command variability. Other stochastic point processes, such as Gamma point processes, have an additional degree of freedom (dependence on a second parameter) that allows the control of spike train variability independent of mean rate (Gabbiani and Cox 2010). Therefore, a corollary to the second hypothesis was that a Poisson point process input would not result in good enough fittings to all the experimentally obtained relationships, and a Gamma point process would.

This study aimed at investigating the neurophysiological mechanisms underlying torque variability during sustained isometric plantar flexion contractions, as well as the variability of the triceps surae (TS) muscles' EMGs, using both human experiments and computer simulations. The choice of the TS muscle group was due to its relevance in several motor tasks, such as upright standing and gait. Interestingly, few studies have systematically analyzed torque and EMG relationships with mean torque level for this muscle group. Computer simulations were carried out on a biologically based mathematical model of the neuromuscular system, which encompasses three pools (associated to each muscle of the TS) of biophysically based MN models subjected to stochastic synaptic inputs. The model parameters were constrained by known physiology from cats or humans (e.g., ionic channel dynamics, synaptic conductance dynamics, axonal conduction velocity distributions). Two different types of stochastic point processes (Poisson and Gamma with variable shape factor) were used to represent the premotoneuronal commands. The resulting discharge rate variability and correlated activity of the MUs emerged from the dynamic behavior of the whole neuronal network in response to the specific stochastic input spike trains to the MN pools. The results obtained from the computer simulations of the neuromuscular model should enable new inferences to be made on what mechanisms may be involved in the human neuromuscular system that controls isometric plantar flexion torque. Preliminary results have been presented as an abstract (Watanabe et al. 2012).

## MATERIALS AND METHODS

### Human Experiments

#### Participants.

The experiments were carried out in 11 healthy and right-footed volunteers (five women and six men) with age 28 (5) yr, height 1.69 (0.1) m and weight 67.5 (17.0) kg. All subjects gave informed consent, and all procedures were approved by the Human Ethics Committee of the Institute of Psychology at the University of São Paulo. The experiments were conducted in accordance with the Declaration of Helsinki.

#### Experimental setup.

The subjects were seated on a customized armchair designed for measuring ankle torque during isolated isometric plantar flexion contractions. The subject's right foot was firmly strapped to a rigid pedal connected to a strain gauge force transducer (N320, Transtec) that measured the compression force with a sensitivity of ∼96 N/V and with an output range of ±5 V. The force signals were low-pass filtered with a 500-Hz cut-off frequency and stored for later offline analysis. The knee was fully extended, while the hip and ankle were maintained at 130° and 90°, respectively.

#### Procedures.

The first step was to determine the torque value associated with each subject's maximal voluntary contraction (MVC) of the plantar flexor muscles. Subjects were asked to perform three MVCs (lasting 3 to 4 s each), with verbal encouragement and visual feedback of the exerted torque. The maximum torque value achieved across these trials was taken as the MVC force value.

After a period of rest, subjects were asked to keep the exerted isometric plantar flexion torque amplitude as constant as possible at target levels that ranged from 10% up to 80% MVC, in increments of 10% MVC. A custom-written program in LabView (National Instruments) provided the visual feedback of the exerted force on an LCD monitor. The target torque was provided as a horizontal line in the middle of the monitor, and the force exerted by the subjects as a thinner line progressing with time from left to right. Subjects were instructed to maintain their force on the target as accurately and as consistently as possible for 9 s.

Each trial began with a time interval of 5 s when the subjects visualized the output torque in a computer screen followed by a 4-s interval without visual feedback of the exerted torque (i.e., subjects were instructed to close their eyes after 5 s). The sequence of torque levels was randomized, and a resting period (duration depended on the previous force exertion) was allowed before the next trial. The gain of the signal display during the visual feedback phase as well as the distance between the subject and the display were kept constant for all subjects during the experiments. Three trials were executed for each particular torque level.

#### Data acquisition.

The force signals measured by the strain gauge were acquired in one of the channels of an A/D board (PCI-6221, National Instruments). The corresponding torque was computed from the measured force value and its moment arm (based on the dimensions and angles of the rigid structure associated with the pedal). The EMGs of the soleus (SOL), medial gastrocnemius (MG) and lateral gastrocnemius (LG) were obtained from surface electrodes positioned at standard Surface Electromyography for the Noninvasive Assessment of Muscle positions, except for the SOL in which one electrode was positioned 4 cm below the MG/LG junction and the second 2 cm distal to the first. The EMGs were amplified and filtered (5–1,000 Hz) by a Nihon-Kohden MEB-4200 system and fed to the same A/D board mentioned above. Both force and EMG signals were acquired at 2,000 samples/s.

#### Data analysis.

For the purpose of this study, only the data obtained from the periods without visual feedback were analyzed. The mean and the SD of the torque were measured from the torque signals (low-pass filtered with a 25-Hz cutoff frequency, zero-phase 4th order Butterworth). To eliminate slow drift components from the torque data, trend removal was done before SD computations. The detrending procedure consisted in computing the linear regression to the data and subtracting it from the signal. The EMG envelopes were estimated by passing the EMG absolute values through a zero-phase fourth-order Butterworth filter with a 5-Hz cutoff frequency. For both torque and EMG signals, a window of 3-s duration was used (i.e., the first and last 500 ms of the periods with closed eyes were discarded to avoid transients due to subject's adaptation and signal filtering). Hereafter, data will be presented as mean (SD) for all measured variables. When applicable, lower and upper limits of the 95% confidence interval are reported to allow appropriate comparisons between experimental and simulated data.

### Model of the Neuromuscular System

A model of the neuromuscular system was used to evaluate potential mechanisms underlying the results obtained experimentally. The model is freely available for access at http://remoto.leb.usp.br/ and follows the same structure proposed elsewhere (Cisi and Kohn 2008). Here, only part of the general structure was used (see Fig. 1). Briefly, the system adopted in the present study encompasses three MN pools, representing the motor nuclei of the SOL, MG, and LG muscles. Each motor nucleus includes mathematical models of individual type-specified MNs (see mathematical description below), with each MN commanding a muscle unit that generates MU twitches and MU action potentials (MUAPs). The number of MUs in each motor nucleus, estimated from human data (Johnson et al. 1973; McComas 1991), are as follows: 900 [800 slow MU (S)-type, 50 fast fatigue-resistant MUs (FR)-type and 50 fast fatigable MU (FF)-type MUs] for the SOL; 600 (300 S-type, 150 FR-type and 150 FF-type MUs) for the MG; and 260 (130 S-type, 65 FR-type and 65 FF-type MUs) for the LG. The sums of modeled twitches and MUAPs from all MUs of a muscle generate the whole muscle force (F) and surface EMG, respectively (see Fig. 1, *A* and *B*; mathematical models are presented further). The combination of the spike trains which are presynaptic to the MN pools driving the TS muscles are modeled by a set of independent homogeneous renewal point processes and will be referred to as the “global input” (GI). This GI is a superimposition of the spike trains from *1*) interneurons activated by descending drive, MNs, other interneurons and peripheral input; *2*) putative direct input from descending tracts; and *3*) activity from peripheral afferent axons. Each presynaptic spike train connects randomly to the MNs into the pool. The mean number of point processes that each MN receives is set by the connectivity ratio between the whole set of presynaptic spike trains and the MN pool. The interspike interval (ISI) statistics of the input stochastic point processes depend on the simulation protocol (see below).

#### MN model.

Conductance-based integrate-and-fire MN models of different types (i.e., S-, FR-, and FF-type) were adopted in each motor nucleus. The morphology of each model was simplified by adopting two compartments to represent the MN soma and its dendritic tree. Figure 1*C* shows the electrical circuit used to represent the passive membrane elements along with the active ionic channels. The axon is simply associated with a conduction velocity (varying along the pool), which causes a delay of the MN spikes to reach the muscle end-plate (the distance between the MN pool and the muscles is equal to 0.86 m). In the model, all the properties known to vary along the MN pool were parameterized following a piecewise linear distribution, and they are in accordance with physiological and anatomical data. Examples of such parameters are as follows: MN geometric and electrotonic properties, axon conduction velocity, and the voltage threshold for spike initiation. In this parameterization approach, the values of a given property were linearly distributed within a group of type-specified MNs, but the slopes of the linear segments could be different for different types of MNs. For instance, the parameterization adopted for the voltage threshold and for the passive membrane properties resulted in a rheobase distribution compatible with experimental findings (see Gustafsson and Pinter 1984; Powers and Binder 1985), i.e., many MNs have low rheobase values, and few have high ones. To add a nondeterministic aspect to the parameter values used in each simulation, zero mean Gaussian-distributed numbers were added to the mean values adopted for each MN's spike threshold (CV equal to 1%) and axon conduction velocities (CV equal to 5%). The same parameterization approach was adopted in Elias et al. (2012), and their Fig. 10 illustrates the rheobase distribution.

In each compartment, the membrane time course is described by *Eq. 1*, with *x* representing the evaluated compartment (i.e., soma or dendrite), and *y* the adjacent compartment (i.e., dendrite when the evaluated compartment is the soma or soma when the evaluated compartment is the dendrite); *C*_{x} is the membrane capacitance (in nF); *g*_{L,x} is the leakage conductance (in μS); *g*_{c} is the coupling conductance (in μS); *E*_{L} is the leakage Nernst potential (equal to the resting potential, in mV); *I*_{ion,x} is the active ionic current (in nA); and *I*_{syn,x} is the synaptic current (in nA). Passive elements are calculated using MN's geometric and electrotonic parameters, estimated from anesthetized cats (Fleshman et al. 1988; Zengel et al. 1985) as previously described by Cisi and Kohn (2008). For convenience, the resting potential is adopted as 0 mV, as done frequently (Hodgkin and Huxley 1952).(1)
where *V̇* is the membrane potential time derivative, and *V* is the membrane potential at time *t*. The somatic compartment comprises Na^{+} and fast K^{+} conductances, responsible for the genesis of spikes or action potentials, along with a slow K^{+} conductance, responsible for the after-hyperpolarization (AHP). The active somatic ionic current (*I*_{ion,S}) is given by *Eq. 2*, in which *ḡ*_{Na}, *ḡ*_{Kf}, and *ḡ*_{Ks} are the maximal ionic conductances (in mS/cm^{2}); *E*_{Na} and *E*_{K} are Na^{+} and K^{+} equilibrium potentials (in mV), respectively; *m*(*t*), *n*(*t*), and *q*(*t*) are state variables for Na^{+}, fast K^{+}, and slow K^{+} conductance activations, respectively; and *h*(*t*) is the state variable for the Na^{+} conductance inactivation. In this study, the MN models have a passive dendrite (i.e., *I*_{ion,d} = 0), representing a situation in which the spinal cord has no influence of descending monoaminergic pathways from the brain stem (Elias et al. 2012; Elias and Kohn 2013; Heckman et al. 2009).
(2)
where *V*_{s} is the membrane potential of the soma compartment. The state variables (*m*, *h*, *n*, and *q*) follow the pulse-based formalism proposed by Destexhe (1997), which reduces the computational load required to numerically solve the differential equations of thousands of individual MN models. Briefly, the first-order differential equations, in the form φ̇ = α_{φ}(1 − φ) − β_{φ}φ (with φ corresponding to a given state variable), can be analytically solved since the time course of forward and backward rates (α_{φ} and β_{φ}) can be approximated by rectangular pulses (with amplitudes α_{Φ} and β_{Φ}), triggered when the membrane potential crosses a given threshold (*V*_{th}). More details can be found in the appendix or elsewhere (Cisi and Kohn 2008; Elias et al. 2012; Elias and Kohn 2013). The active ionic channel parameter values (see Table 1) were chosen so that behaviors of a single model match experimental data from type-specified MNs of anesthetized cats (e.g., AHP amplitude, AHP duration, and *f*-*I* curves) (Schwindt and Crill 1984; Zengel et al. 1985). Differently from the previous model (Cisi and Kohn 2008), the voltage threshold (*V*_{th}) is chosen to match the experimental range of spike thresholds reported by Gustafsson and Pinter (1984). It is worth mentioning that although *V*_{th} is a parameter of the model, the time-varying conductances (see appendix for details) produce a dynamic behavior quite similar to that of a full Hodgking-Huxley model (Destexhe 1997), in which *V*_{th} emerges from the dynamics of the ionic channels.

Synaptic conductances are modeled by the kinetic model proposed by Destexhe et al. (1994). The total synaptic current is given by the sum of the currents generated by each of the *N*_{e} excitatory synapses acting on the dendritic compartment (*I*_{syn,d}) (*Eq. 3*). In this equation, *g*_{e,j}(*t*) is the conductance activated by the *j*-th excitatory synaptic input, given by the model mentioned above (see appendix for details), and *E*_{e} is the reversal potential. The value of maximal conductance (*ḡ*_{ē,l̄}) is 600 nS, and the reversal potential is set at +70 mV (with respect to the resting potential, set to 0 mV for convenience). These values are adopted so that individual excitatory postsynaptic potentials match experimental data from cats (Finkel and Redman 1983). To reduce the computational load when simulating a complex neuronal network, an algorithm (Lytton 1996) was used to account for the effects of multiple synaptic inputs acting on a given MN. Basically, this algorithm relies on the linearity and time-invariance of the kinetic model that describes a synaptic action on the dendritic compartment. A single lumped variable is used to represent the superposition of the arriving synaptic activations on a given compartment. Instead of treating all individual synaptic activations, the algorithm uses the weight of each synapse at a given time instant to update the state variable. More details may be found in the referred paper (Lytton 1996).
(3)
where *V*_{d} is the membrane potential of the dendritic compartment.

#### Muscle force model.

The twitch model used in the present study was based on a model developed by Raikova and Aladjov (2002) and has been recently used elsewhere (Contessa and De Luca 2013). This model was chosen over the classical linear second-order twitch model because its time course follows better some experimentally obtained muscle twitches. The impulse response for the twitch model is given by *Eq. 4*.
(4)
in which , *m* = *kT*_{c}, and *p* = *F*_{MU}^{MAX} exp[−*kT*_{C}(ln *T*_{C} − 1)], where *T*_{C} is the MU contraction time, *F*_{MU}^{MAX} is the MU twitch amplitude, and *T*_{hr} is half relaxation time of the MU twitch.

A muscle unit receives a train of action potentials coming from a MN, represented by a train of Dirac impulses, and generates a force (*F*_{MU}) that is given by the convolution between the train of impulses and the impulse response.

The ranges of values of the MU twitch amplitudes *F*_{MU}^{MAX} and parameters *T*_{C} and *T*_{hr} were based on human data (Bellemare et al. 1983; Buchthal and Schmalbruch 1980; Garnett et al. 1979) and are listed in Table 2.

The MU force *F*_{MU} generated by the twitch model was transformed by a nonlinear function (*Eq. 5*) that represents the muscle force saturation mechanisms (e.g., Ca^{2+} release saturation in the sarcoplasmic reticulum).
(5)
where the parameter *c* was set so that the MU force saturated at the desired frequency. The MU saturation frequencies were estimated from Enoka and Fuglevand (2001) (15–25 Hz for the S-type, 25–55 Hz for the FR-type and 55–65 Hz for the FF-type).

The muscle torque signal *M*_{MUS} (*Eq. 6*) was obtained by summing the force signals from *Eq. 5* for each of the *N* MUs and then multiplying the result by the value of the cosine of the pennation angle α_{MUS} (α_{SOL} = 25°, α_{MG} = 17°, and α_{LG} = 8°), the value obtained in the force-length curve *fl*_{MUS} (*fl*_{SOL} = 0.6, *fl*_{MG} = 1, and *fl*_{LG} = 1) and the muscle moment-arm *m*_{MUS} (*m*_{SOL} = 0.0413 m, *m*_{MG} = 0.0418 m, and *m*_{LG} = 0.0429 m). The values of these biomechanical parameters were calculated from equations given in Menegaldo et al. (2004). The total TS torque is given by the sum of the respective *M*_{MUS} of each muscle that comprises the TS.
(6)

#### EMG model.

As stated above, the simulated EMG was given by the sum of MUAPs from the active MUs. These MUAPs were modeled as first- (*HR*_{1}) or second-order Hermite-Rodriguez functions (*HR*_{2}) (*Eqs. 7* and *8*, respectively) (Lo Conte et al. 1994), which are randomly attributed to a given MU.
(7)
(8)
with *A*_{MUAP} representing the MUAP amplitude, λ_{MUAP} the MUAP duration, *t*_{spike} the time of occurrence of a MU spike at the endplate, and *u*(*t*) the step function. The adopted values of *A*_{MUAP} and λ_{MUAP} are shown in Table 2. These values were chosen to represent intramuscular MUAP time courses recorded from humans. In addition, nonlinear amplitude attenuation and time scaling were applied to these signals according to the position of the MU within the muscle cross section relative to a pair of surface electrodes located over the muscle belly (Cisi and Kohn 2008; Fuglevand et al. 1992). The positions of the MUs were chosen randomly within a muscle cross section. Simulated EMG signals were band-pass filtered (1st-order Butterworth, cutoff frequencies of 50 and 500 Hz) and had a low-level white Gaussian noise added for more realistic visualization (see e.g., Fig. 2, *D*–*F*). This EMG model is a simplification of the known mechanisms underlying MUAP generation and volume conductor filtering (Farina and Merletti 2001; Keenan et al. 2005; Negro and Farina 2012). However, the simplification does not appear to be a critical issue here since the power spectra of simulated and experimental EMG time series (rectified and low-pass filtered) resulted quite similarly (see Endnote).

#### Simulation protocol.

All simulations were performed with a 0.05-ms time step using a Runge-Kutta integration method. The GI was represented by 400 independent renewal point-processes. Each renewal point-process represents the superposition of the activity of a collection of premotoneuronal neurons. Each of these processes activated 30% of the MNs, randomly chosen among the three MN pools. Two different types of GI spike train statistics were used to drive the MN pools to evaluate possible differences in the resulting torque and EMG statistics.

In a first step, the maximum torque of the neuromuscular system, equivalent to the MVC, was determined. This was done using Poisson point processes in the GI with 4-ms mean ISI. With this mean ISI value, the maximum firing rate achieved by a given MN was about 30 spikes/s, a typical value obtained from human MNs during MVC (Enoka and Fuglevand 2001). All the torque and EMG data obtained from the simulations at lower torque values than MVC were processed in the same way as in the experiments with humans and were normalized by the respective values obtained from the simulations at MVC.

Different mean ISI values of the GI were chosen to obtain simulated torque levels similar to the experimental ones (10% to 80% MVC). The same mean GI ISI values were used (range: 10.5 ms at 10% MVC to 5 ms at 80% MVC) for each point process type: Poisson and Gamma process. Poisson point processes are renewal processes with exponentially-distributed ISIs, while Gamma point processes are renewal processes with ISIs following a Gamma probability density function, with a given shape factor (or order). In our simulations, the shape factor of the Gamma point processes varied from 2, corresponding to 80% MVC, to 7, corresponding to 10% MVC. For the intermediate values, the orders were 5 (20% MVC), 4 (30% MVC), 4 (40% MVC), 4 (50% MVC), 3 (60% MVC) and 2 (70% MVC). These values were obtained on the basis of several simulation runs and comparisons with experimental data from individual MUs reported by Person and Kudina (1972).

Each simulation resulted in signals with 5-s duration, but only the last 3 s (60,000 data points) were used to compute the mean and SD of the simulated torque and EMG envelope. These data were processed in the same manner as the experimental data. At each torque level, the mean ISI and the ISI SD of each MN were computed from the times series of the ISIs using Matlab (The MathWorks). In addition, for each MN and for the GI the mean firing rate and its SD were estimated from extra simulations (each lasting the equivalent of 30 s of “neural time”) using the asymptotic firing rate variability given by *Eq. 9* (Cox 1962).
(9)
where, SD_{ISI} is the ISI SD, μ is the mean ISI and μ_{3} is the skewness of the ISI distribution.

The population synchrony index (PSI) presented in Yao et al. (2000) was modified to account for a large number (>21) of active MUs and used to compute the synchronization of the MNs spikes at each torque level for each MN pool. The common input strength (CIS) (Nordstrom et al. 1992) was computed for a single long-lasting simulation (200 s duration) at 10% MVC to compare with experimental data (Keen et al. 2012; Mochizuki et al. 2005).

### Statistical Analysis

Correlation and linear regression analyses were performed to assess the relationship between EMG/torque computations and mean torque levels (for both experimental and simulated data). An analysis of covariance (ANCOVA, with the mean torque levels treated as the covariate) was used to compare the slopes of the regression lines among the experimental and simulated data (Zar 1999). When the null hypothesis of equality between slopes was rejected, the procedures of Dunnett's test were used to perform multiple comparisons to identify which slopes differed from the experimental one. Correlation tests and contrasts were performed using a 0.05 significance level.

## RESULTS

Experimental results are presented in the first part of this section. These are followed by simulated data covering torque and EMG envelope fluctuations as a function of mean torque level for different input statistics. The analyses of the simulated data also include MN and GI spike train statistics that reflect aspects of the stochastic behavior of a whole MN pool. In the simulations, the TS torque is considered as a representative for plantar flexion torque, since some authors report that this muscle group is largely responsible for ankle extension (Delp et al. 1990; Fukunaga et al. 1996; Murray et al. 1977). When applicable, experimental and simulated data (e.g., mean torque vs. torque SD) are presented in the same graphs to facilitate the comparisons and analyses. Examples of experimental and simulated TS plantar flexion torques and EMG envelope signals (in gray) from the SOL muscle are shown in Fig. 2. These are being shown here merely as qualitative results: the more quantitative analyses will be presented in what follows. Also, a more detailed presentation of the model parameterization can be found elsewhere (see Endnote) to permit a more quantitative appraisal of the appropriateness of the models employed.

### Experimental Results

The gray lines in Fig. 3*A* show the torque SD as a function of the exerted mean torque levels plotted in a log-log scale, for each of the 11 subjects. The thicker dark line is the linear fit for the relationship between the group average values of torque SD (crosses) and the mean torque level. Table 3 presents the mean torque and torque SD computed for the group data (with their respective SD). Figure 3 shows the CV of torque as a function of mean torque (the gray lines are data from the subjects, the thickest line is the sample average). For all subjects, the SD of force output was positively and significantly correlated with the exerted force values (*P* values ranging from 10^{−5} to 0.03), while, for most of the subjects, the CV of torque was relatively constant (mean value around 1.30%) across the range of exerted torques (no significant correlation, i.e., *P* > 0.05, was found between the mean CV of torque and the mean exerted torque values). The slope (and lower and upper limits of the 95% confidence interval) of the regression line relating the torque SD and the mean exerted torque values (in a log-log scale) resulted 0.80 (0.63–0.97). Or, in other words, torque SD (σ_{T}) is related by a power law to mean torque (μ_{T}), i.e., σ_{T} = α·(μ_{T})^{0.80}, where α is a positive constant value.

Figure 4, *A*–*C*, shows the log-log relationships between the mean EMG envelopes (of the SOL, MG and LG muscles) and mean plantar flexion torque level, whereas Fig. 4, *D*–*F*, shows the relations between the EMG envelope SDs of the three muscles and the mean plantar flexion torque level. Gray lines are individual data, and the thickest continuous line is the best linear fit. Regardless of the muscle, all subjects had the mean and the SD of the EMG envelope significantly correlated with the mean exerted torque value (*P* values ranging from 1.4 × 10^{−8} to 0.012). Table 4 presents the slopes (and lower and upper limits of the 95% confidence interval) of the regression lines of the mean and SD of the EMG envelope as a function of mean exerted torque (in a log-log scale).

### Simulation Results

#### Influence of GI type on torque and EMG variability.

Figure 3, *A* and *B*, shows the simulated torque SD and CV, respectively, vs. mean torque data (triangles for Poisson inputs and circles for Gamma inputs). The dashed and dash-dotted lines are the regression curves for the torque SD as a function of the mean torque for Poisson and Gamma input point processes, respectively. The slopes of the torque SD vs. mean torque curves and the respective confidence intervals are 0.45 (0.16–0.64) for Poisson point process and 0.74 (0.62–0.80) for Gamma point process inputs. The ANCOVA showed significant differences among slopes obtained experimentally and by computer simulations [and *F*_{(2,21)} = 8.0104, *P* = 0.0026]. Multiple comparisons (Dunnett's tests) showed that when the GI was modeled by Gamma point processes, no significant differences were observed between the simulated slopes and those obtained experimentally [*q′*_{(21,2)} = −0.8716, *P* > 0.05]. Note that the slopes of the relations are within the experimental 95% confidence intervals. Conversely, when the GI was modeled by Poisson point processes, the resulting slopes from the regression relating torque SD and mean torque were significantly lower than those obtained experimentally [*q′*_{(21,2)} = −3.7867, *P* < 0.01]. Note that the slopes of the regression curves of the simulated data are outside the boundaries of the experimental counterpart. Additionally, when Poisson point processes were used to model the GI acting on the MN pool, one may qualitatively observe that the magnitude of torque variability (measured as the torque SD) was higher than in the experimental data, mainly at the lower mean torque levels. This is more dramatically observed in the relation between the torque CV and mean torque (Fig. 3*B*). The simulated CV values for the Poisson point process were always higher than the experimental data; in addition, there was a clear reduction in the torque CV as the mean torque increased when the GI was modeled as Poisson processes, differently from what was observed experimentally (see *Experimental Results*, above). On the other hand, for the Gamma point processes, the CV vs. mean torque curve resulted approximately flat with a mean CV within the experimental range (mean ∼1.3%).

It may be concluded that the simulations with GI modeled by Gamma point processes resulted in a significantly better fit to the average (between-subjects) experimental result (Fig. 3) (see discussion for further explanations).

Figure 4, *A*–*C*, shows that the relations between mean EMG envelope and mean torque for the three muscles resemble reasonably well the experimental data using either Poisson or Gamma point processes to represent the GI. Similarly, for both GI models the relationships between the SD of the EMG envelopes and mean torque levels (Figure 4, *D*–*F*) resulted within the experimental range. A more quantitative comparison is given in Table 4, which shows the slopes of all the regression lines, for both experimental and simulated data. The ANCOVA showed significant differences among slopes obtained experimentally and by simulations only for the relationships regarding the mean EMG envelope of the SOL and MG muscles [*F*_{(2,20)} = 5.3778, *P* = 0.0135 and *F*_{(2,20)} = 4.1218, *P* = 0.00317 for the SOL and MG, respectively]. In these cases, multiple comparisons (Dunnett's tests) showed that, when the GI was modeled as Poisson point processes, the regression slope from the relationship between SOL EMG and mean torque level was significantly higher than that found experimentally [*q′*_{(20,2)} = 3.2209, *P* < 0.01], whereas the slope from the relationship between MG EMG and mean torque level was significantly lower than that obtained experimentally [*q′*_{(20,2)} = −2.5552, *P* < 0.01]. When Gamma processes were used as GI, no significant differences were found between the simulated and experimental slopes of SOL and MG EMG [*q′*_{(20,2)} = 1.1204, *P* > 0.05 and *q′*_{(20,2)} = −1.6979, *P* > 0.05 for SOL and MG, respectively]. Note that, except for the mean EMG envelope of the SOL and MG muscles obtained when the GI was modeled as Poisson processes, all other simulated values were within the 95% confidence intervals of the experimental data (Table 4), suggesting that the relations between EMG MAV, EMG SD and mean torque are less influenced by the statistics of the neuronal input to the MN pool than the relations between torque variability (SD and CV) and mean torque level, as described above.

#### Power spectrum and scaling of the GI.

It is known that the output variability is influenced by the spectral content of the GI (Negro and Farina 2011). The mean spectral content of the GI was estimated from the dendritic conductance fluctuations of a given MN due to the activations from all the presynaptic inputs. Figure 5 presents normalized power spectra (Welch's method) of net dendritic synaptic conductance signals of a representative MN when either Poisson or Gamma point processes are used to model the GI and at three different mean torque levels. The choice of a single MN dendritic conductance signal as a representative for this analysis is valid since the 400 presynaptic axons are randomly selected to make synapses to each MN of the pool at each simulation run. Three different point process intensities were used in each case, associated with three contraction levels (10%, 50% and 80% MVC). For 10% to 50% contraction levels the power of the synaptic input resulted lower for the Gamma point processes compared with the Poisson, especially in the frequency range compatible with the operation of the whole system (i.e., lower than 50 Hz). For high contraction level (80% MVC), both stochastic GIs generated synaptic drives to the representative MN with similar spectral distributions.

Mean firing rate and firing rate SD for a single point process of the GI were calculated from the respective mean ISI and ISI SD (*Eq. 9*) at each torque level. Figure 6 shows the relation between firing rate SD and mean firing rate (lower axis), as well as firing rate SD and mean torque (upper axis). Data points corresponding to a Poisson point process (triangles) resulted in a linear scaling (in log-log scale) with a slope equal to 0.50. Conversely, the linear regression for a Gamma point process (data in circles) resulted in a slope equal to 1.34.

#### Variability and synchronism of MN spike trains.

All the following analyses were applied to the hundreds of spike trains associated with the firings of the MNs of the LG muscle and the firings of the GI. Figure 7, *A* and *B*, shows relations between different MN pool spike train statistics corresponding to all eight target mean torque levels. When the MNs were excited by Poisson point processes (Fig. 7*A*, triangles) higher levels of MN ISI SD were observed compared with the case of GI modeled as Gamma processes (Fig. 7*A*, circles). The difference was more evident for higher values of mean ISIs (i.e., at lower MN firing rate). For lower values of mean ISIs (i.e., higher discharge rates), both Poisson and Gamma point processes resulted in MN pool ISI SD values around 5 ms (see *inset* in Fig. 7*A*). For visualization purposes, the above-cited relation between the ISI SD and mean ISI of MN spike trains were adjusted by a quadratic function (*R*^{2} > 0.99), as drawn in Fig. 7*A* in continuous and dashed lines. Irrespective of the GI model, ISI histograms (data not shown) for high mean ISI MN spike trains were asymmetric with positive skewness, whereas ISI histograms for low mean ISI MN spike trains were approximately Gaussian. Additionally, Fig. 7*B* shows that the MNs' firing rate SD (computed using *Eq. 9*) resulted higher when the GI were Poisson point processes (triangles). Continuous and dashed curves plotted in Fig. 7*B* for visualization purposes are based on the quadratic fittings shown in Fig. 7*A* transformed to rate SD using expression (16). Interestingly, the overall firing rate SD decreases with the mean firing rate. Figure 7*C* shows MN ISI CV as a function of the respective mean ISI for both Poisson (triangles) and Gamma (circles) point processes as input.

Figure 8*A* shows the firing rate SD average values computed from all of the recruited MNs at each torque level. At each level, the mean values of firing rate SD when the GI followed a Poisson point process were higher than those when the GI followed a Gamma point process. Additionally, for both GI models there was a decreasing trend in the firing rate SD as the torque level increased. This is in accordance with Fig. 7*B*, because the MNs discharge faster for higher intensity contractions. Figure 8*B* shows the average values for the ISI CVs calculated from all of the recruited MNs at each torque level. The mean values of ISI CVs when Poisson point processes were used were higher at all torque levels than when the GI was modeled as Gamma point processes.

The values of the overall synchronization index PSI (Yao et al. 2000) as a function of target torque values are shown in Fig. 9. The PSI attained very similar values for both types of point processes used to model the GI for the three MN pools (the slopes and intercepts are given in the figure caption). In addition, the CIS (Nordstrom et al. 1992) value found for 20 pairs of MNs was 0.26 spikes/s for Poisson GI and 0.18 spikes/s for Gamma GI, which is in accordance with the data available in the experimental literature for the SOL muscle (Keen et al. 2012; Mochizuki et al. 2005). In part, this result justifies the use of a 30% connectivity (see *Simulation protocol*) between the independent random processes that constitute the GI and the MNs. The CIS was also computed for the whole TS, resulting in 0.30 spikes/s for Poisson GI and 0.08 for Gamma GI.

## DISCUSSION

In this study, two different statistical descriptions (homogenous Poisson and Gamma point processes) were used to model the overall premotoneuronal input (GI) acting on a multiscale neuromuscular model that describes mathematically the MN pools and muscle units of the TS. Computer simulations provided the corresponding relationships between plantar flexion torque SD and mean torque levels, as well as EMG parameters (mean and SD of EMG envelopes) as a function of the mean torque levels. These were compared with the corresponding relations obtained from experiments on healthy young subjects.

### Influence of GI Type on Torque and EMG Variability

Poisson point processes have been largely employed to model the firing statistics of independent neural ensembles (Burkitt 2006; Elias et al. 2012; Elias and Kohn 2013; Taylor and Enoka 2004) since they approximately represent a superposition of independent point processes, at least locally in time (Lindner 2006). Lee et al. (1998) and Stein et al. (2005) found that for single cortical neurons in monkey motor cortex the mean slope of the relation firing rate SD vs. mean firing rate (in log-log scale) is 0.50. Stein et al. (2005) reported that the same value was also obtained for a population of cortical neurons. These experimental data suggest that a plausible approximate model for the descending spike train statistics is a Poisson point process.

However, the descending commands may not act directly on the MN pool. In fact, some experimental findings suggest that, at least for lower limb muscles, part of the descending cortical command passes through indirect oligosynaptic pathways (Nielsen and Petersen 1995). Simulation results from our laboratory (unpublished observations) showed that, when conductance-based interneuron models (see Cisi and Kohn 2008) receive synaptic inputs from independent Poisson point processes, the resulting output spike train ISI histogram for each interneuron follows approximately a Gamma distribution, and consecutive ISIs are uncorrelated (i.e., as a first approximation, the output spike trains of interneurons may be modeled as a Gamma renewal point process). Another point to emphasize is that the presynaptic command arriving at the MN pool (the GI in our model) originates not only from direct or indirect cortical sources, but also from, for instance, afferent signals from sensory receptors (e.g., muscle spindles, Golgi tendon organs, joint and cutaneous receptors). These afferent signals cannot be considered independent sources since different afferent fibers may discharge with similar or correlated patterns in response to a common mechanical input (Inbar et al. 1979; Windhorst 1977). This would make the Poisson point process description less appropriate to model the GI spike train statistics. Synthesizing, if there is a strong cortico-motoneuronal effect on the TS MN pools, then the Poisson point process model approximation could be appropriate. If the descending commands act preferably through spinal cord interneurons and/or the effect of afferent inputs is significant, then the Gamma point process description may be more adequate. Nonetheless, the adoption of Gamma point processes was somewhat arbitrary in this study. Maybe other point processes (e.g., ISIs with log-normal distributions) could lead to similar results, but we have not evaluated this issue. Another relevant observation is that the GI was modeled exclusively by excitatory commands to the MNs. This simplification is justifiable from previous results (Elias and Kohn 2013; Rudolph and Destexhe 2006), which showed that models such as those used here (passive-dendrite MNs) integrate inhibitory postsynaptic potentials linearly with excitatory postsynaptic potentials to produce the global synaptic current to the MNs. In other words, with the inclusion of inhibitory presynaptic inputs to the MNs, the intensity of the excitatory inputs would have to be increased to keep the MNs firing at their reported rates, but the resulting spike train statistics of all the MNs would remain similar to those described here.

From a modeling point of view, one of the limitations of the Poisson point process as a model for representing spike trains (or their superposition) is that it has a constant ISI CV and a constant Fano factor, both being equal to 1. This means that for any intensity (mean rate) of the point process, its relative variability is always the same. On the other hand, lower ISI CVs may be obtained for Gamma point processes by adjusting the Gamma distribution shape-parameter.

Figure 3 shows that the torque SD vs. mean torque relation obtained when the GI was modeled as Gamma point processes (with shape-parameter decreasing for higher GI firing rates) was within the experimental range, whereas when the GI was modeled as Poisson point processes the relation indicated higher variability, mainly at lower-intensity contractions. This suggests that the premotoneuronal GI ISI CV cannot be independent of the mean ISI of the GI. Actually, the GI stochastic processes have to present a low ISI CV when the mean ISI value is large (i.e., at lower torque levels) and high ISI CV when the mean ISI decreases (i.e., at higher torque levels). Actually, Fig. 6 synthesizes this from a different viewpoint: the GI firing rate SD has to vary in a steeper way as a function of GI mean firing rate (or torque level) than what a Poisson process may provide. Therefore, all these results suggest that the scaling of plantar flexion torque SD with the mean torque is also dependent on the GI statistics at the premotoneuronal level. It is worth mentioning that even though a relatively complex twitch model was used (see materials and methods), qualitatively similar conclusions (superiority of variable-order Gamma point processes as the GI) were reached (data not shown) when the twitch model was the classic impulse response of a critically damped second-order system (Milner-Brown et al. 1973).

The experimentally obtained increasing trend of torque SD as a function of the mean torque level varying from 10% to 80% MVC (Fig. 3) is in accordance with previous reports for a smaller range of torques [2.5% to 10% MVC (Shinohara et al. 2006)] and for a similar range but measured with less resolution (Tracy 2007b). It is interesting to note that, despite the use of 3-s contraction periods in the present research, the values of torque SD were similar to the values found by Tracy (2007) that analyzed 5- to 10-s contractions. Our study was the first to provide an explicit mathematical relation for TS plantar flexion torque variability scaling: torque SD is proportional to (mean torque)^{φ}, with φ = 0.80 (0.63–0.97). Data from other muscles resulted in different exponents for the same power function (Chu and Sanger 2009; Hamilton et al. 2004) or in different nonlinear relations (e.g., see Christou et al. 2002). The shape of this relation must depend on many neurophysiological and muscular features, such as the statistics of the GI (as shown in the present work), the muscle twitch waveform, the number of MUs (Hamilton et al. 2004), among other things.

As can be seen in Table 4, all the EMG envelope measures (except for the SOL and MG EMG envelope means when the GI was modeled as Poisson processes) are within the 95% confidence intervals of the experimental data, which provides additional validation of the models employed (especially when Gamma point processes were used to represent the GI). These findings suggest that the scaling observed between the EMG of each muscle of the TS and the plantar flexion torque may be associated with the different features of these muscles, e.g., proportions of slow- and fast-type MUs, numbers of MNs and ranges of parameter values of the MU models along the pools. However, further investigations are needed to clarify the exact mechanisms. It is worth noting that both Poisson and Gamma point processes provided a fine representation of the scaling between EMG (mean and SD) and torque, suggesting that the statistics of the GI do not dramatically affect the relations between EMG envelope parameters and mean torque level. This may be related to the filtering process used to obtain the EMG envelopes, which attenuates high-frequency components of the absolute value of the EMGs, possibly masking differential effects from each of the stochastic point processes used to represent the GI.

It is interesting to note that in the simulations performed by Keenan and Valero-Cuevas (2007), it was very difficult to satisfy at the same time the relations between mean EMG vs. mean force and force SD vs. mean force. Even when both relations were in agreement with the experimental data, they were obtained using model parameters known to be less realistic from a biological standpoint. Therefore, they proposed the usage of a more refined model of the neural mechanisms. In this vein, Heckman and Enoka (2012) suggested that a model that allows a representation of the synaptic integration by the MNs would lead to a more appropriate representation of the neural drive from the motor pool to the muscles and hence would accurately reproduce the scaling of motor variability found experimentally. Indeed, a special feature of the modeling approach of the present paper is that each MN in the pool is based on dynamic models (ionic channels, membrane capacitance, etc.), the synapses act through time-varying conductances and the GI to the MN pool is modeled as stochastic point processes with adjustable statistics, all of which were probably important to achieve the matching between simulation and experimental data, in terms of both torque and EMG.

### Power Spectrum and Scaling of the GI

Besides having a constant ISI CV and a constant Fano factor, an additional limitation of the Poisson point process for representing spike trains is that it has a flat power spectrum, regardless of its intensity. On the other hand, different ISI CVs and power spectra may be obtained for Gamma point processes. Figure 5 exemplifies the differences in estimated synaptic conductance power spectra associated with Poisson (approximately constant for all frequencies) and Gamma point processes (a valley at low frequencies for Gammas with shape-parameters 7 and 4, corresponding to Fig. 5, *A* and *B*). This allowed the control of the degree of variability of the input spike trains representing the GI. Consequently, distinct motor output statistics (MN spike train ISIs and rate SDs; torque and EMG SD and CV) could emerge from the adoption of those two different types of stochastic point processes. In this vein, Fig. 7*A* shows that when the MNs were excited by Poisson point processes, the ISI SD of MNs showed higher values than those observed experimentally in humans by Person and Kudina (1972) (see their Fig. 6a), especially at the higher values of mean ISI (i.e., lower MN firing rate). In contrast, the use of Gamma processes with shape-parameter varying according to the target torque led to a reasonably good replication of the experimental relation between the mean vs. SD of ISIs reported in the literature (e.g., Person and Kudina 1972 and Clamann 1969).

The scaling of the firing rate SD with the mean firing rate of the GI is a square root when a Poisson point process model is used for the different contraction levels (Fig. 6). However, using *Eq. 10*, the Gamma point processes with decreasing order for higher firing rates produced a relation between firing rate SD and mean firing rate of the GI with a higher slope (1.34) and with a smaller variability for lower-intensity contractions compared with the Poisson case. In the power spectral analysis, this is reflected as a lower power of the net excitatory synaptic conductance observed at low- to mid-intensity contractions when the GI was modeled as a Gamma point process (see Fig. 5). This might have some impact on optimal control studies that usually rely only on the experimentally measured torque SD vs. mean torque level behavior.

### Variability and Synchronism of MN Spike Trains

From *Eq. 9* we were able to estimate the firing rate variability of several MNs discharging at different rates and for different contraction levels. These results (depicted in Figs. 7*B* and 8*A*) show that when the GI is modeled as Poisson point processes the MN pool rate variability is higher than when Gamma processes are used instead. Taylor et al. (2003) and Moritz et al. (2005) stated in their studies that the firing rate variability (although they did not show the firing rate explicitly) of the MNs commanding a given muscle is a fundamental characteristic to reproduce the motor output variability (i.e., torque or force variability) obtained experimentally. What we found here was in a similar vein but with an added view of causality, i.e., when the MN pool discharge variability (both in terms of ISI or firing rate) is generated by appropriate presynaptic MN inputs, the relationships between torque variability (either SD or CV) and mean torque are reasonably well reproduced by the neuromuscular model. Jones et al. (2002) showed that the MN pool discharge variability had little influence on the proportional scaling of motor output variability, albeit marked changes were found in the magnitude of the torque variability. These results may be somewhat expected since they adopted a constant ISI CV as contraction intensity increased. They also found that the MU organization (i.e., the recruitment range and the orderly recruitment of MUs) was the key feature to obtain experimentally observed responses. In fact, in our study the structure and dynamics of the MN pool model and the muscle units they innervate tried to capture the main biophysical features of the biological system: for instance, the range of twitch amplitudes and contraction times; the passive membrane properties of MNs; the spike thresholds along the pool. Using this biologically-based structure we were able to better reproduce the experimental data when the GI was represented by Gamma point processes with shape-parameters dependent on target torque level. Therefore, these results suggest that the discharge variability across a MN pool is an important factor contributing to torque variability, and that it is greatly influenced by the statistics of the presynaptic inputs that drive the MN pool. It is noteworthy that the biologically based structure adopted in the model was essential to achieve a torque variability scaling compatible with the experimental data, which could only be observed when the GI was modeled by Gamma point processes with shape-parameter dependent on mean GI firing rate.

What we discussed so far in terms of differences between GI statistics has a major implication for lower-intensity contractions. The results of Fig. 8 show that, as the simulated contraction level increases, the MN pool firing rate SD and ISI CV tend to decrease, achieving a plateau approximately at 40% MVC. However, in spite of this decrease (which is more evident with Poisson point processes as the GI), the torque variability starts to increase for mid-to-high intensity. This raises the question of what mechanisms are responsible for the increase in torque SD as a function of mean torque. One mechanism is associated with the higher peak twitch amplitudes of higher-threshold MNs that are recruited at higher forces. At the same time, several authors discussed the relevance of the synchronism of MUs for torque variability (Taylor et al. 2002; Taylor and Enoka 2004; Yao et al. 2000). Figure 9 shows that the population synchrony of each MN pool increases as a function of the mean torque level, but with no important differences between the models adopted for the GI. Therefore, these two mechanisms (i.e., recruitment of high-threshold MNs and increased synchronism) seem to appropriately counterbalance the decreased firing rate and ISI variability of the MNs in the pool at higher force levels.

### Concluding Remarks

To conclude, this study suggests that the experimental relations linking plantar flexion torque and EMG variability with mean torque level are associated with quantitative requirements on both the GI statistics and the structural/biophysical features of the MU pools involved. Probably the increased biological realism, together with the appropriate statistics in the GI, led to a good reproduction of the experimental findings, for both the torque SD and EMG envelope relations, besides providing interesting hypotheses as to the mechanisms involved.

It is worth mentioning, however, that this study was focused solely on the TS neuromuscular system. Therefore, the modeling and the simulated results are applicable only to this particular neuromuscular group. For more general conclusions on the interplay of the dynamics of the MN pool/biomechanics and the GI statistics, a theoretical approach or a wide-scale Monte-Carlo investigation would be needed. In addition, other aspects of the biological system [e.g., viscoelastic muscle properties, low-frequency oscillations due to proprioceptive feedback loops and descending commands (Dideriksen et al. 2012; Negro et al. 2009)] should be explored in further studies, since they may also influence the torque and EMG relations associated with motor behavior.

## GRANTS

The research was funded by grants from Fundação de Amparo a Pesquisa do Estado de São Paulo (FAPESP) 2006/54185-8 and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) (303313/2011-0; Brazilian Funding Agencies) to A. F. Kohn. R. N. Watanabe, F. H. Magalhães, L. A. Elias, and E. M. Mello were supported by scholarships from FAPESP (nos. 2010/12934-0, 2011/13222-6, 2009/15802-0, and 2006/60399-0, respectively), V. M. Chaud was supported by scholarship from CNPq (no. 132776/2011-1). R. N. Watanabe and L. A. Elias are now supported by scholarships from FAPESP (nos. 2011/21103-7 and 2013/10433-1).

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

## ENDNOTE

At the request of the authors, readers are herein alerted to the fact that additional materials related to this manuscript may be found at the institutional website of the authors, which at the time of publication is: http://remoto.leb.usp.br/remoto/Publications/Supplemental_01.pdf. These materials are not a part of this manuscript, and have not undergone peer review by the American Physiological Society (APS). APS and the journal editors take no responsibility for these materials, for the website address, or for any links to or from it.

## AUTHOR CONTRIBUTIONS

Author contributions: R.N.W. and A.F.K. conception and design of research; R.N.W., L.A.E., and E.M.M. performed experiments; R.N.W., F.H.M., and L.A.E. analyzed data; R.N.W., F.H.M., L.A.E., V.M.C., and A.F.K. interpreted results of experiments; R.N.W., L.A.E., and V.M.C. prepared figures; R.N.W., F.H.M., L.A.E., V.M.C., and A.F.K. drafted manuscript; R.N.W., F.H.M., L.A.E., V.M.C., and A.F.K. edited and revised manuscript; R.N.W., F.H.M., L.A.E., V.M.C., E.M.M., and A.F.K. approved final version of manuscript.

## ACKNOWLEDGMENTS

We are grateful to S. A. Miqueleti for the technical assistance during human experiments.

## Appendix

#### Pulse-Based Method for Computation of Ionic Channels State Variables (Destexhe 1997)

When the threshold is crossed, the rectangular pulses representing the forward and backward rate constants turn on and remain active during a given time duration. While the pulse is on, *Equations A1* and *A2* give the time course of activation (*m*, *n* and *q*) and inactivation (*h*) variables, respectively. In these equations, *α*_{Φ} and *β*_{Φ} are forward and backward rate constants, respectively; *t*_{on} is the time at pulse activation, and φ_{a}^{0} and φ_{i}^{0} are the values of φ_{a}(*t*) and φ_{i}(*t*) at that time. The pulse duration is set at 0.60 ms.
(A1)
(A2)

Before and after the pulse, *Eqs. A3* and *A4* give the time course of activation and inactivation variables, respectively. In these equations, α_{Φ} and β_{Φ} are forward and backward rate constants, respectively; and *t*_{off} is the time at pulse deactivation.
(A3)
(A4)

#### Computation of Synaptic Conductance (Destexhe et al. 1994)

The time course of a synaptic conductance may be represented by a first-order kinetic model described by *Eq. A5*, with *r*(*t*) representing the fraction of bound receptors; α and β are the forward and backward constants for transmitter binding, respectively; and [*T*] is the concentration of transmitter released into the synaptic cleft.
(A5)

Destexhe et al. (1994) assumed that [*T*] occurs as a pulse, since the concentration of neurotransmitter in the cleft rises and falls rapidly. With this assumption, *Eq. 5* may be analytically solved, so that during the pulse [*T*] = *T*_{max} (set to 1 mM) and *r*(*t*) are given by *Eq. 2*. After the pulse, [*T*] = 0 and *r*(*t*) is given by *Eq. A6*.
(A6)
with *r*_{∞} = α*T*_{max}/(α*T*_{max} + β); τ_{r} = 1/(α*T*_{max} + β); and *t*_{0} is the occurrence time of a presynaptic event. Here, α and β were set to 0.50 ms^{−1} and 2.50 ms^{−1}, respectively.
where *t*_{1} is the time at the end of the pulse. Pulse duration (*t*_{1}–*t*_{0}) was set to 0.20 ms, so that the synaptic saturation occurs for input frequencies higher than 5 kHz.

The total synaptic conductance is given by the product between *r*(*t*) and the maximum synaptic conductance (*ḡ*_{e,l}).

- Copyright © 2013 the American Physiological Society