Jones, Kelvin E. and Parveen Bawa. Computer simulation of the responses of human motoneurons to composite 1A EPSPS: effects of background firing rate. J. Neurophysiol. 77: 405–420, 1997. Two compartmental models of spinal alpha motoneurons were constructed to explore the relationship between background firing rate and response to an excitatory input. The results of these simulations were compared with previous results obtained from human motoneurons and discussed in relation to the current model for repetitively firing human motoneurons. The morphologies and cable parameters of the models were based on two type-identified cat motoneurons previously reported in the literature. Each model included five voltage-dependent channels that were modeled using Hodgkin-Huxley formalism. These included fast Na+ and K+ channels in the initial segment and fast Na+ and K+ channels as well as a slow K+ channel in the soma compartment. The density and rate factors for the slow K+ channel were varied until the models could reproduce single spike AHP parameters for type-identified motoneurons in the cat. Excitatory synaptic conductances were distributed along the equivalent dendrites with the same density described for la synapses from muscle spindles to type-identified cat motoneurons. Simultaneous activation of all synapses on the dendrite resulted in a large compound excitatory postsynaptic potential (EPSP). Brief depolarizing pulses injected into a compartment of the equivalent dendrite resulted in pulse potentials (PPs), which resembled the compound EPSPs. The effects of compound EPSPs and PPs on firing probability of the two motoneuron models were examined during rhythmic firing. Peristimulus time histograms, constructed between the stimulus and the spikes of the model motoneuron, showed excitatory peaks whose integrated time course approximated the time course of the underlying EPSP or PP as has been shown in cat motoneurons. The excitatory peaks were quantified in terms of response probability, and the relationship between background firing rate and response probability was explored. As in real human motoneurons, the models exhibited an inverse relationship between response probability and background firing rate. The biophysical properties responsible for the relationship between response probability and firing rate included the shapes of the membrane voltage trajectories between spikes and nonlinear changes in PP amplitude during the interspike interval at different firing rates. The results from these simulations suggest that the relationship between response probability and background firing rate is an intrinsic feature of motoneurons. The similarity of the results from the models, which were based on the properties of cat motoneurons, and those from human motoneurons suggests that the biophysical properties governing rhythmic firing in human motoneurons are similar to those of the cat.
The biophysical properties of a motoneuron structure the computational process in which the cell is engaged. These biophysical characteristics include both active and passive parameters resulting from a host of ligand- and voltage-dependent ionic channels, the resistive and capacitive elements of the membrane, and the ionic concentration inside and outside the cell. None of these factors can be assessed in human studies, and for this reason, results from animal studies, especially on cat motoneurons, are cautiously extrapolated to the human spinal cord.
In repetitively firing spinal motoneurons of the cat, intracellularly recorded membrane voltage trajectories between spikes undergo predictable changes with increases in firing rate (Baldissera and Gustafsson 1974a,b; Schwindt and Calvin 1972). Within the primary range of firing, an increase in firing rate is accompanied by a decrease in the maximum depth of the membrane voltage trajectory (the scoop) following a spike with relatively little change in the latter portion of the trajectory (the ramp) in which the membrane potential approaches threshold in a ramplike fashion (Schwindt and Calvin 1972). The changes in membrane voltage at different firing rates are accompanied by firing rate dependent changes in membrane conductance (Schwindt and Calvin 1973). These changes in scoop depth and conductance during the interspike interval (ISI) are expected to alter the “excitability” of the motoneuron and, because of their rate dependence, motoneuronal excitability is predicted to vary with firing rate. In addition to the effects of membrane voltage and conductance, the excitability of the motoneuron is affected by the level of membrane noise produced by background synaptic input (Calvin and Stevens 1968; Matthews 1996).
Recent results on human motoneurons have illustrated that motoneuron excitability, measured as the response to composite excitatory postsynaptic potentials (EPSP) inputs, varies with background firing rate, excitability being higher at slower firing rates (Jones and Bawa 1995; Jones et al. 1995; Olivier et al. 1995). These results contradicted some earlier studies in which motoneuron excitability was reported to be independent of background firing rates (Ashby and Zilm 1982; Miles et al. 1989). The interpretation of the results from the earlier studies led to speculations about the biophysical properties of human motoneurons and the proposal of a model for the repetitive firing of human motoneurons distinct from the behavior of cat motoneurons (Ashby and Zilm 1982). Alternatively, the results showing a relationship between motoneuronal excitability and background firing rate have been interpreted in the context of the biophysical properties of cat motoneurons. These authors have suggested that the greater response to an excitatory input at lower firing rates is due to the curved shape of the membrane voltage trajectory at the end of the interspike interval during slow repetitive firing (Jones and Bawa 1995; Olivier et al. 1995).
The objectives of the following simulations were: to develop motoneuron models based on the biophysical properties of cat motoneurons and use these to simulate the responses of repetitively firing human motoneurons to transient input volleys; and to explore the relationship between excitability and background firing rate in these models and evaluate the importance of various biophysical factors in determining the relationship between background firing rate and the response to composite 1A EPSP input.
To examine whether the results, relative to fast and slow motoneurons, were similar, the motoneuron models in the present study have been fashioned architecturally and electrotonically after a S and a FR type motoneuron innervating the medial gastrocnemius muscle of the cat (Cullheim et al. 1987; Fleshman et al. 1988). The conductances used to simulate active properties included a sodium and fast potassium conductance in the initial segment as well as a sodium, fast potassium, and slow potassium conductance in the soma (Dodge and Cooley 1973; Traub 1977). The slow potassium conductance was described as a time- and voltage-dependent process rather than a calcium-mediated process for simplicity. The results of the simulations using these models support the idea that extrapolation of properties of cat motoneurons to human motoneurons is appropriate for the explanation of previous human results (Jones and Bawa 1995). This work has been presented as an abstract (Jones and Bawa 1996a,b) and forms part of a PhD thesis (Jones 1995).
All simulations were done with the Nodus software (De Schutter 1989) on a Macintosh Centris 660AV computer. Once the characteristics of cat motoneurons were reproduced by the models, experiments previously carried out on human motoneurons were simulated with the model motoneurons.
Motoneurons were modeled in a compartmental fashion with partitions representing the soma, initial segment, and equivalent dendrite. The measurements of the compartments were based on the reconstruction of cat lumbosacral motoneurons reported in the literature (Cullheim et al. 1987; Fleshman et al. 1988). The specific cells used from these publications were a S type soleus motoneuron (identified as cell 35/4) and a FR type MG motoneuron (cell 43/5). The structure of the two motoneuron models is illustrated in Fig. 1. The somata are represented by spheres of diameters 50.9 and 48.8 μm for the S and FR motoneuron, respectively (Cullheim et al. 1987). The initial segment compartments are represented by identical cylinders with a diameter of 10 and length of 100 μm (Traub 1977) in both motoneurons. The dendrites for each motoneuron are represented by an equivalent cable whose diameter changes as a function of distance from the soma (Clements and Redman 1989; Fleshman et al. 1988). Each equivalent dendrite is based on measurements taken from Fleshman et al. (1988) (Fig. 9, Step Model). The diameter of the stem compartment of the equivalent dendrite (deq(stem) ) was estimated from the equation Equation 1where d j is the mean diameter of the jth stem dendrite for n number of stem dendrites (Fleshman et al. 1988). For the S type motoneuron, the equivalent dendrite with a stem diameter of 25 mm, extended for 17 compartments to a final diameter of 0.63 μm. The FR type motoneuron had an equivalent dendritic stem diameter of 40 μm; this gradually narrowed to a diameter of 1 μm in the 19th and final compartment. Both dendrites terminated in a closed end. The length of each dendritic compartment was varied so that no compartment had an electrotonic length >0.2 λ (Segev et al. 1985). The total lengths of the equivalent dendrites were 7,000 and 6,675 μm for the S and FR type motoneurons, respectively.
Passive and active parameters
The passive parameters determining the cable properties of the neurons were set according to the values reported by Fleshman et al. (1988) for these two cells. Specific membrane capacitance, C m, was 1.0 μF/cm2, specific cytoplasmic resistivity, R i, was 70 Ω⋅cm and the specific membrane resistivity, R m, was set according to the step model of membrane resistivity (Clements and Redman 1989; Fleshman et al. 1988; Segev et al. 1990). The S type motoneuron had a value of R m at the soma of 0.7 kΩ⋅cm2, which was stepped to a value of 20 kΩ⋅cm2 for the dendritic compartments. The FR type motoneuron had values for R m of 0.225 and 11 kΩ⋅cm2 for the soma and dendrites, respectively. The resulting electrotonic lengths of the equivalent dendrites were 2.29 and 2.448 λ for the S and FR type neurons, respectively. The R m of the initial segment in both neurons was set equal to that of the soma.
Active properties were associated with the initial segment and soma compartments only and consisted of five voltage-gated ionic currents (equations in appendix) simulated by standard Hodgkin-Huxley (H-H) equations (Hodgkin and Huxley 1952) as described and adapted from Traub (1977). The initial segment contained a fast sodium conductance (Eq. EA5 ) and a fast delayed rectifier type potassium conductance (Eq. EA10 ) (cf. Barrett and Crill 1980). The soma compartment contained a fast sodium conductance (Eq. EA13 ), a fast potassium conductance (Eq. EA18 ), and a slow potassium conductance (Eq. EA21 ), which gave rise to the afterhyperpolarization. The “fast” and “slow” potassium conductances are analogous with those described by Barrett et al. (1980) for cat motoneurons. It is well known that the conductance underlying the afterhyperpolarization (AHP) is unaffected by subthreshold membrane potential displacements (Baldissera and Gustafsson 1974a). Therefore it was important that the slow potassium conductance expressed in H-H formalism was not sensitive to membrane potential displacements in the subthreshold region. This was accomplished by an activation time course that had a steep sigmoidal quality, i.e., the activation parameter, q 2 (Eqs. A22 and A23), remained near 0 for all voltages below threshold then quickly rose to full activation. Thus the activation of the conductance is triggered by a spike and then decays in a strictly time-dependent manner. The time course of this decay was exponential throughout and did not reproduce the S-shaped time course reported by Baldissera and Gustafsson (1974a). The sodium and “fast” potassium conductances in the initial segment were activated at a lower threshold relative to the equivalent conductances in the soma (Dodge and Cooley 1973). All sodium and fast potassium conductances were the same for the two types of motoneurons; however, the slow potassium conductance parameters were adjusted differently in each model as described below. The density of sodium and fast potassium conductance was set according to previous models of cat spinal motoneurons (Traub 1977; Traub and Llinás 1977) and were the same for both models.
For the firing behavior of the motoneuron, the following equation was solved for membrane potential, Vm[x] , at the soma Equation 2where the subscript [x] refers to the soma compartment and [y] to the compartments linked to [x] (i.e., d0 and the initial segment); C m, the total somatic capacitance; E R, the resting membrane potential; g ion[x], the conductance for a given channel in the soma; t, time; E ion, the equilibrium potential; E [y], the membrane potential in a linked compartment; R [xy], the axial resistance between compartments; and I inj, the injected current. The total somatic resistance is given by Equation 3with a diameter, D = 50.9 and 48.8 mm for the S and FR type motoneurons, respectively. R [x] is taken to be the inverse of the passive leak conductance g l in our model. The first term in Eq. 2 represents the leak current, the next term represents the sum of the ionic currents located in the soma compartment. The contribution of currents from other compartments is represented by the third term and the injected current is represented by the final term. Toobtain V m[x] as a function of time, Eq. 2 was integrated using the forward Euler method (Moore and Ramon 1974) with a variable time step according to the equation Equation 4where V 0 is the initial value of membrane voltage and the initial slope. The time step was varied so that the change in membrane voltage in a single time step was never >0.2 mV. With this method, time steps ranged between 0.04 and 2.0 ms in the simulations reported. Similar equations to Eq. 2 were solved for all compartments in the model, taking in to account compartments closed on one end, according to standard compartmental modeling techniques (Segev et al. 1989).
Noisy membrane simulations
During voluntary isometric contractions in which human subjects are asked to maintain a constant firing rate, the motor unit fires, not with a constant ISI as above, but exhibits some level of variability (Clamann 1969; Matthews 1996; Person and Kudina 1972). The source of this variability is the ongoing bombardment of the motoneuron by synaptic inputs resulting in fluctuation of the motoneuron membrane potential (Calvin and Stevens 1968). Variability in the ISIs was incorporated into the simulations to determine whether the behavior of the model motoneurons approximated the human experimental data (Jones and Bawa 1995). This was accomplished by injecting noisy current into various compartments of the motoneuron models. White noise current with a Gaussian distribution of amplitudes was used which had a mean of 0.0 nA and a variable standard deviation. Figure 2 illustrates the effect of injecting the noisy current into various compartments on the membrane voltage of the soma compartment of the models. The top three panels show the membrane voltage computed for the soma compartment when the noisy current is injected in to the soma, first dendritic (d0) and second dendritic (d1) compartments, respectively. In each case, the standard deviation of the noisy current injection was 20 nA. It is clear that some filtering occurs as the current injection is moved from the soma along the dendrite. Injection of noisy current into d1 resulted in membrane potential fluctuations in the soma compartment that were similar to those reported for cat motoneurons during muscle stretch (Calvin 1966; Calvin and Stevens 1968; Cope et al. 1987; Gustafsson and McCrea 1984).
The statistical properties of the membrane voltage fluctuations were analyzed and compared with results from cat motoneurons. The fourth panel in Fig. 2 shows the amplitude distribution histogram of membrane voltage in the soma compartment when noisy current (with0 ± 20 nA; mean ± SD) was injected into d1. The histograms illustrate the probability of finding the membrane potential at a particular voltage during the 1-s sampling period as well as the peak-to-peak amplitude of membrane noise. The mean of the distribution is slightly negative, which reflects the true resting potential of the model neurons being slightly <0 mV. The shapes of these histograms are similar to those reported for cat motoneurons under conditions of increased synaptic noise (Calvin 1966). To evaluate the rate of change of membrane potential fluctuations due to the noisy current injection, the autocorrelation, R xx(τ), of the membrane potential was calculated from a 1-s sample of data. The results are shown in the bottom panel of Fig. 2 for current injection into the soma, d0 and d1 (thick line) compartments. Even when the current was injected in the soma, the membrane voltage shows an exponential shape indicating filtering due to the RC properties of the soma. However, this correlation fell off more rapidly than that expected from background synaptic input to a motoneuron (Calvin 1966). This was remedied by injection of the noisy current into the dendritic compartments, which showed autocorrelations typical for membrane potential fluctuations in cat motoneurons (Calvin 1966). The conclusion from this analysis was that the noisy currents, with the characteristics as outlined above, should be injected into the d1 compartment to match characteristics of noise in the soma similar to those reported for cat motoneurons.
The model motoneurons were made to fire rhythmically by constant current injection at physiological rates observed in human motoneurons (range 5–20 impulses/s) whereas synaptic EPSPs or pulse potentials (PPs) perturbed the firing of the motoneuron. The protocol used during the simulations was similar to that used during human experiments (Jones and Bawa 1995). Stimuli, either a current pulse input or synchronous synaptic input, were given randomly with respect to motoneuron spikes with an interstimulus interval between 0.5 and 2 s. Model data, equivalent to 4 min of human data, were simulated in as much as 30 h of computation depending on the model, type of input (synaptic or current pulse), and the presence of noisy current injection.
Spike train data resulting from the simulations were analyzed using the SPIKE2 software (Cambridge Electronic Design) on a 486 personal computer as described by Jones and Bawa (1995). In short, peristimulus time histograms (PSTHs) were constructed between the stimulus (current pulse giving rise to PPs or EPSPs) and the spikes of the model motoneuron. The onset and offset of an excitatory peak in a PSTH were identified as the first and last bins of a consecutive series of bins in which the number of counts was greater than the mean baseline counts per bin plus 2 SD. The effect of the input on the firing probability of the motoneuron was quantified in terms of response probability, P = the area of the peak above the background over the number of stimuli (Bawa and Lemon 1993). Background firing probability was determined from the average firing probability in a 50-ms window preceding the stimulus at time 0 in the PSTH. The cumulative summation (CUSUM) was calculated by summating the bins in the PSTH after subtraction of background firing probability.
The results will be divided into three sections. The first section will describe the properties of the motoneuron models. The purpose of this section is to compare the results from the models with the data describing the behavior of cat motoneurons on which the models are based. The second section will describe the results when experimental paradigms carried out on human motoneurons were simulated using noisy model motoneurons. The purpose of this section is to determine whether the models show a similar relationship between background firing rate and response probability as that seen for human motoneurons. The third and final section, explores the biophysical factors contributing to the relationship between response probability and background firing rate.
Properties of the motoneuron models
To obtain the input resistance, R N, of the two motoneuron models, hyperpolarizing current pulses of 1 nA amplitude and 50 ms duration were applied to the soma compartment (Fig. 3 A). The linear steady state current/voltage relationships (Fig. 3 B) resulted in R N = 3.9 and 2.0 MΩ for the S and FR type motoneurons, respectively. The time constants, τ0, for the two motoneurons were obtained by injecting 20-nA, 0.2-ms hyperpolarizing pulses into the soma compartments. The resulting voltage transient was plotted on a semilogarithmic plot and τ0 was computed by curve fitting with the equation Equation 5(Fleshman et al. 1988). The equation was fit to the transient record during a 5-ms period centered at 1.5 times the experimental value of τ0 reported in Fleshman et al. (1988). The τ0 values computed in this manner were 16.4 and 8.8 ms for the S and FR type motoneurons, respectively.
A single action potential was elicited by injecting a 0.5-ms, 10-nA depolarizing current pulse into the initial segment and the resulting AHP recorded from the soma is plotted in Fig. 3 C for each motoneuron. The parameters describing the AHP of the two motoneurons are reported in Table 1. The depth and duration of the AHPs were adjusted to match the values reported by Zengel et al. (1985) for type identified motoneurons in the cat by adjusting parameters of the “slow” potassium conductance. Depth of the AHP was adjusted by changing the density of the slow potassium conductance on the soma, and duration was adjusted by changing the numerator in the backward rate function Bq (Eqs. A23 and A26).
After adjustment of the slow K+ conductance in the two models, repetitive firing was elicited by injecting long-duration current steps into the soma compartments of the two motoneuron models. The interspike intervals (ISIs) were measured after adaptation and the inverse, firing rate (f), was plotted with respect to the amplitude of current injected into the soma to produce the f/I curves plotted in Fig. 3 D. The S type motoneuron exhibited both a primary and secondary range (Kernell 1965a) with slopes of 4.4 and 12.3imp⋅s−1⋅nA−1, respectively. The FR type motoneuron was characterized by a single linear regression line with a slope of 5.4 imp⋅s−1⋅nA−1. The f/I curves in the primary range agree with the data reported by Kernell (1965b) where small (i.e., S type) motoneurons are shown to have a less steep slope and begin repetitive firing at lower values of injected current. The secondary range in the S type motoneuron was unexpected but was not explored further as this range of firing rates is seldom reported in normal human motoneurons. The motoneurons also show an inverse relationship between AHP duration and minimal rhythmic firing rate characteristic of cat motoneurons studied intracellularly (Kernell 1965b).
MEMBRANE VOLTAGE TRAJECTORIES DURING RHYTHMIC FIRING AT DIFFERENT FIRING RATES.
In cat motoneurons, an increase in firing rate is accompanied by a decrease in the depth of the membrane voltage trajectories, i.e., the scoop, between spikes (Baldissera and Gustafsson 1974a,b; Schwindt and Calvin 1972). To examine whether the models showed similar changes in membrane voltage trajectories, each of the models was fired at a number of different firing rates and the trajectories are illustrated in Fig. 4. The figure shows changes in the membrane potential trajectories, which are similar to those reported in cat motoneurons, namely the maximum depth of the AHP decreases as firing rates increase (Baldissera and Gustafsson 1974a,b; Schwindt and Calvin 1972). In the S motoneuron, the scoop of the AHP decreased by 2 mV between the fastest and the slowest rates shown and in the FR motoneuron, the scoop decreased by 1 mV. Figure 4 also reveals changes in the latter part of the AHPs as the firing rates increased. Initially at the minimal rhythmic firing rate, 5.6 imp/s for the S motoneuron and 10.2 imp/s for the FR motoneuron, the recovery of the membrane potential toward threshold showed an exponential region. At the faster firing rates, this exponential region gradually decreased and the membrane potential returned to threshold in a more ramplike fashion.
SYNAPTIC AND CURRENT IMPULSE INPUTS TO THE MOTONEURONS.
Synaptic conductances were added to the equivalent dendrites of the two motoneurons according to Segev et al. (1990; Fig. 4, Step Model). Instead of assigning 300 individual synapses of the same peak conductance to the motoneuron models (Segev et al. 1990), the number of synapses in a particular compartment was determined and mod eled as a single equivalent synapse with a peak conductance, g max, equal to Equation 6where g peak(ss) is the peak conductance of a single synapse (5 nS) (Finkel and Redman 1983) and n is the total number of synapses assigned to that compartment. Each equivalent synapse was modeled as a synaptic conductance transient defined by an alpha function with a time to peak conductance of 0.2 ms (Eqs. A27 and A28). The distribution of synaptic input on the equivalent dendrites is shown in Fig. 1. Both motoneurons were assigned 300 synapses; however, these were represented by 12 or 10 equivalent synaptic conductances in the S and FR type motoneurons, respectively. The synaptic conductances were activated synchronously to produce a composite EPSP that was larger and slower in the S type motoneuron (Fig. 4 B). To prevent the motoneurons from firing an action potential in response to an EPSP at resting potential, the value of g peak(ss) was reduced to 1 nS for each synapse. This had the advantage of maintaining the same relative distribution of synaptic input over the equivalent dendrite.
As a means of speeding up the simulations, it was explored whether current pulses of 1.0 ms duration could be injected into the dendrite to produce PPs in the soma that resembled the EPSPs produced by synchronous synaptic activation. This was pursued based on the successful results of Walmsley and Stuklis (1989). Current pulses were injected into sequential dendritic compartments starting with the compartment d0 connected to the soma. The criteria used was to match the profiles of the PPs to the amplitude and rise times of the corresponding synaptically produced EPSPs. These are the two critical factors that affect changes in motoneuron firing probability (Fetz and Gustafsson 1983). A qualitative match was found when current impulse inputs were given to the third dendritic compartment, d2, in both models (Fig. 4 B). These EPSPs and PPs were used for all the simulation results reported.
COMPARISON BETWEEN EPSPS AND PPS AND THEIR RELATIONSHIP TO CHANGES IN FIRING PROBABILITY.
In these simulations, both model motoneurons were made to fire repetitively at rates of ∼10 imp/s. The synaptic inputs, arriving randomly with respect to the motoneuron spikes, increased the firing probability of the rhythmically firing motoneurons as shown by the excitatory peaks in PSTHs in the top panels of Fig. 5. The bottom panels illustrate the CUSUMs of the PSTHs around the peak superimposed on the respective composite EPSPs measured at rest. As has been shown for cat motoneurons (Fetz and Gustafsson 1983), there was a good relation between the shape of the CUSUM and the rising phase of the EPSP.
Current pulses of 9 and 10 nA injected into the third dendritic compartment of the S and FR motoneuron models, respectively, resulted in PPs that qualitatively resembled the EPSPs (Fig. 4 B). Although these PPs did not exactly match the EPSPs, they altered the firing probability of the model motoneurons during rhythmic firing in a similar fashion to the EPSP input. This is illustrated in Fig. 6 for simulations with PP input during repetitive firing at a rate of ∼10 imp/s for both models. The shape of the CUSUMs closely resembled the shape of the rising phase of the PPs.
Comparison of the responses with EPSPs and PPs show that the P values are quite similar with the two inputs in both the S and FR motoneurons. Differences between the two inputs are more evident upon comparing peak widths. For the S motoneuron, the peak width resulting from EPSP input was 0.4 ms longer than that produced by PP input. For the FR motoneuron, this difference was 0.7 ms. These differences are likely due to differences in the slope of the rising phase of the input potentials, which are lower for the PPs. Also the strict definition of peak onset and offset excluded bins at the peak offset resulting from PP input; these bins may have been included in a less strict definition. To a first approximation, the PP input produced similar changes in response probability as that produced by EPSP input, and there was good agreement between the shape of the CUSUM and the rising phase of the PP (Fig. 6, bottom). These results suggested that PPs might be used to simulate changes in firing probability caused by EPSPs.
Simulation of human motoneuron experiments
It seems reasonable that synaptic noise plays a major role in the excitation of human motoneurons (Matthews 1996); therefore it was imperative that the simulations replicated the type of repetitive firing behavior seen in human motoneurons during voluntary activation. Figure 7 illustrates the interspike interval histograms when noisy current with a standard deviation of ±20 nA was added to the current necessary to elicit repetitive firing at different rates (see Methods). For each of the two models, as the mean ISI decreased, i.e., firing rate increased, variability of the ISIs decreased as is typical of the human data. It also was confirmed that the falling phases of the ISI interval histograms were exponential by plotting the surviving intervals transform of the data (Matthews 1996) on a semilog plot. The resulting plots showed that the data points from the falling phases of the ISI interval histograms fell along a straight line, indicating an exponential decay similar to that seen in human motoneurons at various firing rates (see Matthews 1996 for details). However, there were differences between the two neurons. Variability in ISI did not depend on the absolute ISI but rather on its value relative to the single spike AHP (Table 1), i.e., on the intrinsic property of each motoneuron. At the fastest firing rate shown for the two neurons, the coefficient of variation (c.v.) of ISI was 0.08 for both when ISIs were 45% of respective AHP duration (72/160 for S and 41/90 for FR). On the other hand, when ISI was the same for S and FR (72 ms), c.v. was greater for the FR (Fig. 7, middle right) compared with that for the S (Fig. 7, bottom left).
Once it was confirmed that the models could replicate repetitive firing variability similar to that seen in human data, the effects of changes in background firing rate on response probability were investigated. In these simulations, PPs were produced randomly with respect to the motoneuron spikes while the motoneurons were fired with the ISI characteristics shown in Fig. 7. The resulting PSTHs are illustrated in Fig. 8 and show many similarities to those computed for human motoneurons in similar experimental paradigms. For example, the excitatory peaks are followed by a period of decreased firing probability. Unlike the case in the absence of noise, this region is not completely devoid of spikes (cf. PSTHs in Figs. 5 and 6). This feature is typical of PSTHs computed from human motoneurons (Jones and Bawa 1995). The excitatory peaks had an average onset of 1.2 ms and width of 1.8 ms in the S motoneuron model and 1.2 and 1.4 ms in the FR motoneuron model at the three different firing rates. These values for peak onset and width are slightly larger than those reported for the models in the absence of noise input (Figs. 5 and 6) as would be expected. Thus in the presence of noise, the time course of the integral of the peaks, i.e., the CUSUM, becomes longer, diverging from its approximation of the rising phase of the underlying PP. However, the main feature of interest was the response probability calculated from the excitatory peaks. Figure 8 shows that as background firing rates increased, the response probability to a constant PP input decreased. Thus when the motoneuron models were firing at the slowest rates, with the greatest ISI variability, they exhibited the greatest probability of responding to the input. This is the same behavior as has been shown for human motoneurons to both H-reflex and transcranial magnetic stimulation volleys (Jones and Bawa 1995; Olivier et al. 1995).
Factors contributing to the relationship between background firing rate and response probability
That the models exhibited the same disputed relationship between firing rate and response probability as human motoneurons was encouraging and further simulations were done to investigate the factors underlying this relationship. Previously it had been predicted that the trajectory of the membrane voltage between spikes would play a key role in mediating this “rate effect” (Jones and Bawa 1995; Olivier et al. 1995). In the presence of noise, it would be difficult to compute deterministic values of parameters responsible for the rate effect. Therefore the noisy current was removed from these parts of the simulations. However, it first remained to be determined whether the models would exhibit the inverse relationship between response probability and firing rate in the absence of noise.
RELATION BETWEEN BACKGROUND FIRING RATE AND RESPONSE PROBABILITY IN THE ABSENCE OF NOISE.
Each model motoneuron was made to fire repetitively at three different firing rates by injection of current step inputs to the soma compartment. The membrane voltage trajectories at these different rates are shown in Fig. 4 during steady-state firing. Unlike simulations done in the presence of noise, the interspike intervals showed almost no variability and the only variability present was due to effects of the stimulus. During rhythmic firing, current impulses were randomly applied to the motoneurons to generate PPs and the resulting PSTHs were computed. The response probability calculated from the resulting excitatory peaks is plotted with respect to background firing rate in Fig. 9. It is clear from the figure that in the absence of noise, the inverse relationship between response probability and firing rate persists. The data from simulations done in the presence of noise are also shown in Fig. 9 for comparison. The figure shows that in the presence of noise at a given firing rate, the response probability of the model motoneurons to the same input is decreased. This effect seems to be diminished at the fastest firing rate illustrated for the FR model motoneuron, 24 imp/s. At this rate, response probability computed in the absence of noise is only slightly greater than that computed in the presence of noise.
NOISELESS MEMBRANE: FACTORS CONTRIBUTING TO THE EFFECT OF FIRING RATE ON RESPONSE PROBABILITY.
Having established that the feature of interest in the present study was evident in simulations done in the absence of noise, the factors giving rise to this relationship were sought. In the absence of noise, the response probability of a motoneuron to an impulse input during the ISI is a step function. After the spike, and below a certain postspike delay, d e, the response probability is 0. For postspike delays greater than or equal to d e, the response probability is unity. We have called d e the effective delay, and a number of parameters were evaluated with respect to this delay (Table 2). The computations showed that as the firing rate increased the effective delay decreased. That is, the same current pulse input was able to trigger a spike at a shorter value of d e with increasing firing rate of the motoneuron. The effective delay was then usedto calculate the percentage of the ISI in which the PP was within reach of threshold As firing rate increased the percentage of the ISI within reach of threshold to the PP decreased (Table 2). Therefore, although d e decreased at the faster firing rates, it did not decrease in proportion to the reduction in ISI at the faster firing rates. Regression analysis of the relationship between the percentage of the ISI within reach of threshold (%ISI) and response probability suggested that changes in %ISI explain, in large part, the relation between firing rate and response probability (r 2 values were 0.98 for the S and 0.82 for the FR motoneuron).
The next question that was addressed was whether action potential threshold changed with firing rate and whether the firing threshold varied during the ISI (Calvin 1974; Reyes and Fetz 1993). The membrane potential at which the action potential began its rapid rise was evaluated from the derivative of the voltage record at the different firing rates, and defined as the firing level (Reyes and Fetz 1993). Simulation data showed that the firing level over the range of firing rates tested did not vary much and had a mean value of10.0 ± 0.04 mV (mean ± SD) in the S motoneuron and 11.3 ± 0.10 mV in the FR motoneuron. Therefore the firing level did not vary as a function of firing rate, yet it was still not clear whether the firing threshold varied during the ISI.
As an alternative means of estimating voltage threshold to dynamic inputs, a 1.0-ms current pulse was injected into the soma compartment under resting conditions. The amplitude of the current pulse was adjusted so that it was just able to trigger an action potential, and the voltage at which the action potential began its rapid rise was considered as threshold. These threshold values tested under dynamic conditions were similar to the values for firing level evaluated during repetitive firing: 10.4 mV for the S and 11.3 mV for the FR motoneuron model.
To begin to address the issue of firing threshold changes during the ISI, the membrane potential at the effective delay was measured at the different firing rates to estimate the distance between membrane voltage and the firing level (Table 2). In the S motoneuron, the average membrane potential at d e was 3.4 ± 0.25 mV. However, there was a tendency for the membrane potential at the effective delay to increase with increasing firing rates. In the FR motoneuron, the average membrane potential was 4.9 ± 0.06 mV at d e. When the peak amplitude of the PP measured at rest (3.4 and 2.5 mV for the S and FR motoneuron, respectively) is added to the membrane potential at the effective delay, the sum falls short of the firing level by 3.2 mV in the S motoneuron and 3.9 mV in the FR motoneuron. This would imply that either the PP amplitude was significantly greater during rhythmic firing or that the threshold for an action potential decreased during the ISI.
To investigate the first option, i.e., whether the amplitude of the PP elicited by a constant current pulse input varied during rhythmic firing, the size of the PP was measured at sequential times during the ISI at different firing rates. This was done for both model motoneurons with similar results and the results from the S motoneuron are illustrated in Fig. 10. The motoneuron was fired repetitively at 10.3 imp/s. and a PP was induced at sequential 10-ms delays beginning with a delay of 10 ms after the spike (Fig. 10, top). At a postspike delay of 10 ms, the amplitude of the PP was reduced with respect to that at rest and progressively increased in size at longer delays after the motoneuron spike until a postspike delay of 50 ms at which the amplitude of the PP had returned to resting values. At a delay of 60 ms, the PP was of sufficient magnitude to bring the motoneuron to threshold eliciting an action potential. This same pattern of increasing amplitude of the PP was found in the FR motoneuron during rhythmic firing. Thus an affirmative answer to the option posed above could be given, the PP amplitude varied as a function of time following the spike. However, it was not yet clear whether this variation in PP amplitude was affected by firing rate.
To pursue an answer to this question, the amplitude of the PP at different postspike delays was evaluated during rhythmic firing at a number of different firing rates using the same current pulse input to the S motoneuron. The data are plotted in Fig. 10 (bottom), where the ordinate scale is a measure of the PP amplitude with respect to that obtained at rest. The figure shows that PP amplitude is approximately equal for all firing rates at the shortest delay evaluated (10 ms) and reduced considerably as compared to that at rest. As the postspike delay increased, the rate of change of PP amplitude varied with firing rate in such a way that the faster the firing rate of the motoneuron, the earlier the PP amplitude returned to resting values. The final values for relative PP amplitude were evaluated at a postspike delay 1 ms before the effective delay (the absolute values for both the S and FR motoneurons are given in Table 2). The slower the firing rate, the greater the relative amplitude of the PP prior to the effective delay. The results from the FR motoneuron were similar to those from the S motoneuron. That is, at a postspike delay of 10 ms, the amplitude of the PP with respect to rest was similar at all rates tested. However the decrease in amplitude at a postspike delay of 10 ms was significantly greater for the S motoneurons, 0.62 ± 0.02, than for the FR motoneuron, 0.78 ± 0.02 (t-test, P < 0.001).
The shapes of the PPs during the ISI also were compared with that evoked at rest. PPs evoked at earlier delays during the ISI not only had smaller amplitudes but also slightly shorter rise times compared with the PP evoked under resting conditions. The average PP evoked during the ISI was calculated from the PPs elicited at 10-ms intervals (as illustrated in Fig. 10, top). The rise time of the average PP was the same as the PP evoked under resting conditions. However, the amplitude of the average PP was dependent on the background firing rate. At the slowest firing rates in the S and FR motoneuron models, the amplitude of the average PP during repetitive firing was the same as that evoked at rest. At faster firing rates the amplitude of the average PP was less than that evoked at rest. This relationship is obvious upon closer examination of Fig. 10, bottom. If the relative amplitudes of the PPs at a given rate are averaged, the average amplitude is greater the slower the firing rate.
Three points may be made here: first, at the same postspike delay, PP is larger at the faster firing rate; second, at d e the PP is larger at the slower firing rate; third, the average amplitude of the PP during the ISI is larger at the slower firing rate. These observations are useful for interpreting the results from human motoneurons. It remained to be determined whether the variation in PP amplitude during the ISI could account completely for the relative distance between the membrane voltage and the firing level. Adding the amplitude of the PP immediately before the effective delay to the membrane potential at the effective delay gave an average value of 7.8 ± 0.3 mV for the S motoneuron and 8.2 ± 0.3 mV for the FR motoneuron. These values are still, on the average, >2–3 mV shy of the control firing levels of 10.0 and 11.3 mV for the S and FR motoneurons, respectively. This result implies that the firing level varied with time during the ISI. Thus in answer to our original query, it appears that in the models both changes in PP amplitude and decreases in spike threshold during the ISI contribute to bridging the gap between membrane voltage and firing level.
The primary aim of these simulations was to construct motoneuron models based on the properties of cat motoneurons and then simulate experiments done on human motoneurons to test whether the results from the models and human motoneurons were similar. If so, then the biophysical parameters that account for the observed results could be identified in the models and implicated in the behavior of human motoneurons. In general, the results were encouraging.
Validity of the models
The first question that must be asked, and answered, is whether the models realistically simulated the properties of cat motoneurons upon which they were based? The underlying philosophy behind these models was that they should be simple and yet retain some of the essential input/output characteristics of motoneurons (Segev 1992). The models succeeded in doing this by incorporating some of the morphological, electrotonic, synaptic, and single spike AHP features of adult cat spinal motoneurons (Cullheim et al. 1987; Fleshman et al. 1988; Segev et al. 1990; Zengel et al. 1985). After the tuning of the models to reproduce single spike AHPs, the models were able to simulate f/I curves with steady-state characteristics similar to cat motoneurons. Furthermore, the relative differences between the S and the FR model in minimal rhythmic firing rates, thresholds to current step injection and slopes of the f/I curves were appropriate relative to S and FR-type cat motoneurons (Fig. 3) (Kernell 1965b). The trajectories of the AHPs during repetitive firing showed similar changes to those described for cat motoneurons with increasing firing rates (Fig. 4). That is, increased firing rates were accompanied by a decrease in the depth of the AHP (Baldissera and Gustafsson 1974a,b; Schwindt and Calvin 1972). The models also demonstrated variation in spike threshold during the ISI that has been reported in cat motoneurons (Calvin 1974) and also has received previous modeling attention (Powers 1993).
RELATION BETWEEN SHAPES OF EPSPS AND CHANGES IN FIRING PROBABILITY.
The models successfully simulated a key feature of cat spinal motoneurons that has been addressed in many previous studies, i.e., the matching of the rising phase of the EPSP to the CUSUM of the corresponding peak of the PSTH. Earlier simulations using a single compartment model, similar to that reported by Powers (1993), failed to reproduce the relationship between the excitatory peak in the PSTH and the underlying EPSP reported for cat motoneurons (Cope et al. 1987; Fetz and Gustafsson 1983; Gustafsson and McCrea 1994). The current models, with multiple compartments and Na+ channels, successfully simulated the relationship between the shapes of the PSTH and the corresponding composite EPSP (Figs. 5 and 6).
ISI VARIABILITY AT DIFFERENT FIRING RATES IN THE PRESENCE OF NOISE.
Another important feature of the models was the successful simulation of ISI variability. The injection of white noise current into the dendritic compartment, d1, of the motoneuron models resulted in variation of the ISIs to a constant depolarizing step current injected into the soma (Fig. 7). The relative variations in the ISIs at different firing rates, measured by the coefficient of variation, are compared for the models and human soleus and FCR motoneurons in Table 3. The changes in the coefficient of variation with firing rate were similar for human soleus motoneurons and the S-type model and for FCR motoneurons and the FR-type model. However, at any one firing rate, the coefficients of variation were notably different between the slow (model S or human soleus) and the fast (model FR or human FCR) motoneurons.
Similar differences in coefficients of variation between fast and slow human motoneurons have recently been reexamined by Matthews (1996). For any one motoneuron, the difference between the coefficients of variation at different firing rates is most likely due to differences in the exponential regions of the AHP trajectories reflected by their minimal rhythmic firing (Fig. 3). This exponential region with low input conductance (Fig. 8) is very sensitive to fluctuations in membrane potential. In the model motoneurons, as firing rates increased, the exponential region decreased making the cells less sensitive to membrane noise, leading to less variability in the ISIs (Matthews 1996).
It also should be pointed out that a lower minimal rhythmic firing rate could be achieved in both models in the presence of noise. In the S model, the minimal rhythmic firing rate in the absence of noise had an ISI of 180 ms, which was increased to a mean of 187 ± 50 ms in the presence of noise. In the presence of noise, the amount of current injected into the soma was 0.45 nA less than that used to elicit minimal rhythmic firing in the absence of noise. It is clear that the addition of noisy current did not increase the mean level of driving current in the soma. This means that the S model was producing rhythmic firing with a mean level of driving current in the soma compartment below that necessary to elicit repetitive firing in the absence of noise. Thus the repetitive firing behavior was brought about by threshold crossings that occurred as a result of the synaptic noise. This would fit with Matthews (1996) hypothesis that during voluntary activation of human motoneurons at slow firing rates, the processes mediating the AHP have been completed before the next spike and that the membrane potential remains below threshold. This places a stronger emphasis on the role of synaptic noise in the production of spikes at these firing rates.
Limitations of the models
While the models simulated many features of motoneurons that were of interest, there were other details of motoneuron response properties that were not incorporated in the models. The main point of divergence between the active properties of the models and real motoneurons lies in the simulation of the conductances responsible for the AHP. This process is likely mediated by a Ca2+-activated K+ channel(s) in motoneurons (Krnjevic et al. 1978; Zhang and Krnjevic 1987) as in many other neurons (Sah 1996). However, there are very little experimental data available to realistically model a calcium-dependent K+ current and the associated Ca2+ currents in adult cat spinal motoneurons. Previous modeling studies have been able to reproduce many of the characteristics of motoneuron firing by using simple time-dependent (Baldissera and Gustafsson 1974a,b) or voltage- and time-dependent (Powers 1993; Traub 1977) slow K+ conductances similar to that used in the present simulations.
Also, the models failed to reproduce in full detail the dynamic discharge behavior of motoneurons. While there was some degree of adaptation present in the models to step current inputs, the first interval f/I relations in the models were lower than those obtained in real motoneurons (Kernell 1965a,b). This may, in part, be due to the higher than average slope of the steady-state f/I relation in the models (Fig. 2) (cf. Kernell 1965b). As pointed out by Powers (1993), a model's early time course of the AHP trajectory may play a role in the discharge characteristics of motoneurons. This early phase of the AHP in the present models is controlled primarily by the outward currents generated by the fast and slow K+ channels. This phase of the AHP in real cat motoneurons depends on the time course of the outward current (Baldisserra and Gustafsson 1974a,b) but also may depend on inward currents producing the delayed depolarization (Nelson and Burke 1967) which has been implicated in doublet discharges (Calvin 1974). We did not attempt to model these parameters of motoneuron behavior in the present simulations as they did not seem critical to the questions addressed; this certainly is an area for improvement in future models.
Simulation of human results
Once the models were able to simulate many of the features of cat motoneurons, they were used to simulate experiments done on human motoneurons. It was important to determine whether the results obtained in the human motoneuron studies could be confidently attributed to the properties of the motoneuron or, if not, whether a different motoneuron model (Ashby and Zilm 1982), or some other spinal mechanisms such as presynaptic inhibition or Renshaw inhibition would help explain the effects seen with voluntary changes in background firing rate (Jones and Bawa 1995).
INVERSE RELATION BETWEEN BACKGROUND FIRING RATE AND RESPONSE PROBABILITY.
When the stimulus arrived randomly with respect to the motoneuron spikes, the motoneuron models showed an inverse relationship between background firing rate and response probability to a randomly elicited transient input. This effect of firing rate on response probability, or “rate effect,” was present in the models in the presence and absence of noise (Fig. 9).
In the noiseless motoneuron, the important factors that determine the magnitude of the rate effect are the depth of the AHP, the duration of the exponential region, and the membrane conductance. With the same input, the PP increased monotonically during the ISI starting with a relatively small value near the spike. At the same postspike delay, the PP was smaller at the slower firing rate (Fig. 10). The smaller PP added to the lower membrane potential at a specific postspike delay would result in a lower response probability at the slower firing rate. However, during random stimulation, the motoneuron is responsive for a larger percentage of the ISI when it is firing at a slow rate than when it is firing at a faster rate, leading to a higher response probability at the slow rate. This higher response probability at slow firing rates may be attributed to the exponential nature of the membrane voltage trajectory and low values of conductance in the latter part of the ISI. As a result of the membrane conductance changes, the amplitude of the average PP is greater at the slower firing rates contributing to a higher overall response probability during random stimulation. We submit that this reasoning on the origin of the rate effect also holds for simulations in the presence of noise.
The lack of rate effect reported by Ashby and Zilm (1982) for tibialis anterior motoneurons may be due to the firing rates used relative to the minimal rhythmic firing rates of these motoneurons or, in addition, due to use of very weak peripheral stimuli resulting in low values of response probability. When higher P values were used with corticospinal volleys, results from the same laboratory (Brouwer et al. 1989) reported the inverse relationship seen in the present modeling and previous human results (Jones and Bawa 1995; Olivier et al. 1995).
Another difference between the human results reporting lack of rate effect (Ashby and Zilm 1982; Miles et al. 1989) and results from our laboratory (Jones and Bawa 1995; Olivier et al. 1995; and present results) is the difference in quantification of response probability. In the original work by Ashby and Zilm (1982), differences in background firing probability levels of the PSTHs during fast and slow firing rates were not accounted for when quantifying the short latency excitatory peak. Similarly, the normalization procedure used by Miles and coworkers (1989) is questionable and may actually mask the rate effect (see Jones and Bawa 1995 for discussion). Similar studies done on cat motoneurons (Fetz and Gustafsson 1983) would suggest that the relative size of the peak above baseline should be further normalized for background firing rate. That is, response probability as calculated in the present study should be normalized by dividing by firing rate. This was not done in the present experiments to allow for a comparison to previously published human motoneuron data (Jones and Bawa 1995). However, it should be said that such normalization would enhance the inverse relationship between firing rate and response probability seen in the human data and the present results.
RESULTS OF TRIGGERED STIMULATION.
As mentioned above, when stimuli were applied randomly with respect to motoneuron spikes, Jones and Bawa (1995) and Olivier et al. (1995) reported the inverse relationship between response probability and firing rate with composite 1a EPSPs and complex corticospinal EPSPs, respectively. However, when stimuli were applied at a fixed time with respect to the motoneuron spike, response probability was higher at the faster firing rate. Two factors contribute to this increase in firing probability with firing rate at a given postspike delay. First, at a fixed time following the spike, the membrane voltage trajectory of the AHP is always lower at the slower rate; this would lead to a lower P value. Second, at a fixed time following the spike, the same synaptic current would produce a smaller EPSP at the slower rate (Fig. 8). When comparing the response probability of the motoneuron at the same postspike time at two different firing rate, the two factors, the depth of AHP and EPSP size, both favor a higher response probability at the faster firing rate.
These simulations support the idea that the repetitive firing behavior of human motoneurons may be represented by biophysical processes that are similar to those that underlie repetitive firing in cat motoneurons. Thus the need for a model of human motoneurons that is different from cat motoneurons is supplanted. The results of this study also have implications for the estimation of EPSP size in human studies (Palmer and Ashby 1992). This usually is done by measuring the excitatory peak of a PSTH. The amplitude of this peak is not only dependent on the size of the EPSP, but also will be affected by the background firing rate of the motoneuron and perhaps the type of motoneuron being recorded. For example, in the present study, a smaller EPSP was generated in the FR motoneuron by a given number of synapses compared with that produced in the S motoneuron. Yet, at a similar firing rate, the response probability was greater in the FR motoneuron. This would imply that to compare responses to a given input in different motoneurons, one would have to account for the underlying post-spike afterhyperpolarization (Matthews 1996).
That the rate effect previously described for human motoneurons is likely to be mediated by the intrinsic properties of the motoneuron is supported by these simulations. However, the contribution of spinal interneuronal circuits can not be dismissed in considering the human results. During voluntary contractions as firing rate is increased, there may be concomitant changes in descending and peripherals volleys along with changes in the excitability of the spinal cord.
The authors are grateful to Å. Vallbo and B. Gustafsson for helpful comments and discussions about this work.
This work was supported by grants from the British Columbia Health Research Council (BCHRF) and the Natural Sciences and Engineering Research Council of Canada (NSERC) to P. Bawa and a NSERC postgraduate scholarship to K. E. Jones.
Address reprint requests to K. E. Jones.
- Copyright © 1997 the American Physiological Society
In Eq. 2, the second term represents the sum of the ionic currents in the soma compartment. These currents can be written as Equation A1where ḡ ion are the maximum conductances for the sodium, fast potassium, and slow potassium channels; m, h, n, and q are time- and voltage-dependent parameters (represented as η in the following 3 equations), which determine the magnitude and time courseof the activation and inactivation of the respective currents. The time and voltage dependence is given by Equation A2where Equation A3 Equation A4The time course and magnitudes of the activation/inactivation parameters used in the simulations of Eq. 2 are given below.
Initial segment currents
Na+ current Equation A5 Equation A6 Equation A7 Equation A8 Equation A9Fast K+ current Equation A10 Equation A11 Equation A12
Na+ current Equation A13 Equation A14 Equation A15 Equation A16 Equation A17Fast K+ current Equation A18 Equation A19 Equation A20Slow K+ current: S type motoneuron Equation A21 Equation A22 Equation A23Slow K+ current: FR type motoneuron Equation A24 Equation A25 Equation A26A larger βq for the FR motoneuron makes τq shorter, resulting in a shorter duration AHP.
Dendritic synaptic current
Equation A28 Equation A27