Motor neurons are often assumed to generate spikes in proportion to the excitatory synaptic input received. There are, however, many intrinsic properties of motor neurons that might affect this relationship, such as persistent inward currents (PICs), spike-threshold accommodation, or spike-frequency adaptation. These nonlinear properties have been investigated in reduced animal preparation but have not been well studied during natural motor behaviors because of the difficulty in characterizing synaptic input in intact animals. Therefore, we studied the influence of each of these intrinsic properties on spiking responses and muscle force using a population model of motor units that simulates voluntary contractions in human subjects. In particular, we focused on the difference in firing rate of low-threshold motor units when higher threshold motor units were recruited and subsequently derecruited, referred to as ΔF. Others have used ΔF to evaluate the extent of PIC activation during voluntary behavior. Our results showed that positive ΔF values could arise when any one of these nonlinear properties was included in the simulations. Therefore, a positive ΔF should not be considered as exclusive evidence for PIC activation. Furthermore, by systematically varying contraction duration and speed in our simulations, we identified a means that might be used experimentally to distinguish among PICs, accommodation, and adaptation as contributors to ΔF.
- motor neuron
- computer model
mammalian motor neurons have long been considered as prototypical neurons (e.g., Eccles 1950) that integrate synaptic inputs from many sources and are thought to generate action potentials at rates proportional to the overall excitatory input (Granit et al. 1963). There are, however, a number of properties of motor neurons that may disrupt a simple relationship between synaptic input and firing rate output. In this report, we consider three such properties and their effects on motor neuron output, namely, persistent inward currents (PICs), spike-threshold accommodation, and spike-frequency adaptation.
First described in the 1970s (Schwindt and Crill 1977), PICs have now been established as a robust intrinsic property of motor neurons (Bennett et al. 1998, 2001b; Conway et al. 1988; Hamm et al. 2010; Hounsgaard et al. 1984, 1988a; Hounsgaard and Kiehn 1989; Hultborn et al. 2003; Lee and Heckman 1998a, 1998b, 1996; Schwindt and Crill 1980, 1982; Turkin et al. 2010). PICs are thought to be mediated by voltage-activated dendritic (Bennett et al. 1998; Booth et al. 1997; Carlin et al. 2000b; Hounsgaard and Kiehn 1993; Lee and Heckman 1996) and somatic channels (Ballou et al. 2006; Moritz et al. 2007) that require neuromodulators, such as serotonin, to be enabled (Hounsgaard et al. 1988a). PICs may last for many seconds (Lee and Heckman 1998a) and are composed predominantly of an L-type calcium current (Carlin et al. 2000a; Hounsgaard and Kiehn 1989), but persistent sodium currents may play a role in specific circumstances (Li et al. 2004a; Manuel et al. 2007). PICs have been implicated in synaptic input amplification (Hultborn et al. 2003; Lee et al. 2003; Prather et al. 2001), as well as self-sustained firing, whereby motor neuron activity continues even after the cessation of excitatory synaptic drive (Conway et al. 1988; Crone et al. 1988; Hounsgaard et al. 1984, 1988a; Lee and Heckman 1998b). As such, PICs have been proposed to facilitate prolonged activity in motor neurons needed for maintenance of posture (Heckman et al. 2003; Kiehn and Eken 1998), to boost excitability of motor neurons during locomotion (Heckman et al. 2003), and to provide the primary source of depolarizing current to motor neurons during normal motor behaviors (Heckman et al. 2005).
Spike-threshold accommodation refers to a progressive increase in the amount of depolarizing current required to bring a neuron to action potential threshold as the rate of rise of current decreases (Hill 1936; Wigton and Brink 1944). In motor neurons, the current required to initiate activity can increase by as much as twofold with slow rates of rise of injected (Araki and Otani 1959; Bradley and Somjen 1961; Burke and Nelson 1971; Sasaki and Otani 1961, 1962; Schlue et al. 1974a, 1974b) or synaptic current (Eccles 1946; Homma et al. 1970). The mechanisms underlying accommodation have not been completely elucidated but appear to be primarily related to sodium channel inactivation (Bradley and Somjen 1961; Hodgkin and Huxley 1952; Kernell 2006; Schlue et al. 1974c; Vallbo 1964).
Finally, spike-frequency adaptation represents a time-dependent diminution in firing rate in response to a constant current input, and in motor neurons can be divided into early and late phases (Gorman et al. 2005; Powers et al. 1999; Sawczuk et al. 1995; Spielmann et al. 1993). Early adaptation involves the first few spikes at the outset of activity in response to a step increase in depolarizing current and is responsible for a marked decay in firing rate from initial high values, often in excess of 100 impulses/s (e.g., Sawczuk et al. 1995). Late adaptation, on the other hand, is characterized by a slow exponential drop in firing rate proportional to the magnitude of the injected current leading to decreases in firing rate by as much as 40–60% over a period of 30 s of continuous activity (Button et al. 2007; Kernell and Monster 1982; Sawczuk et al. 1995). Although a variety of mechanisms have been proposed to account for early (Granit et al. 1963; Miles et al. 2005; Powers et al. 1999; Sawczuk et al. 1995) and late adaptation (Barraza et al. 2009; Partridge and Stevens 1976; Sawczuk et al. 1997), definitive mechanisms underlying adaptation are yet to be fully identified.
Because intrinsic properties such as those described above can markedly affect the way neurons transform input current into spiking output (e.g., Hamm et al. 2010), it is essential to characterize these properties to fully understand how neurons process information. Whereas such intrinsic properties have been studied extensively in motor neurons of reduced preparations, little is known about how these properties actually interact with synaptic input to shape the output response of a motor unit pool during natural behaviors. This is because it is presently impossible to quantify the extent of synaptic input to motor neurons in awake behaving animals (Fuglevand et al. 2006; Heckman et al. 2009).
In an attempt to overcome this difficulty, Kiehn and Eken (1997) and Gorassini and colleagues (1998, 2002a, 2002b) developed a clever method that involved using changes in firing rate of a low-threshold motor unit to represent changes in synaptic drive to a population of motor neurons. This method was then used to assess the presence and magnitude of PICs in motor neurons using the following approach. The activities of pairs of human motor units were recorded during voluntary ramp increases and decreases in muscle force. The firing rate of the first recruited, “control” unit was used as an indicator of synaptic drive. Because PICs are thought to be activated rapidly and fully near the time a motor neuron reaches spiking threshold (Bennett et al. 1998) and to remain relatively stable for several seconds thereafter (Lee and Heckman 1998a), any subsequent changes in firing rate that occur in the control unit after recruitment were taken to reflect changes in the distributed synaptic input acting on the motor neuron pool.
As a consequence, recruitment and derecruitment of the second, higher threshold, “test” unit could be registered in terms of changes in firing rate of the control unit. Once the test unit is recruited, and if a PIC is then instigated, there are two sources of depolarizing current sustaining activity in the test unit: that arising from synaptic input and that due to the intrinsic PIC. Subsequently, to terminate activity in the test unit, the total depolarizing current must drop below that associated with spiking threshold. Accordingly, the derecruitment threshold of the test unit should be associated with a lower overall synaptic drive (reflected in a lower firing rate of the control unit) compared with that at recruitment because of the added contribution of the PIC. As such, a difference in firing rate of the control unit at the times of recruitment and derecruitment of the test unit, referred to as ΔF, should be related to the magnitude of PIC activation.
An expanding number of reports consider positive ΔF measurements (i.e., higher firing rate of control unit at recruitment compared with at derecruitment of a test unit) as prima facie evidence for PIC activation in motor units (ElBasiouny et al. 2010; Gorassini et al. 1998, 2002a, 2002b, 2004; Heckman et al. 2005, 2008, 2009; Kiehn and Eken 1997; Mottram et al. 2009; Oya et al. 2009; Powers et al. 2008; Stephenson and Maluf 2010; Udina et al. 2010). It seems feasible, however, that other nonlinear properties of motor neurons, such as accommodation or adaptation, might also lead to positive ΔF values. To evaluate this possibility, we used an updated version of a motor unit pool model (Fuglevand et al. 1993) in which we selectively included or omitted each of these nonlinear properties in simulations of muscle contractions. Because these individual properties depend differentially on rate of rise of excitation and on time, we also carried out simulations in which we varied contraction speed or duration while selectively enabling each nonlinear property. For each simulated condition, the magnitudes of ΔF were measured across populations of active motor units. Our results indicate that positive ΔF values can arise due to any one of these nonlinear properties, and therefore, a positive ΔF should not be considered exclusive evidence of PIC activation. However, we go on to show that changing the duration of ramp contractions may provide a means to distinguish among the probable causes of ΔF.
Motor unit spike times and isometric muscle force were simulated using a motor unit population model (for details, see Fuglevand 1993). To this model, we added representations of PICs, spike-threshold accommodation, and spike-frequency adaptation. Our overall approach involved simulating the phenomena, rather than the detailed biological mechanisms that shape the output behavior of a population of motor units. In this way, each motor neuron and motor unit could be represented individually with relatively few parameters. After a brief description of the original model, the physiological justification and procedures used to simulate the nonlinear features are described. The simulations presented in this report involved motor unit activities associated with isometric forces up to about 15% of the maximum force of the simulated muscle. This was done to simulate the conditions typically used in human motor unit experiments to assess ΔF. As such, the simulations primarily involved activities of low-threshold motor units only. The model was implemented in the MATLAB environment (The MathWorks, Natick, MA), and the software is available on request.
Motor neuron pool.
A population of 120 motor neurons was modeled to represent the pool innervating a typical human muscle. Each motor neuron received the same level of excitatory drive. Excitatory drive can be considered more or less equivalent to the net synaptic current reaching the soma and spike-initiating zone. Motor neuron spike thresholds were determined by an exponential function where many motor neurons had low thresholds and few had high thresholds for recruitment. The range of recruitment thresholds (RTE) across the pool was 50-fold. The minimum firing rate (MFR) was set to 8 impulses/s (imp/s) for all motor neurons (Monster and Chan 1977; however, see Moritz et al. 2005). If excitatory drive exceeded the assigned threshold for a particular motor neuron, i, then that motor neuron generated action potentials based on the following firing rate (FR) equation: (1) where g is a constant gain factor for all motor neurons and E(t) is a time-varying excitatory drive function. The units of excitatory drive are normalized such that 1.0 unit of excitation represents the approximate level of excitation needed to recruit the lowest threshold motor neuron. For all simulations reported, g was set to 1.0 imp/s per excitation unit. According to this equation, motor neuron firing rate increased linearly until a predetermined maximum rate was reached: 35 imp/s for the first recruited motor neuron and down to 25 imp/s for the last recruited motor neuron. The normal variability in motor neuron discharge was emulated by adjusting the spike times to simulate a Gaussian random process with a 20% coefficient of variation (however, see Moritz et al. 2005).
Motor unit twitch force was modeled as an impulse response of a second-order critically damped system. The lowest threshold motor unit innervated the smallest number of muscle fibers and generated the smallest twitch force, whereas the last recruited motor unit innervated the largest number of muscle fibers and produced the largest twitch force. Motor unit twitch forces varied exponentially as a function of recruitment order over a 100-fold range such that most motor units were assigned relatively small twitch forces and relatively few, high-threshold units produced very large forces. Twitch contraction times were inversely related to twitch force and ranged from 30 to 90 ms such that the first recruited motor unit had the longest contraction time. A variable gain factor amplified motor unit twitch forces based on discharge rate to mimic the well-known sigmoidal relationship between discharge rate and isometric force. The total force produced by the muscle was a linear sum of the individual motor unit forces.
To this basic model, three modifications to the spiking behavior of the neurons were added: PICs, spike-threshold accommodation, and spike-frequency adaptation. These parameters could be selectively activated for each simulation.
Persistent inward currents.
A PIC provides an intrinsic source of excitation to motor neurons that is thought to boost the synaptic drive the motor neuron receives (Lee et al. 2003; Prather et al. 2001). Because calcium-mediated PICs appear to be the dominant form of PIC in motor neurons (Hounsgaard and Kiehn 1989), we focus our simulations primarily on this type and hereafter refer to them simply as PICs. In the simulations, PICs were added as an excitatory “current” that summed linearly with excitatory drive once a motor neuron was active. Because the activation threshold for PICs occurs near the time that motor neurons reach spiking threshold (Lee and Heckman 1998a; Li and Bennett 2003; Li et al. 2004b), particularly if motor neurons are activated synaptically (Bennett et al. 1998), PICs were therefore modeled to activate in each motor neuron at the moment the neuron was recruited (Fig. 1).
Intracellular recordings from motor neurons indicate that PICs may remain active for ten s or more (Carlin et al. 2000a; Lee and Heckman 1998a; Li and Bennett 2003), and a hyperpolarizing pulse may be needed to stop the self-sustained firing attributed to plateau potentials mediated by PICs (Hounsgaard et al. 1988a). Whereas PICs recorded in high-threshold motor neurons may decrement with time, PICs in low-threshold motor neurons remain relatively stable over several seconds (Lee and Heckman 1998b). Therefore, for purposes of the present simulations involving mostly low-threshold motor neurons, PICs were initially simulated as a steady level of excitatory current that remained active until the motor neuron was derecruited (see Fig. 1).
The reported magnitudes of PICs (in nanoamperes) vary widely. Certainly, this is likely due to differences in experimental preparations used, particularly in terms of the neuromodulatory state of the preparation. However, this can also occur for motor neurons within a pool under similar experimental conditions. For example, Lee and Heckman (2000) showed only a weak tendency for higher threshold motor neurons supplying cat triceps surae to exhibit larger PICs. On the other hand, Hamm et al. (2010) found no systematic relationship between PIC magnitude and indicator of recruitment threshold in rat hindlimb motor neurons. However, when PIC magnitude is expressed relative to recruitment threshold current (Hamm et al. 2010), low-threshold motor neurons systematically exhibit larger normalized PICs than do high-threshold motor neurons (Hamm et al. 2010; Lee and Heckman 2000). For example, depending on how PICs were measured, the lowest threshold motor neurons in the Hamm et al. study had PICs on the order of 1.5–2.5 × threshold current, whereas PICs for the highest threshold motor neurons were ∼0.2 × threshold current. Likewise, assuming rheobase currents of 3 and 30 nA for the lowest and highest threshold motor neurons, respectively (Gustafsson and Pinter 1984), and a roughly uniform PIC current of 7.8 nA across all motor neurons in the standard neuromodulatory state in the Lee and Heckman study (2000), the normalized PIC would be ∼2.6 × threshold current for the lowest threshold neurons and 0.26 × threshold current for the high-threshold neurons.
To partially mimic these values in our simulations, the absolute magnitude of the PIC was assigned the same value of 2.0 excitation units for all motor neurons. Therefore, for the lowest threshold motor neuron (with a recruitment threshold of 1.0 excitation unit), this value represented a normalized PIC amplitude of 2.0 × threshold excitation, whereas for a motor neuron with a 10-fold higher threshold, normalized PIC amplitude was 0.2 × threshold excitation. These values, therefore, are similar to those associated with the experimental findings of Hamm et al. (2010) and Lee and Heckman (2000).
It is not entirely clear as to how quickly a PIC activates. Some evidence suggests PICs activate in more or less an all-or-none manner (Bennett et al. 1998; Heckman et al. 2005; Hounsgaard and Kiehn 1989; Lee and Heckman 1999), whereas other findings indicate that PICs are engaged more slowly, sometimes requiring a second or more to activate fully (Bennett et al. 2001a, 2001b; Li and Bennett 2003; Li et al. 2004a). Furthermore, activation characteristics may depend on the magnitude of excitation (Hounsgaard and Kiehn 1989; Hultborn et al. 2003). Since there is no established consensus in the literature, our nominal simulations modeled PICs in the simplest way, i.e., as being fully active at their recruitment (Fig. 1). This assumption has also been adopted by Gorassini et al. (2004) associated with the ΔF measurement.
The threshold for spiking may accommodate, that is, it may vary based on the rate of excitation delivered to the motor neuron (Araki and Otani 1959; Bradley and Somjen 1961; Burke and Nelson 1971; Homma et al. 1970; Sasaki and Oka 1963; Sasaki and Otani 1961, 1962; Schlue et al. 1974a, 1974b, 1974c). Although some reports have shown higher threshold motor neurons to exhibit the most accommodation (Burke and Nelson 1971; Sasaki and Otani 1961), other reports have shown low-threshold motor neurons to respond with the greatest degree of accommodation (Bradley and Somjen 1961) or for there to be no systematic difference in accommodation across motor neurons of differing thresholds (Schlue et al. 1974a). The magnitude of accommodation, namely, the current needed to bring a neuron to threshold during ramp current injection compared with that during a step current injection (i.e., the rheobase current, Irh), has been shown to range between 1.1 and 3.3 × Irh with an average value of ∼1.8 × Irh during 1-s ramps (Schlue et al. 1974a). It should be pointed out that accommodation has also been reported based on an increase in the actual voltage threshold (Homma et al. 1970) at the onset of spiking. For our model, we have defined accommodation based on the amount of additional current needed to bring a neuron to action potential threshold.
Based on the above, we simulated accommodation in the following way. For every motor neuron, threshold increased above the nominally assigned recruitment threshold excitation (equivalent to rheobase current) as a simple inverse function of the rate of rise of excitatory drive (Fig. 2). Because there is no clear consensus as to whether there are systematic differences in accommodation with threshold, we assigned the relative increase in threshold with increased rate of rise of excitation to be the same for all motor neurons. As described below, we tested three different rates of rise of excitation in these simulations: ∼1.6, 3.2, and 6.4 excitation units/s. As shown in Fig. 2 for two example motor neurons (MN 5 and MN 40), the slowest rate of rise tested in this study would lead to about a 63% increase in threshold for both these neurons. This level of accommodation is within the limits of those often reported experimentally and may even partially underestimate the extent of accommodation because only high rates of rise of excitation have been tested experimentally. Finally, it should be pointed out that accommodation as modeled in this study only affected the first spike in a train of spikes. Presumably, the rapid repolarization and afterhyperpolarization immediately following the first spike relieves most fast sodium channels from inactivation and thereby returns the spike-initiating zone to its near full complement of activatable sodium channels for the subsequent spike (however, see discussion).
Firing rate adaptation.
Motor unit activity rarely exhibits the dramatic drop in firing rate at the outset of muscle contractions indicative of early adaptation. Because of this, and also because our simulations involved gradual rather than abrupt increases in excitatory drive, early adaptation was not included in the model. Instead, we focused on late adaptation, namely, a progressive, and relatively slow, reduction in firing rate in response to a steady level of depolarizing current. Late adaptation, which we refer to hereafter simply as adaptation, is a prominent feature of motor neurons (Button et al. 2007; Gorman et al. 2005; Hounsgaard et al. 1988b; Kernell and Monster 1982; Sawczuk et al. 1995; Spielmann et al. 1993; Turkin et al. 2010; Viana et al. 1995). The magnitude of adaptation tends to be larger with greater levels of depolarizing current (Gorman et al. 2005; Kernell and Monster 1982; Sawczuk et al. 1995). In addition, although some research suggests that the degree of adaptation is more pronounced in larger compared with smaller motor neurons (Kernell and Monster 1982), other work suggests that there is no influence of recruitment order on adaptation magnitude (Gorman et al. 2005; Sawczuk et al. 1995). Therefore, adaptation magnitude was simulated in this study to be independent of excitation threshold.
Although it is not clear what leads to late adaptation in neurons (Brownstone 2006), some work suggests adaptation may be due to a non-calcium-dependent, slowly increasing potassium conductance (Barraza et al. 2009; Partridge and Stevens 1976; Sawczuk et al. 1997). Partridge and Stevens (1976) simulated this conductance as a simple time-dependent rising exponential with a relatively long time constant. We adopted their approach in this study by including an exponentially rising outward “current” that was subtracted from the excitatory drive function to yield the net excitation acting on a modeled neuron (Fig. 3A). The extent of this adaptation current, A, for any neuron i was a function of both the time since motor neuron recruitment, rti, and the excitation level, E(t), namely, (2) where e is the base of the natural logarithms and τ is the time constant. On the basis of experimental observations of Sawczuk et al. (1995) and Gorman et al. (2005), we assigned τ to have a value of 22 s. The parameter q in Eq. 2 designates the maximum value of the adaptation current, and it depended on excitation level above threshold [i.e., E(t) − RTEi] as (3) The parameter ϕ was selected to match the magnitude of adaptation for different levels of excitation as reported by Kernell and Monster (1982) and was assigned a value of 0.67. Finally, the parameter d was included to account for the observation that the minimum firing rate at recruitment appears not to be the absolute minimum firing rate that a motor unit can sustain (Barry et al. 2007). Since a lower firing rate at derecruitment relative to recruitment in human motor units may be due to adaptation (Barry et al. 2007; De Luca et al. 1982), the firing rate was allowed to decay with time below the initially specified minimum firing rate by a small amount determined by d. In the present simulations, d was set to a value of 2 imp/s, similar to the values reported by Moritz et al. (2005) and Barry et al. (2007). Equation 3 was evaluated at each time step of the simulation (see below) to solve for q based on the level of excitatory drive at each moment. The resulting instantaneous value of q was then substituted into Eq. 2 to estimate the adaptation current associated with that time step. Figure 3B shows examples of the simulated instantaneous firing rates of a motor neuron in response to three different levels of steady excitatory drive when adaptation, as represented by Eqs. 2 and 3, was included in the model.
A time resolution of 1 ms was used for all simulations. Isometric force and motor unit firing were initially simulated under four different conditions: control (no nonlinearities added to the model), inclusion of PICs, inclusion of accommodation, and inclusion of adaptation. In each case, the excitatory drive function linearly increased for 5 s and then decreased linearly back to baseline in 5 s, giving rise to a triangular force profile. The peak value of excitatory drive was selected such that the peak force under the control condition reached a value of about 15% of maximum force and was associated with recruitment of about 80 motor units. Maximum force (MF) was determined as the average force generated when all motor units in the pool were driven at supramaximal excitation levels.
These four conditions were further tested at two other rates of rise and fall of excitation such that peak force (i.e., 15% MF) was achieved in 2.5 s (fast condition) or in 10 s (slow condition). Furthermore, we also carried out simulations for which we 1) varied the parameters characterizing PICs, 2) varied contraction duration independently of rate of rise of force, and 3) concurrently included all three nonlinear features in the model.
The dependent measure of primary focus here was ΔF (see Introduction), which is calculated as the difference in firing rates of a low-threshold motor unit between the instances at which a higher threshold unit is recruited and then derecruited during a triangular contraction. For each simulation condition tested in this study, 10 low-threshold motor units were used as control units (the first 10 recruited units) and 10 higher threshold motor units were used as test motor units (the 31st through 40th units recruited). By comparing each lower threshold control unit with each higher threshold test unit, we created 100 motor unit pairs for each simulation condition. Firing rates of control units were calculated as five-point moving averages. Linear interpolation of the five-point moving average was then used to identify the firing rate associated with the specific times at which higher threshold units were recruited and derecruited. The firing rate value at derecruitment was subtracted from that at recruitment to determine ΔF for each low-threshold–high-threshold motor unit pair. In addition, we also determined the recruitment and derecruitment forces for every unit as the muscle force at the instant of the first or last spike, respectively, during each triangular contraction. Last, motor unit firing rates at recruitment and derecruitment were calculated as the average firing rate over the first and last five interspike intervals, respectively, for each simulated contraction.
A one-way analysis of variance (ANOVA) was used to determine whether ΔF differed across the four primary simulation conditions (control, PICs included, accommodation included, adaptation included). In addition, ANOVA was carried out to determine whether ΔF varied with simulated rates of rise or duration of contraction. Paired t-tests were performed to determine whether recruitment and derecruitment firing rates and forces were different within simulation conditions. Significance was set at P < 0.05.
The goal of these simulations was to determine whether properties of motor neurons other than PICs could lead to a positive ΔF and, should this be the case, to ascertain a way to distinguish which of these properties might best account for ΔF observed experimentally. Figure 4A shows results from the primary set of simulations for each of the four test conditions. For each condition, in addition to the predicted whole muscle force, the firing patterns associated with a representative motor unit pair (a low-threshold control unit and a higher threshold test unit) are shown. These units indicate the firings of unit 5 and unit 40 in the recruitment order.
Even though the excitatory drive function was identical for each trial, the peak force varied somewhat across the four conditions (15.2, 17.8, 13.8, and 14.4% maximum force for the control, PIC, accommodation, and adaptation conditions, respectively; Fig. 4A). As would be expected, the PIC condition produced the largest force, because every motor neuron, once recruited, received an additional source of excitation in the form of the simulated PIC. On the other hand, the lowest peak force was produced under the accommodation condition. This occurred because for the intermediate rate of rise of excitation used in these simulations, the recruitment threshold of all motor neurons was slightly increased above Irh when accommodation was enabled (see Fig. 2). Therefore, for the accommodation condition, fewer total motor units were recruited (75) compared with the other three conditions (84).
For the example pair of units shown in Fig. 4A, positive ΔF values, although small, were evident for the PIC, accommodation, and adaptation conditions but not for the control condition. Across all 100 pairs of motor units, ANOVA indicated a significant effect (P < 0.001) of simulation condition on ΔF (Fig. 4B). Post hoc analyses (Holm-Sidak method) indicated that individual mean ΔF values for each condition were significantly different (P < 0.05) from that of all other conditions. Furthermore, a simple location t-test indicated that whereas ΔF for the control condition was not significantly different from zero (P > 0.05), ΔF was positive and significantly different (P < 0.001) from zero for the PIC, accommodation, and adaptation conditions.
We also examined firing rates and force levels at recruitment and derecruitment of the motor unit population for each simulation condition. In the control condition, firing rates at recruitment and derecruitment were not different from one another (P > 0.05, Fig. 5A). By comparison, the firing rates of simulated motor units under PIC, accommodation, or adaptation conditions were all lower at derecruitment compared with recruitment (Fig. 5A, P < 0.001 for all comparisons). Experimentally, a lower firing rate at derecruitment compared with recruitment is often observed in human motor units (e.g., Clamann 1970; De Luca et al. 1982; Milner-Brown et al. 1973b; Pascoe et al. 2011) and has been explained as a result of adaptation (De Luca et al. 1982) or as a result of PIC activation (Gorassini et al. 2002a). As far as we are aware, the role that accommodation might play in differences in firing rates at recruitment and derecruitment has not been discussed in the literature. In terms of the present simulations, higher firing rates at recruitment occurred with accommodation because the onset of spiking was delayed and, as a consequence, was associated with a higher level of excitatory drive compared with the control condition. Because firing rate was directly related to excitatory drive (Eq. 1), this led to higher firing rates at the outset of spiking. Accommodation had no effect, however, on motor unit firing during the descending phase of excitation.
In addition, we compared recruitment to derecruitment force thresholds. Under control simulations, the average force at recruitment was lower than that at derecruitment (P < 0.001, Fig. 5B). A similar result was also found for the simulation that included adaptation. On the other hand, recruitment and derecruitment forces were not different from one another when PICs or accommodation were included in the simulations (P > 0.05, Fig. 5B). These disparate results parallel to some degree findings in experimental work. For example, some investigators (De Luca et al. 1982; Fuglevand et al. 2006; Milner-Brown et al. 1973a; Patten and Kamen 2000) found smaller forces at recruitment of human motor units compared with derecruitment, similar to the results reported in this study for the control and adaptation conditions. This type of finding has been explained primarily on the basis of muscle mechanics (De Luca et al. 1982; Fuglevand et al. 2006; Milner-Brown et al. 1973a). Other studies, however, have shown recruitment force to be larger than derecruitment force (Gorassini et al. 2002a; Person and Kudina 1972; Romaiguere et al. 1989, 1993). This type of result has been provisionally explained as being due to PIC activation (Gorassini et al. 2002a) or to enhanced activity of antagonist muscles during the descending phase of triangular contractions (Fuglevand et al. 2006; Patten and Kamen 2000).
It might be noted that when examining average recruitment force values across experimental conditions (filled bars, Fig. 5B), there was no difference between the accommodation and control conditions. At first glance, this result might seem unexpected given that accommodation was associated with higher excitation thresholds needed to recruit motor units compared with control. This, however, did not translate into higher forces at recruitment for the following reasons. With accommodation enabled, and at a given level of excitation during the rising phase, fewer motor units were recruited compared with control conditions. With fewer units recruited, the overall muscle force was diminished. Indeed, for this reason, peak muscle force was lowest for the accommodation condition (see Fig. 4A) despite identical excitatory drive across conditions. Therefore, lower muscle force at a given level of excitation in the accommodation condition tended to offset the higher level of excitation needed to recruit a unit such that the net effect was little average change in recruitment thresholds (in terms of force) across the population.
Effect of changing simulated contraction speed on ΔF.
Although the PIC condition yielded the largest ΔF values (see Fig. 4B), the physiological significance of differences in the absolute magnitudes of ΔF across simulation conditions should be considered with caution. This is because the specific parameter values selected to represent the different nonlinear properties, while influenced by experimental findings, were chosen somewhat arbitrarily. Perhaps the most salient finding thus far is that positive ΔF values were observed for all three conditions involving nonlinear properties. Therefore, we next sought to determine whether it might be possible to distinguish among these properties as contributors to ΔF by varying characteristics of the muscle contraction. Firing rate adaptation, for example, is a time-dependent property of motor neurons, so if the total duration of the contraction changes, adaptation-related ΔF should also vary as a function of total contraction duration. Spike-threshold accommodation, on the other hand, depends on the rate of rise of excitation. As such, if the excitation rate of rise changes, accommodation-related ΔF should also be affected. PICs, however, as they were simulated in this study, should be relatively insensitive to changes in contraction speed or contraction duration. Therefore, to test these ideas, we reran the simulations using slower and faster rates of rise of excitation than used in the initial simulations. The “slow” simulations involved symmetric triangular excitation profiles that peaked at 10 s (rate of force rise = 1.5% MF/s) and “moderate” speed peaked at 5 s (rate of force rise = 3.0% MF/s), whereas for the “fast” speed, peak excitation occurred in 2.5 s (rate of force rise = 6.0% MF/s). In all cases, the value of the peak excitation was the same. Although the moderate-speed conditions were the same as those used for Fig. 4, we reran those simulations nevertheless. Because of the imposed randomness in the specific discharge times for individual motor units, the average values of ΔF varied to some degree for each simulation under identical conditions.
Figure 6A illustrates simulated muscle force and activities of an example motor unit pair for the control condition at the three different contraction speeds. For the example motor unit pair in Fig. 6A, there was a small positive ΔF for the slow condition and negligible ΔF for the moderate and fast conditions. For all 100 motor unit pairs, there was no significant effect of contraction speed on ΔF for the control condition (P > 0.05, Fig. 6B). Likewise, when PICs were enabled in the model, ΔF was not significantly different across contraction speeds (P > 0.05). However, when accommodation was included in the model, there was a significant effect of contraction speed on ΔF (P < 0.001). Post hoc analysis (Holm-Sidak method) indicated that ΔF decreased significantly for each increase in contraction speed in the accommodation condition. Likewise, simulations with adaptation enabled also showed a significant effect of contraction speed on ΔF (P < 0.001) with a significant diminution with each increase in contraction speed (Fig. 6B). Overall, these results suggest an experimental means to distinguish the possible contribution of PICs vs. accommodation or adaptation to ΔF. If ΔF remained stable across contraction speeds, then this would be consistent with a primary role of PICs. If on the other hand, ΔF diminished with contraction speed, then this would be more likely due to the contributions of accommodation or adaptation. Such a reduction in ΔF with increased contraction speed has recently been demonstrated experimentally in human motor units (Stephenson and Maluf 2011).
Effects of varying PIC parameters.
Although these results seemed promising to distinguish among possible contributors to ΔF during voluntary contractions, it is possible that PIC behavior is more complex than our initial assumptions. In particular, PICs might decrease in magnitude over time (Hounsgaard and Kiehn 1989; Lee and Heckman 1998a) or might require some time to be fully expressed (Bennett et al. 2001a, 2001b; Li and Bennett 2003; Li et al. 2004a). Therefore, we varied these PIC parameters to test their effect on ΔF with two further sets of simulations involving variations in contraction speed.
In the first, PICs were activated instantaneously but then decayed linearly such that there would be a 50% reduction in PIC amplitude over 10 s, similar to the maximum rate of decay shown by Lee and Heckman (1998a) for partially bistable motor neurons. In this case, ANOVA indicated a significant (P < 0.05) effect of contraction speed on ΔF (Fig. 7A). Although post hoc comparisons were not significant, overall there was a tendency for ΔF to increase with contraction speed, opposite to the effect when accommodation or adaptation were activated (Fig. 6B).
For the second set of simulations, we addressed the possibility that PICs do not activate instantaneously. For these simulations, PICs increased linearly from zero to their assigned maximum value in 500 ms and then decayed at the same rate as above. Interestingly, this modification also resulted in a significant effect of speed on ΔF (P < 0.05), yet in this case ΔF decreased as simulation speed increased (Fig. 7B).
The results of these simulations indicated that PICs could have a variable effect on ΔF depending on the specific form of PIC activation. In particular, the results of the latter simulation demonstrated that such PIC activation could lead to a similar ΔF profile as seen with accommodation and adaptation (Fig. 6B). This result could be seen to complicate the interpretation of ΔF as obtained experimentally. Therefore, we next sought to determine whether there was an additional test we could use to better distinguish among the motor neuron properties most likely contributing to ΔF observed experimentally.
Effect of varying contraction duration.
Because accommodation is dependent on rate of rise of excitation whereas adaptation is a time-dependent phenomenon, we carried out simulations for which the total duration of the simulations was varied without changing the rate of rise. It seemed reasonable to expect, therefore, that varying the duration of the contraction in this way might lead to differential effects on ΔF when each of these nonlinear properties were activated.
To vary the total duration of the simulation without changing the rate of rise, we inserted a plateau of constant excitation between the rising and falling phases of excitation. We used the moderate rate of rise (i.e., 5 s to reach 15% maximum force) and tested three simulation durations (Fig. 8A): brief (no plateau, 10-s total simulation time), medium (5-s plateau, 15-s total simulation time), or long (10-s plateau, 20-s total simulation time). For PICs in the simulations shown in Fig. 8, we used the same characteristics as used for the simulations associated with Fig. 7B, namely, a 500-ms activation time and 50% decay rate in 10 s. In addition, we also carried out a set of simulations using the simpler form of PIC with instantaneous onset and constant amplitude.
As shown in Fig. 8B, for the control condition, there was no significant effect (P > 0.05) of contraction duration and ΔF values were small. Likewise, when the simple form of PIC was used, there was no effect of contraction duration on ΔF (P > 0.05, Fig. 8B). However, when the more complex form of PIC was enabled, there was a significant effect of contraction duration (P < 0.05) and post hoc analysis indicated a significant decline in ΔF between the brief and long contractions (Fig. 8B). When accommodation was included in the simulations, ΔF was not different across simulation durations (P > 0.05, Fig. 8B). Finally, when adaptation was activated, ΔF increased significantly (P < 0.001) as simulation duration increased (Fig. 8B). Therefore, by using contraction patterns that vary in duration only, it might be possible to differentiate among the contributions of these properties to ΔF recorded in vivo.
Effects of multiple nonlinearities on ΔF.
Ultimately, motor neurons may operate with all of these nonlinearities activated concurrently. Therefore, to gain insight into how these properties might interact, we ran simulations with all three properties (PICs, accommodation, and adaptation) simultaneously activated. This was done for simulations that varied in contraction speed (Fig. 9A) and contraction duration (Fig. 9B). The excitatory drive functions associated with slow, moderate, and fast speeds (Fig. 9A) were the same as those used in the simulations shown in Fig. 6A. Likewise, the excitatory drive functions for brief, medium, and long contractions (Fig. 9B) were the same as those associated with Fig. 8A. It should be pointed out that the moderate speed condition in Fig. 9A is the same as the brief duration condition in Fig. 9B. In all cases, we used the more complex form of PIC with a 500-ms activation time and 50% decay in 10 s.
The overall magnitudes of ΔF were larger in these simulations compared with previous simulations, indicating a general additive effect of these properties on ΔF. Furthermore, as shown in Fig. 9A, there was a significant effect of contraction speed on ΔF (P < 0.001) with a significant decrement (Holm-Sidak test, P < 0.05) in ΔF for each tested increment in speed. Similarly, there was a significant effect of contraction duration (Fig. 9B) on ΔF (P < 0.001) with a significant increase (Holm-Sidak test, P < 0.05) in ΔF for each increase in duration. In both cases, the pattern of change in ΔF with speed or duration was dissimilar to that occurring when PICs alone were enabled (cf. Fig. 6B and Fig. 8B) but was generally comparable to that arising from simulations in which accommodation or adaptation operated (also, Fig. 6B and Fig. 8B).
Experiments involving reduced preparations in many species have shown motor neurons to exhibit a variety of nonlinear spiking behaviors. In this study, we have focused on three possible contributors to nonlinear spiking responses: PICs, which provide an intrinsic source of depolarizing current to boost firing rate when activated; spike-threshold accommodation, which alters the current needed to recruit a neuron depending on the rate of rise of current; and spike-frequency adaptation, which leads to a decrease in firing rate over time. Because it is experimentally difficult to characterize how these intrinsic properties influence motor neuron activity in awake behaving animals, we therefore opted to address this issue using computer simulations. Furthermore, we restricted our focus in this study to the intriguing ΔF phenomenon, which is seen experimentally (Gorassini et al. 2002a, 2002b, 2004; Kiehn and Eken 1997; Mottram et al. 2009; Oya et al. 2009; Powers et al. 2008; Revill et al. 2010; Revill and Fuglevand 2009; Stephenson and Maluf 2010; Udina et al. 2010) and has been accounted for as primarily arising due to one of these intrinsic properties, namely, to PICs (Gorassini et al. 2002a; Mottram 2009; Oya 2009). Our results indicate that a positive ΔF can arise only by the inclusion of nonlinear intrinsic properties, and as such, previous models of motor unit populations (Fuglevand et al. 1993; Heckman and Binder 1991) do not account for this physiological feature of motor unit behavior. Furthermore, our results suggest that any one of these nonlinear properties can lead to positive ΔF. Nevertheless, by examining the pattern of change in ΔF across contractions of varying duration, we have shown that this approach may provide an experimental means to identify the predominant nonlinear property contributing to ΔF.
ΔF estimation of PICs: confounding effects of adaptation and accommodation.
The original motivation underlying the measurement of ΔF was that it might provide an indicator of the presence and magnitude of PICs influencing firing of human motor units during natural behaviors (Gorassini et al. 1998, 2002a, 2002b; Kiehn and Eken 1997). The ΔF calculation relies on the firing rate of a lower threshold, “control” motor unit to faithfully reflect changes in synaptic drive to the entire motor unit pool. If, however, adaptation affects firing of the motor unit pool, then control-unit firing rate will not necessarily provide an accurate reflection of changes in synaptic drive. Nevertheless, it has been argued that if adaptation magnitudes were similar across motor units, then adaptation would have little influence on ΔF (Gorassini et al. 2002a; Udina et al. 2010). In our simulations, adaptation values were matched to experimentally observed values (Kernell and Monster 1982), and all motor units adapted in the same way, although the absolute magnitude of adaptation varied based on the extent of excitation above threshold for each motor unit. Despite adaptation being implemented in an identical way across the simulated motor unit population, adaptation still had a marked effect on ΔF. Therefore, in the absence of PICs, a test unit can be recruited and derecruited at the same level of synaptic drive but will appear to be derecruited at a lower level of synaptic drive because the control-unit firing rate has diminished due to adaptation. Furthermore, our simulations suggest that this effect should be magnified by increasing the duration of the contraction.
An additional argument has been put forth that adaptation might not be present during voluntary contractions (Gorassini et al. 2002a; Udina et al. 2010) or during fictive locomotion (Brownstone et al. 2010). Evaluating the presence of adaptation in human motor units is difficult, since it is not possible to measure directly the synaptic drive to motor units. Nonetheless, indirect evidence for firing rate adaptation has been suggested in human studies. For example, Moritz et al. (2005) and Barry et al. (2007) demonstrated that the initial firing rate achieved on recruitment is not the lowest level of steady firing motor units can achieve. Indeed, during threshold contractions, the minimal firing rate of human motor units decreased by 2–3 imp/s over a period of 5–20 s. These results were interpreted as being consistent with the effects of firing rate adaptation (Barry et al. 2007). Likewise, Johnson et al. (2004) showed a significant increase in surface EMG amplitude while subjects attempted to hold constant the firing rates of individual motor units. The increase in EMG was thought to reflect an overall increase in synaptic drive to the motor unit pool needed to offset the tendency for firing rate to decline caused by adaptation in the target motor unit. Therefore, it seems probable that firing rate adaptation is playing some role in shaping the activities of motor units during voluntary contractions.
As far as we are aware, accommodation has not been previously considered as a factor that could influence ΔF. Bawa and colleagues, however, have suggested that spike-threshold accommodation might partly underlie rotation in the activities of low-threshold motor units during prolonged contractions (Bawa and Murnaghan 2009; Bawa et al. 2006). Furthermore, accommodation might partially be responsible for the reduction in motor unit recruitment threshold as speed of contraction increases (Budingen and Freund 1976; Desmedt and Godaux 1978). As implemented in this study, accommodation would shift the first spike of the test unit toward a higher threshold, thus leading to a greater level of synaptic drive (as assessed by the firing rate of the control unit) at test unit recruitment compared with test unit derecruitment. This, in turn, led to significant positive values of ΔF when accommodation was enabled. Therefore, this result also suggests that positive ΔF values can arise for reasons other than the presence of PICs.
Although the model used in this study simulated each nonlinear property as distinct and noninteracting, this is unlikely to be the case physiologically (Hamm et al. 2010; Iglesias et al. 2011). For example, sodium channel inactivation may contribute to both accommodation (e.g., Schlue et al. 1974c) and adaptation (e.g., Miles et al. 2005). Therefore, differential effects of these properties on ΔF may not be distinguished as readily as implied by the simulations.
Also, the effects of accommodation were simulated to be restricted to the first spike and to depend only on rate of rise of excitation. It should be recognized, however, that with increasing stimulus strength, membrane potential between spikes may achieve relatively depolarized values (Schwindt and Crill 1982; Turkin et al. 2010) and could lead to sodium channel inactivation. As such, accommodation might affect motor neuron discharge at time points after spiking onset. Whether this happens or not will depend critically on the magnitude and perseverance of afterhyperpolarization (AHP) conductances needed to alleviate sodium channel inactivation over the course of motor neuron discharge (Iglesias et al. 2011).
Another limitation with the present study has to do with the representations of PICs used in the simulations. For example, the absolute amplitude of PICs was assigned the same value for all motor neurons in the pool. Although this formulation is more or less consistent with experimental data obtained in cat motor neurons (Lee and Heckman 1998a), it is at variance with findings in the rat (Hamm et al. 2010). Also, the magnitudes of PICs relative to rheobase for the simulated population were modest: approximately twice rheobase for the lowest threshold neurons and close to zero for the highest threshold neurons. Such low values of assigned PICs might partly account for the relatively small values of ΔF (∼2 imp/s) in simulations during which PICs were enabled. Experimentally, average values of ΔF in human subjects are ∼4 imp/s (Gorassini et al. 2002a; Mottram et al. 2009; Stephenson et al. 2011; Stephenson and Maluf 2010). It should be noted, however, that when multiple nonlinearities were included in the model (Fig. 9), simulated ΔF values did approach 4 imp/s.
In addition, only relatively simple temporal patterns of PICs were simulated, and the same patterns were imposed across the entire motor neuron population. The permutations tested ranged from the simplest case of instantaneous onset and fixed magnitude (Gorassini et al. 2004) to more complex forms of PIC behavior where the PIC needed time to activate fully and then decayed over time (Bennett et al. 2001a, 2001b; Lee and Heckman 1998a; Li and Bennett 2003; Li et al. 2004a). PIC magnitudes, however, may also vary with the amount of current injected (Hultborn et al. 2003) and amount of inhibition present (Bui et al. 2006; Hultborn et al. 2003; Hyngstrom et al. 2008; Kuo et al. 2003). Furthermore, the model did not take into account the diverse nature of PIC activation observed across motor neurons within a pool. For example, PICs appear to decay more markedly in high-threshold compared with low-threshold motor neurons (Lee and Heckman 1998a). Ultimately, addressing each of these limitations will require more mechanistic and spatially complex neural models (e.g., Booth et al. 1997; Bui et al. 2006; Carlin et al. 2000b; ElBasiouny et al. 2006; Manuel et al. 2007; Miles et al. 2005; Powers 1993; Taylor and Enoka 2004). Nevertheless, the model as implemented in this study provided a coarse view as to how various nonlinear properties influence motor unit activity and the development of muscle force.
Other intrinsic properties.
Our consideration of three intrinsic motor neuron properties (PICs, accommodation, and adaptation) should not be considered complete. For example, Wienecke and colleagues (2009) recently showed AHP in motor neurons to elongate over the course of a ramp current injection. The mechanisms underlying AHP elongation were not elucidated, but a functional outcome of this elongation was that firing rate at derecruitment was lower than at recruitment. Similarly, Li and Bennett (2007) and Hamm et al. (2010) provided evidence for a voltage-activated, sustained outward current during triangular voltage-clamp experiments in rat motor neurons. Furthermore, Iglesias et al. (2011) recently demonstrated that slow inactivation of sodium channels can lead to hysteresis in the current-frequency relationship. Because all these phenomena will likely contribute to ΔF, future investigations of ΔF should incorporate these findings.
The results presented in this study ultimately lead to testable predictions for human subject experiments. If motor unit pairs are compared across different contraction speeds and ΔF remains relatively constant, then PICs are likely a primary contributor to ΔF. If, however, ΔF is inversely modulated with contraction speed, then ΔF could be a result of accommodation, adaptation, or complex PIC behavior. To better distinguish among these three properties, varying contraction duration should be useful. If ΔF decreases as contraction duration increases, then PICs that activate slowly and show inactivation are likely to be the dominant property. If ΔF remains stable across contraction durations, PICs of fixed magnitude or accommodation are the most likely explanations. Finally, if ΔF increases as contraction duration increases, then adaptation is probably the dominant property contributing to ΔF. Indeed, recent experimental work by Stephenson and Maluf (2011) showed a progressive increase in the magnitude of ΔF as contraction duration increased up through 20 s (the longest duration tested in the present study). This finding, therefore, is consistent with a primary role of adaptation rather than PICs in mediating ΔF.
In summary, modifications were made to an existing motor unit pool model to test the effect of PICs, accommodation, and adaptation on the output of a motor pool. The simulation results suggest that any of these properties could lead to positive ΔF values. Caution, therefore, should be exercised when interpreting positive ΔF as an exclusive indicator for the presence of PICs. Varying contraction duration and speed could provide experimental means to partially distinguish among these properties as contributors to ΔF.
This study was supported by National Institute of Neurological Disorders and Stroke Grant NS61146 (to A. J. Fuglevand) and ARCS Fellowship (to A. L. Revill).
No conflicts of interest, financial or otherwise, are declared by the author(s).
- Copyright © 2011 the American Physiological Society