Recruitment in retractor bulbi muscle during eyeblink conditioning: EMG analysis and common-drive model.
Free full text in Europe PMC
Abstract
Free full text
Recruitment in Retractor Bulbi Muscle During Eyeblink Conditioning: EMG Analysis and Common-Drive Model
1Abstract
To analyze properly the role of the cerebellum in classical conditioning of the eyeblink and nictitating membrane (NM) response, the control of conditioned response dynamics must be better understood. Previous studies have suggested that the control signal is linearly related to the CR as a result of recruitment within the accessory abducens motoneuron pool, which acts to linearize retractor bulbi muscle and NM response mechanics. Here we investigate possible recruitment mechanisms. Data came from simultaneous recordings of NM position and multiunit electromyographic (EMG) activity from the retractor bulbi muscle of rabbits during eyeblink conditioning, in which tone and periocular shock act as conditional and unconditional stimuli, respectively. Action potentials (spikes) were extracted and classified by amplitude. Firing rates of spikes with different amplitudes were analyzed with respect to NM response temporal profiles and total EMG spike firing rate. Four main regularities were revealed and quantified: 1) spike amplitude increased with response amplitude; 2) smaller spikes always appeared before larger spikes; 3) subsequent firing rates covaried for spikes of different amplitude, with smaller spikes always firing at higher rates than larger ones; and 4) firing-rate profiles were approximately Gaussian for all amplitudes. These regularities suggest that recruitment does take place in the retractor bulbi muscle during conditioned NM responses and that all motoneurons receive the same command signal (common-drive hypothesis). To test this hypothesis, a model of the motoneuron pool was constructed in which motoneurons had a range of intrinsic thresholds distributed exponentially, with threshold linearly related to EMG spike amplitude. Each neuron received the same input signal as required by the common-drive assumption. This simple model reproduced the main features of the data, suggesting that conditioned NM responses are controlled by a common-drive mechanism that enables simple commands to determine response topography in a linear fashion.
INTRODUCTION
Classical conditioning of the rabbit's nictitating membrane (NM) response has been extensively studied in two contexts: first, as an exemplar of associative learning (e.g., Dudeney et al. 2007; Gormezano et al. 1983) and, more recently, as a suitable task for investigating cerebellar function in the production of a learned movement (e.g., Hesslow and Yeo 2002; Thompson 2005). Although a major attraction of using the NM response has been its apparent simplicity, accumulating evidence suggests that this simplicity is only relative. For example, although the NM response does indeed have much simpler mechanics than those of limb movement, accurate control of NM position and velocity is still not a straightforward problem. The rabbit NM response is produced by retraction of the globe into the orbit, which displaces Harder's gland to force the NM across the cornea (Eglitis 1964). Globe retraction is produced primarily by the retractor bulbi muscle, which is innervated by motoneurons in the accessory abducens nucleus (Disterhoft et al. 1985). Restoration of globe and NM position depends on relaxation of the retractor bulbi muscle and the elastic properties of the globe, the NM, and surrounding tissues, rather than the action of an antagonist muscle. This quite complex chain of events linking motoneuron firing to NM movement suggests on a priori grounds that the relation between the two will not be simple (e.g., Delgado-García and Gruart 2005).
To examine these relationships explicitly Bartha and Thompson (1992a,b) constructed a detailed model of the rabbit NM response, which incorporated the mechanics of the globe, Harder's gland, the NM, and the motor units of the retractor bulbi muscle. Subsequent reimplementation of this model to examine the relation between motoneuron firing and NM response parameters (Mavritsaki et al. 2007) found that model performance depended critically on the nature of motoneuron recruitment (i.e., on the relationship between the number of motoneurons firing and the muscle force). If there were no recruitment and all motoneurons fired with the same rate (rate coding of muscle force), as investigated by Bartha and Thompson (1992a,b), then maximum NM displacement had a complex, sigmoidal relation to firing rate, reflecting the properties of retractor bulbi motor units (Lennerstrand 1974). If rate coding were combined with appropriate recruitment of motor units, however, peak NM amplitude and total motoneuron pool firing had a simple linear relation over the normal range of NM response amplitudes (Mavritsaki et al. 2007).
The question of how far this linearity is found in practice was addressed by analyzing electromyographic (EMG) records from the retractor bulbi muscle during the production of conditioned NM responses (Lepora et al. 2007). Spikes were extracted from the EMG record and the relationship of their firing rate to the NM response parameters was examined. Two relations were evident from this analysis: 1) the temporal profile of the EMG spike rate was well approximated by a Gaussian function, whose peak amplitude was linearly related to the peak amplitude of the conditioned response; and 2) the NM response profiles could be generated from the corresponding total EMG spike-rate profiles by passing the latter through a first-order linear differential equation equivalent to a simple first-order filter with time constant of the order of 100–200 ms. Over acquisition of the conditioned response, the peak EMG spike rate increased appropriately to drive a larger peak NM response, whereas the time-to-peak of the EMG spike rate and NM response amplitude did not change appreciably.
Therefore the overall effective dynamics that link motoneuron firing to NM movement can be viewed as a “plant simplification” that make conditioned response parameters such as amplitude, duration, and timing much easier to control than the complexities of the plant (i.e., muscle plus globe) would at first suggest. Such a linearized “virtual plant” should result in a lower computational load in the neuronal systems driving the movement and give more accurate motor control. Similar plant simplifications occur frequently in physiological motor control (Angelaki and Hess 2004; Mussa-Ivaldi et al. 1994; Nichols and Houk 1976), so the mechanisms by which they are achieved are of general interest.
Our previous analysis, which treated all EMG spikes as equal, irrespective of their amplitude (Lepora et al. 2007), cast no direct light on the underlying mechanism. Here we extend the analysis to take spike amplitude into account because inspection of retractor bulbi EMG records strongly suggested that spike amplitude and conditioned response parameters were related. Such a relation implies recruitment and, indeed, the analysis of EMG spike amplitudes indicated a simple procedure for recruiting motor units in the retractor bulbi muscle to control conditioned response (CR) amplitude. The same command is sent to all motoneurons in the accessory abducens nucleus pool. Its intensity varies over time with a Gaussian profile, whose peak is related to CR amplitude and width to CR duration. Motoneuron threshold varies with EMG spike amplitude, so that the higher the threshold, the larger the spike. Spike amplitude varies, in turn, with motor unit strength. In this manner, a simple high-level command can generate the forces necessary for the appropriately sized CR, using the common-drive mechanism for motor control previously proposed for a range of movements that use other muscles (De Luca et al. 1982a,b, 1996; Jabre and Spellman 1996; Masakado et al. 1995; Sauvage et al. 2006).
In these previous studies, EMG records allowed identification of the action potentials from individual motor units, so that unit firing rate could be related unequivocally to the amplitude of its action potential. Here, the technical difficulties of recording from the small and mobile retractor bulbi muscle prevent individual motor unit action potentials from being reliably identified. To determine whether the common-drive mechanism could in principle account for the EMG patterns observed in the previous study, we attempted to simulate the data patterns using a computational model. This model was able to characterize various observed features of the EMG records using simple assumptions about the properties of motoneurons receiving common-drive synaptic input. The success of the model is consistent with the view that, because of the technical difficulties of EMG recording and analysis, “computational models are fast becoming the dominant means to establish the mechanism relating EMG to muscle force and neural drive” (Keenan and Valero-Cuevas 2007). Some of the findings previously appeared in abstract form (Lepora et al. 2006).
METHODS
Data collection
Procedures for surgery, training, and data acquisition have been described previously (Lepora et al. 2007). Briefly, data were obtained from two different laboratories, one at University College London (Yeo) and the other at the State University of New York at Stony Brook (Evinger). All procedures carried out at University College London conformed to UK Home Office, Animal Procedures regulations. All procedures carried out at the State University of New York at Stony Brook adhered to federal, state, and university guidelines concerning the use of animals in research.
At University College London, EMG electrodes were implanted in the retractor bulbi muscle of three male Dutch belted rabbits (2.0–2.2 kg), subsequently referred to as RB1–RB3, under halothane general anesthesia. All general surgical procedures were as described in Yeo (1985) and those specific to the EMG electrode implantation are as previously described in detail (Lepora et al. 2007). EMG electrodes were of Teflon-coated stainless steel wire (50-μm diameter bare). In subject RB1, a pair of additional EMG electrodes, similar to those used in the retractor bulbi muscle, were implanted in the orbicularis oculi muscle immediately dorsal to the upper eyelid at a spacing of about 3 mm. Postoperatively, all subjects received buprenorphine twice per day for 3 days to provide postimplantation analgesia and minimize discomfort. Displacement of the NM was measured with an isotonic transducer based on a low-torque potentiometer (Gruart and Yeo 1995) whose output was fed into an A/D converter (CED 1401). For adaptation and training, each subject was placed in a restraining box in a sound-attenuating chamber facing a centrally mounted loudspeaker and the NM transducer was fitted. In the 50-min adaptation session, no stimuli were presented. For acquisition, the conditioned stimulus (CS) was a 1-kHz sinusoidal tone of intensity 85 dbA and duration 560 ms (RB1) or 410 ms (RB2, RB3), delivered through the loudspeaker. The unconditioned stimulus (US) was a 60-ms train of three biphasic current pulses of 2 mA delivered to the periorbital region via a pair of 3-mm Michel clips, attached 1–2 mm below the lower lid and 1–2 mm posterior to the lateral canthus, as previously described (e.g., Gruart and Yeo 1995). The intertrial interval between CS presentations varied randomly between 25 and 35 s and 100 trials were presented each session. On every 10th trial the CS was presented alone. On each trial EMG and NM response data were recorded for a 2-s period starting 1 s before CS onset. EMG data were sampled at 5–40 kHz and NM response data at 1 kHz.
After acquisition, RB1 and RB2 were tested with triplets, which were repeating sets of one trial CS alone, one trial CS + US, and one trial US alone. In addition, RB1 was tested with varying CS intensities (65, 75, 85, and 95 dBA). In the case of RB3, triplet testing was carried out with an airpuff US, which reduced CR amplitude in comparison with the earlier trials using a shock US that had been analyzed previously (Lepora et al. 2007). However, the triplet CRs were more evenly distributed across the amplitude range and were thus more suitable for analysis here.
General surgical procedures at the Stony Brook University laboratory were identical. To insert EMG electrodes, a pair of Teflon-coated stainless steel wires (76-μm diameter bare, 140 μm coated; A-M Systems #791000) were inserted into one slip of the retractor bulbi with a 30-gauge hypodermic needle. Five hundred microns of insulation was removed from the tip of each wire and tip separation between the wires was between 1 and 2 mm. To eliminate any postoperative discomfort, rabbits received buprenorphine twice daily for 2 days after the surgery. Data were analyzed from rabbit RB4, a New Zealand white. After initial acquisition of the conditioned response, RB4 was tested under a variety of conditions, including variations in CS and US intensities and extinction.
Once CR amplitude increased beyond the chosen threshold amplitude of 0.5 mm, it reached asymptotic values in a small number of trials (median 45). Since CS alone trials constituted only 10% of all trials during acquisition, few CRs from acquisition were available for analysis. Thus most of the CRs analyzed here (median 93%) came from the postacquisition phase.
Data selection
The data used for analysis came from trials in which the CS was presented on its own and CR amplitude was >2.5% of the maximum CR amplitude for that subject. About 15% of these trials were subsequently discarded either 1) because there was spontaneous movement of the NM just before CS presentation or 2) there was interference in the EMG or NM response signals or there was evidence from the EMG recording suggesting electrode movement during the CR.
Data analysis
The occurrence of motor-unit action potentials was estimated from spikes in the EMG records (Sanders et al. 1996) using a threshold crossing procedure so that a spike was registered whenever EMG amplitude increased past a minimum spike amplitude threshold θ (Fig. 1A). Given an EMG signal x = (x0, x1, x2, …) sampled at frequency fEMG Hz, [i.e., at discrete time tEMG = (0, 1/fEMG, 2/fEMG …)], this procedure gives a spike sequence s that is a binary time series of zeros and ones
Example of spike extraction and classification for electromyogram (EMG) of retractor bulbi muscle (subject RB4). A: a portion of EMG record shown with expanded time base to illustrate threshold crossing criterion for spike extraction and classification for 5 spike amplitude classes. B: entire EMG record for a conditioned-stimulus alone trial (middle trace), showing extracted and classified spikes (bottom trace) and conditioned nictitating membrane response (NMR; top trace). The record is the same as that used in Fig. 1 of Lepora et al. (2007), with A shifted by −50 ms to show the tallest spikes. (The threshold was erroneously shown at too high a value in that figure, although the correct value of 0.0375 mV was used correctly elsewhere in the study.)
Possible problems with this method were discussed in a previous analysis of the same EMG data set (Lepora et al. 2007). First, high-frequency noise can cause multiple-spike counting artifacts, so the EMG records were downsampled to a frequency of 5 kHz. Second, at high spiking rates there can be significant spike overlap, producing a wavelike interference pattern rather than discrete action potentials (Sanders et al. 1996). The effects of this interference were demonstrated to be minor for the EMG data sets under consideration (Lepora et al. 2007). Third, a suitable threshold θ must be chosen for each subject. A signal-to-noise separation procedure was used to estimate the values for θ, which are reproduced here in Table 1 for each of the four subjects RB1–RB4.
Table 1.
Parameters describing the observed and fitted EMG spike-rate profile and CR profile for subjects RB1–RB4
| Subject | Spike Amplitude | Total Spike Rate Gaussian Fits | NM Response Amplitude Range y, mm | ||||
|---|---|---|---|---|---|---|---|
| Threshold θ0, mV | Maximum εmax, mV | Median Center μ, s | Median Width σ, s | Peak Range h, spikes/s | Median Fit R2 | ||
| RB1 | 0.17 | 2.53 | 0.45 | 0.13 | 0–440 | 0.92 | 0–4.4 |
| RB2 | 0.15 | 0.75 | 0.29 | 0.06 | 0–220 | 0.94 | 0–2.4 |
| RB3 | 0.0275 | 0.67 | 0.36 | 0.10 | 0–450 | 0.83 | 0–3.3 |
| RB4 | 0.0375 | 0.36 | 0.39 | 0.07 | 0–370 | 0.82 | 0–6.1 |
As explained in the introduction, the focus of the present study is EMG spike amplitude. The distribution of spike amplitudes was quantified by placing the peak amplitude ε of each of the spikes found by Eq. 1 in one of Nclass distinct spike amplitude classes. These classes were defined from a set of equally spaced thresholds θ1 < θ2 < … < θNclass, such that a spike of amplitude ε is in class k (where 1 ≤ k ≤ Nclass) if
where θ0 is the spike-detecting threshold used in Eq. 1 and εmax is the maximal spike amplitude across all EMG records under consideration (values for εmax are given for RB1–RB4 in Table 1). A characteristic amplitude ![[epsilon with macron]](/corehtml/pmc/pmcents/x03B5x0304.gif)
k for the spikes in each amplitude class is the midpoint of their range
The above procedure results in Nclass spike sequences (s0(k), s1(k), s2(k) …) with each a binary time series of whether spikes occur in a particular amplitude range—for example, with Nclass = 5 classes, the spike classes k = 1, 3, 5 can be thought of as signaling when the “small,” “medium,” and “large” spikes occur, as illustrated in Fig. 1. It is important to note that these classes are a convenient measure for dividing spikes by amplitude along a continuum of values, rather than representing true clusters of spikes into distinct classes.
Once spikes had been extracted, their firing rates were calculated from the number of spikes in successive 50-ms time intervals. The noise levels of individual records meant that data had to be pooled across trials, based on similarity of CR amplitude. Three-trial batches were typically used (Lepora et al. 2007) to give a suitable number of spikes for analysis, with the trials ordered by peak NM amplitude to ensure that similarly sized responses contributed to the average. All results in this study were checked for robustness with respect to 5, 6, 7, 8, 9, and 10 classes; results were not found to change appreciably apart from a general deterioration of the results as the number of classes was increased because of too few spikes in each class. Thus for simplicity the number of spike classes was always set to Nclass = 5, a value small enough for the results to have enough spikes in each class and large enough for statistical analysis of the EMG records with respect to spike height.
Previous analysis of the EMG data indicated that for CRs a portion of the total spike rate profiles could be well fitted by Gaussian curves (Lepora et al. 2007) according to the equation
where f(t) is instantaneous firing rate calculated over successive 50-ms intervals, μ is the time of the distribution mean, h its maximum amplitude, and σ is a measure of its width. Fitting was restricted to the interval from CS onset to one interstimulus interval (ISI) after US onset (ISI is the time between CS and US onsets), to exclude a tail in the spike rate profile of some records. This range actually covers most of the spike rate variation of each record because the peak frequency typically occurs significantly before US onset (see Fig. 2, for example). Effectively, this procedure characterized each spike distribution by three parameters estimated from the firing rate data set as follows. The mean μ is the time-weighted mean of the spike rate
where fj represents the firing rate of the jth bin from a total of n bins, tj is the time at which that firing rate occurs, and N is the total number of spikes divided by the time width δ of the bin. The maximum amplitude h is the peak value across the spike rate record. The width σ is found by equating the spike total of the record to that for the Gaussian curve (i.e., matching areas under the curves) through the relation
The goodness of fit for the function in Eq. 4 is then the R2 value
which describes the proportion of variance accounted for by the fit of f(t) to the spike rate data. An associated significance level can be calculated using the variable R
Common drive model of motoneuron firing. Input to model is the common-drive synaptic current, which is distributed equally across 100 simplified model neurons. These model motoneurons have firing rates that are linearly proportional to the input synaptic currents above a threshold R (the rheobase current) with gain G. Each model motoneuron is assumed to control a specified population of muscle fibers—the motor unit. Assuming each motor unit generates an EMG spike with characteristic amplitude (here assumed proportional to the rheobase current), these motoneuron firing rates can be used to simulate an EMG record for the muscle in response to the common-drive current.
Here we also applied this Gaussian-fitting analysis to the firing rate profiles of the different EMG spike amplitude classes according to
where 
k(t) is the firing rate of class k. This results in a set of three parameters (μk, σk, hk) describing the spike distribution of each of the Nclass spike amplitude classes, with a fitting parameter rk2 describing the goodness of fit.
Finally, the relationships between the values for instantaneous total firing rate f(t) and the class firing rates k(t) were analyzed for each data set. It proved possible to fit these relations with straight lines of the form
where gk is the gradient of the best-fit line for spike amplitude class k, 
k is the base total-spike rate representing the intercept above which k(t) > 0, and the operation [x]+ truncates its argument to positive values.
Simulation of EMG data by common-drive mechanism
The basic structure of the motoneuron pool model is shown in Fig. 2. Motoneurons in the accessory abducens nucleus were represented individually. Each motoneuron had its own threshold and gain, but was driven by an input common to all the motoneurons. EMG spike amplitude was related to motoneuron threshold, so that EMG records could be simulated and then analyzed in the same way as the actual EMG records. The model is a simplified version of EMG models described previously (e.g., Fuglevand et al. 1993; Keenan and Valero-Cuevas 2007; Keenan et al. 2006; Zhou and Rymer 2004), with the additional features of 1) a Gaussian-profile common drive and 2) a grouping of motoneuron outputs on the basis of EMG spike amplitude.
The model is first described with parameters appropriate for data set RB1. The manner in which parameters were varied to fit the other data sets (RB2–RB4), and then to test model robustness, are subsequently described.
GAUSSIAN-PROFILE COMMON DRIVE.
Analysis of EMG spikes recorded from the retractor bulbi muscle during the production of a set of conditioned responses at various stages through conditioning indicates that their overall firing rate has a temporal profile that is approximately Gaussian in shape (Lepora et al. 2007). Accordingly, the common-drive input current I(t) to each motoneuron was given a Gaussian profile
where μ is the time at which the current reaches its peak amplitude I0 and σ is the width of the Gaussian function. The parameters μ and σ are independent of the amplitude of the conditioned response, whereas I0 is assumed to be directly proportional to the peak amplitude of the conditioned response (Lepora et al. 2007).
Parameter values for μ and σ were taken directly from the EMG spike firing rates for data set RB1, giving μ = 0.43 s and σ = 0.14 s (values given in Fig. 6, A and B). However, the third parameter I0 is not specified by available measurements, and so has to be initially assigned a plausible value. Given the maximum possible NM response amplitude is 15 mm (references in Mavritsaki et al. 2007), a numerically convenient value for the I0 corresponding to that response is 15 nA. The maximum value of I0 (=Imax) for a given data set is then simply determined by the amplitude of the maximum conditioned response in that data set. For RB1 this was 4 mm, giving Imax = 4 nA.
Recorded and simulated EMG spike-rate profiles for different spike-amplitude classes, corresponding to conditioned NM responses of different amplitude for subjects RB1 and RB4. A–F: EMG data. Each trace is a mean of 3 EMG firing rate profiles for spike amplitudes partitioned into 5 classes. An individual profile was calculated every 50 ms as the firing rate of the spikes in a given amplitude class in a surrounding 50-ms bin. In each panel, CS onset is at time 0 and US onset is denoted by a dotted vertical line. Gaussian curves were fitted to traces (dashed lines). G–L: model EMG results. Model results are generated for 10 model trials with peak input currents I0 ranging from 0 to 4 nA. The small, medium, and large responses are the 2nd, 4th, and 9th model trials, respectively. All other simulated conditions are as for the data traces. Note that because the Gaussian fits (dashed lines) well approximate the model results (solid lines) the 2 curves almost overlie.
PROPERTIES OF MOTONEURON POOL.
The method for calculating the firing rate of a motoneuron given its synaptic input current and intrinsic properties was taken from previous simplified motoneuron models (Binder et al. 1993; Dean 1997; Heckman and Binder 1990). These models assume that an individual motoneuron starts to fire when its intrinsic threshold current, termed the rheobase R, is exceeded by the effective synaptic current I(t). Instantaneous firing rate f(t) for a motoneuron is then given by
where Gi is a constant representing the intrinsic gain of the motoneuron, as indicated by the slope of the relation between firing rate and injected current (Binder et al. 1993), and the notation [x]+ corresponds to x for x > 0 and 0 for x ≤ 0; the index i labels the motoneuron and ranges between one and the total number of motoneurons nMN (here set to 100, from observations in cat; Lennerstrand 1974).
In accordance with previous studies (Fuglevand et al. 1993; Keenan and Valero-Cuevas 2007; Zhou and Rymer 2004) the gains Gi were assumed to be the same (=G) for all motoneurons (cf. properties of abducens motoneurons as discussed in Dean 1997) and the thresholds Ri varied, to give many low-threshold motoneurons and progressively fewer high-threshold motoneurons. An exponential distribution is convenient for this purpose
where Rmax is the maximum threshold and the exponent q defines the proportion of small to large rheobase currents.
The combined values for G and Rmax are constrained by the previous assignment of 15 nA for the value of the input current that produces the maximum possible NM extension of 15 mm. For the initial model Rmax was set to 7.5 nA. The last recruited neuron thus starts to fire at this value of input current, and must reach its maximum effective firing rate when the input current reaches 15 nA (i.e., over a range of 7.5 nA). Since the isometric forces produced by motor units in the retractor bulbi muscle saturate at stimulation frequencies of about 150 spikes/s (Lennerstrand 1974), motoneuron gain G must be 20 spikes·s−1·nA−1. This can be compared with measured values of 20–40 Hz/nA for cat abducens motoneurons (Grantyn and Grantyn 1978). Finally, q in Eq. 12 was assigned the value of 3.45 by trial and error. The relation of this parameter to EMG properties was subsequently explored (details in the following text).
The properties of the model motoneuron pool produced by Eqs. 11 and 12 with these parameters are shown in Fig. 3. Because Eq. 12 describes an exponential relation between motoneuron number and threshold (Fig. 3A), the actual distribution of motoneuron thresholds is hyperbolic (Fig. 3B). This distribution produces a relationship between input current and total pool firing rate that is close to linear over the range of input currents for the RB1 data set (Fig. 3D).
Properties of the model of the motoneuron pool with parameters derived from the RB1 data set (details in methods). A: exponential relation (Eq. 12) between motoneuron threshold (rheobase current) and the index i of the motor unit, numbered from 1 (smallest threshold) to 100 (largest threshold). For the RB1 data set the largest input current was 4 nA, so that neurons with a rheobase higher than this value were not recruited. B: distribution of motor units with respect to rheobase current, derived from the relation shown in A. The histogram indicates the number of units in each of 5 equal bins over the range of rheobase currents 0–4 nA. The exponential relation in Eq. 12 gives a hyperbolic distribution of motor unit numbers (best-fit hyperbola for a histogram with 5 intervals y = 41/[1 + 1.86(R − 0.4)] shown as black curve on plot), as expected theoretically by integrating Eq. 12 over motor unit number. C: linear relation between rheobase current and EMG spike amplitude. The x-intercept of this linear relation is the minimum threshold θ0 = 0.17 mV for EMG spike detection in data set RB1. D: relation between common drive input current and total motoneuron firing rate from model Eqs. 11 and 12. Note the weak deviation from linearity for small input currents.
EMG GENERATION.
Simulated EMG records were generated from model firing-rate profiles as follows. The shapes of the individual EMG spikes, and possible interference between spikes, were ignored because analysis of the actual EMG records for the retractor bulbi muscle suggested that for the data under consideration interference effects were not important (Lepora et al. 2007). The relationship between spike amplitude and motoneuron threshold was assumed to be linear
where εi is the height of spike produced by the ith motor unit, θ0 is the threshold height for spikes to be counted, εmax the maximum possible spike height, and Imax is the maximum input current for the RB1 data set. Finally, a sampling parameter 0 < P < 1 was introduced to convert summed motoneuronal firing to summed EMG spike firing
In the absence of any prior evidence on the motor units sampled, it was assumed that the subsample of measured units was unbiased across the motoneuron pool, given the small size of the retractor bulbi muscle, so that p can be regarded as a constant for any given data set of EMG records.
For the RB1 data set θ0 = 0.17 mV, εmax is the maximum observed spike height of 2.53 mV and Imax = 4 nA. These parameters give the linear relation between spike amplitude and motoneuron threshold depicted in Fig. 3C. The sampling fraction p was set to 0.09, determined by the ratio of the observed peak EMG firing rate (440 spikes/s) to the summed model firing rate for RB1's largest conditioned response (peak amplitude 4 mm).
FITTING DATA SETS RB2–RB4.
A subset of model parameters was varied to fit the remaining three data sets RB2–RB4. These were the parameters I0, μ, and σ for the Gaussian-profile common drive (Eq. 10), the parameter q for the distribution of motoneuron thresholds (Eq. 12), and the EMG parameters θ0, εmax, and p (Eqs. 13 and 14). Appropriate values of these parameters gave fits to the data sets RB2–RB4 similar to those shown for data set RB1. Fits for subjects RB1 and RB4 are shown in results (Figs. 66–8) and fits for subjects RB2 and RB3 in Supplemental Figs. S1–S3.1
) Experimentally determined EMG threshold θ0 = 0.15, 0.0275, and 0.0375 mV and maximum EMG spike height εmax = 0.75, 0.67, and 0.36 mV for data sets RB2–RB4 (Table 1). Overall model behavior is unaffected by changing these parameters because the class- and total-firing rates do not depend on these EMG parameters (both the spike amplitudes and class boundaries are scaled identically, so the motor units remain in the same classes).
) The Gaussian mean μ and width σ of the common drive current were set equal to the corresponding values for the mean and intercept width of the data for RB2–RB4, which are μ = 0.28, 0.35, and 0.31 s and σ = 0.10, 0.11, and 0.10 s, respectively (Figs. 7 and and8).8). Because the model Eqs. 10 to 14 depend implicitly on time, the model output remains of the same functional form in Eqs. 16 and 17 but with the precise temporal profiles rescaled by the changes to the model input or common drive current.
) The motoneuron sampling proportion p is fixed to the peak model firing rate; peak firing rates for data sets RB2–RB4 are 220, 450, and 370 spikes/s (see Table 1), giving p values of 0.14, 0.10, and 0.053, respectively. Overall model behavior is again robust to changing the sampling proportion because only the absolute firing rates are affected, not the functional form of Eqs. 16 and 17.
) The parameter q specifies the precise exponential relation between motoneuron number and threshold (Eq. 12) and its value determines the parameters u, v in Eqs. 16 and 17 (see also Figs. 7 and and8)8) in a complicated manner that depends on the EMG threshold θ0 and maximum EMG spike height εmax. Appropriate values of q were determined by trial and error for the data sets RB2–RB4, giving 1.8, 4.4, and 1.9 (see also robustness in the following text).
Parameters of Gaussian fits to recorded and simulated EMG spike-rate profiles as functions of spike amplitude for subjects RB1 and RB4. A–F: EMG data analysis. Data points correspond to the mean (A and D), width (B and E), and peak firing rate (C and F) of the Gaussian fits in A–F of Fig. 6 for each spike-amplitude class, plotted against the central value of spike amplitude for that class. Dashed lines join data points from trials with similar values of CR amplitude (results averaged over 3 trial batches). US onset is denoted by the dotted horizontal line on A and D. The best-fit line and r2 value (where r is the correlation coefficient) are shown on B and E. The best-fit exponential and its r2 value are shown on C and F. The ratio of peak firing rates hk/h1 is used in C to show the dependence of peak class-firing rate on spike amplitude independently of the peak total firing rate. G–L: model EMG results. Results generated similarly to the EMG data, for the simulated EMG results in G–L of Fig. 6. All fits shown here were significant (P < 0.01, calculated as described in methods).
Recorded and simulated change in class spike rate with total spike rate for subjects RB1 and RB4. A and D: class spike rate plotted against total spike rate sampled at 50-ms intervals. For each spike amplitude class, this relation is approximately linear (linear fits denoted by the solid lines; data by the points). B and E: base firing rates (total-spike rate where a class spike rate becomes nonzero) plotted against the central spike amplitude for each of the spike amplitude classes in A and D. The solid line is the best linear fit to these data, passing through zero base firing rate when the spike amplitude equals the EMG threshold θ0. C and F: the gradients of the best fits in A and D plotted against the central spike amplitude for each spike amplitude class. Hyperbolas (see Eq. 17) were found to give good fits to these data. G–L: model EMG results. Analysis of the simulated EMG records is the same as that for the EMG data in A–F. Data fits shown in central and right columns were significant (P < 0.05, calculated as described in methods).
PARAMETER VARIATION: MODEL ROBUSTNESS.
Parameters were subsequently varied systematically to investigate the robustness of the model. For this purpose the appropriate parameters were those describing model structure in Eqs. 11 and 12. Nonlinear relations between EMG spike amplitude and rheobase current were also considered, such as exponential and power laws, but no significant improvement to model performance was achieved for the corresponding increase in model complexity.
RESULTS
Analysis of EMG data
For clarity, the data analysis and modeling results (Figs. 66–8) are presented here for only two of the four subjects: RB1 (a Dutch belted rabbit) and RB4 (a New Zealand white). The corresponding results for the remaining two subjects RB2 and RB3 (both Dutch belted rabbits) are shown in Supplemental Figs. S1–S3. The two subjects shown here are sufficient to illustrate the similarities and differences between the data sets. Results across all four subjects were consistent with the modeling study.
FIRING-RATE PROFILES FOR TOTAL EMG SPIKES AND RELATION TO CONDITIONED NM RESPONSES.
Previous analysis of retractor bulbi EMG (Lepora et al. 2007) showed that during the production of conditioned responses the firing-rate profiles for total EMG spikes have a particular shape that gradually rises from its onset, reaches a peak, then gradually falls to its offset (representative record shown in Fig. 4 A, solid line). The peak times of the firing-rate profiles for individual trials are similar, roughly 50–100 ms before US onset. The shapes of the firing rate profiles across all trials can be well fitted by Gaussian functions (Fig. 4A, dashed line), as explained in methods.
Representative EMG spike-rate profile and conditioned response (CR) profile for subject RB1. A: temporal profile of EMG total-spike rate and its best-fit Gaussian profile. Conditioned stimulus (CS) onset is at time 0 and unconditioned stimulus (US) onset is denoted by the dotted vertical line. The trace is the mean of 3 EMG spike-rate profiles. B: temporal profile of conditioned nictitating membrane (NM) response and the best-fit derived from a linear filter model (Eq. 1). The 2 traces are means over 3 responses corresponding to the records averaged over in A. C: predictions from linear-filter model plotted against data points for all conditioned NM responses for RB1 (sampled every 1 ms for each response). The close match shows the linear model is a good fit to the data. All 3 fits shown here were significant (P < 0.01, calculated as described in methods).
The conditioned NM response also has a stereotyped shape, with a steep rise from onset to its peak and then a gradual fall back toward zero displacement (Fig. 4B, black line). Separate trials have a similar latency-to-peak time that coincides with US onset. The NM response profiles can be generated from the corresponding total EMG spike-rate profiles by passing the latter through a first-order linear differential equation
where FEMG is EMG spike firing rate, y is NM position, dy/dt is NM velocity, and K and R are constants related to the elasticity and viscosity respectively of the overall system. This system is equivalent to a first-order filter with time constant R/K of about 200 ms and with gain fixed across all trials of each data set (Fig. 4, B and C). Thus analysis of total EMG-spike firing rates in relation to conditioned NM response profiles suggests that the nictitating-membrane plant (retractor bulbi muscle, globe, Harder's gland, and membrane) behaves effectively as a simple linear system, even though its complicated anatomical structure would incorrectly suggest more complicated dynamics.
PEAK EMG AMPLITUDE IS PROPORTIONAL TO PEAK CONDITIONED NM RESPONSE.
The analysis referred to earlier does not address the issue of motor-unit recruitment, a process that modeling studies (Mavritsaki et al. 2007) suggest could make an important contribution to plant linearity. To address this issue, we analyzed the properties of EMG spike amplitude. As can be seen from Fig. 5, peak EMG spike amplitude is related to peak conditioned response amplitude in a roughly linear manner over the range of response amplitudes studied (up to ~6 mm). Insofar as EMG spike amplitude is related to motor-unit strength (discussion), these results are consistent with the view that conditioned NM response amplitude is controlled by recruiting motor units of increasing strength, at least for the range of response amplitudes investigated.
Peak EMG amplitudes plotted against peak conditioned NM responses for subjects RB1–RB4. A–D: each data point is a mean of 3 EMG records of similar CR peak amplitude (thus the peak values appear less than those in Table 1). The relationship between peak EMG amplitude and peak NM amplitude is approximately linear for all subjects (best linear fit constrained to pass through EMG-spike threshold θ0 at zero response amplitude). There is no evidence for saturation of EMG spike amplitude over the range of response amplitudes available. All 4 fits shown here were significant (P < 0.02, calculated as described in methods).
A previous study found that peak conditioned response amplitude linearly relates to peak EMG spike rate over the same range of response amplitudes (Fig. 5; Lepora et al. 2007). Therefore peak EMG spike-rate and peak EMG spike amplitude are also related in a roughly linear manner. The precise relation between EMG spike amplitude and EMG spike rate will be qualified over the next few sections, to suggest the pattern of motor unit recruitment involved in controlling the conditioned NM response. Because the relation between temporal profiles of instantaneous EMG spike rate and NM conditioned response amplitude is described by the linear differential Eq. 15, this motor unit recruitment strategy will be more complex than a simple dependence on NM position during the CRs.
FIRING-RATE PROFILES FOR EMG SPIKES OF DIFFERENT AMPLITUDES.
Subsequently more detailed analysis of EMG spike amplitudes showed that spikes of different sizes fired in fixed relation to one another. As described in methods, EMG spikes were partitioned into five amplitude classes and representative firing-rate profiles for those classes are shown in Fig. 6.
Four regularities of the retractor bulbi EMG are illustrated in these records. First, larger spikes only appear in conjunction with larger conditioned responses, consistent with the results shown in Fig. 5. Second, when larger spikes do appear, they are always preceded and succeeded by smaller spikes (see also Fig. 1). Third, when spike firing rates are compared across trials, the most prominent feature is that the firing rates of the spike amplitude classes change in step with each other. Finally, just as for overall spike frequency (Fig. 4A), the frequency profiles for individual spike classes can also be fitted by Gaussian functions. The firing rates for each class increase gradually in between the trials shown (see also relation to total emg).
The Gaussian fits to the spike rates are convenient because they allow the firing-rate profile for each spike-amplitude class k to be described by just three parameters: a mean time of spike occurrence μk (which equals the time of peak firing for a Gaussian profile), a measure of profile width σk, and the peak firing rate hk (methods, Eq. 8). This allows the regularities outlined qualitatively above to be described quantitatively in terms of those parameters (Fig. 7). Thus 1) the mean spike times μk for all the classes are timed about 50 ms before US onset (Fig. 7, A and D), 2) the firing widths σk decrease linearly with spike amplitude to zero value at maximum spike amplitude εmax (Fig. 7, B and E), and 3) the peak firing rates hk decrease exponentially with spike amplitude (Fig. 7, C and F). These relations can be summarized by a three-equation system
where l and u are an empirical gradient and exponent, respectively (values shown in Fig. 7) and k is the midpoint amplitude for the spikes in amplitude class k.
Although the relations in Eq. 16 give good fits to the data shown in Fig. 7, there is still an appreciable amount of scatter around the lines of best fit. This arises partly from random variations in the timing of the conditioned responses themselves (references in Kehoe et al. 2008; Lepora et al. 2007), which give rise directly to variations in the timing of the inferred motor commands (Lepora et al. 2007; Figs. 11 and 12). For example, the time-to-peak conditioned NM response amplitude is known to vary randomly from trial to trial, as is evident in the variation of the Gaussian means in Fig. 7, A and D. However, a second source of scatter comes from pooling data across different CR amplitudes. For example, it can be seen from Fig. 6 that the width of the best-fit Gaussian increases with response amplitude for a given amplitude class. This source of scatter can be reproduced by the deterministic common-drive model described in the following text.
RELATION TO TOTAL EMG SPIKE RATE.
Because of these regularities in the structure of the retractor bulbi EMG, the firing rates of each individual amplitude class are related to total spike rate in a systematic manner. This is shown in Fig. 8, with firing rates calculated over 50-ms intervals, as described in methods.
The firing rate for each individual amplitude class is approximately linearly related to total spike firing rate (Fig. 8, A and D). The lines of best fit vary with amplitude in two ways. First, the intercept with the x-axis, which corresponds to the total spike rate at which spikes of the amplitude class appear, increases approximately linearly with spike amplitude (Fig. 8, B and E). This corresponds to the qualitative observation that spikes appear and disappear in strict order of amplitude (Fig. 6). Second, the gradients of the lines decline with amplitude (Fig. 8, C and F), reflecting that high-amplitude spikes always fire more slowly than smaller ones (Fig. 6). This decline is fitted well by a hyperbolic relation, which is consistent with an exponential distribution of motor unit thresholds (see Modeling in the following text).
These regularities can be summarized in a system of two equations
where the gradient m and v are empirical parameters (values given on appropriate panels of Fig. 8), k is the midpoint spike amplitude for the spike amplitude class k defined in Eq. 3, and gk is the gradient of the class versus total spike rate line.
Simulated EMG records
The simulated EMG records produced by the common-drive model (Fig. 1) were analyzed in the same way as the actual EMG records described earlier. The model was used to generate EMG records for 10 different peak input current amplitudes equally distributed over the range 0–4 nA (methods). To aid comparison with results from the EMG data analysis, the model results are shown for subjects RB1 and RB4 in Figs. 66–8 and for subjects RB2 and RB3 in Supplemental Figs. S1–S3.
PROPERTIES OF SIMULATED AMPLITUDE CLASSES.
Model EMG spikes were partitioned into five amplitude classes in the same way as the recorded spikes. This procedure resulted in a firing rate k(t) for each of the five amplitude classes.
Results of the simulation (Fig. 6, G–L) are in good agreement with the experimental data, both qualitatively, as comparison with the experimental data indicates (Fig. 6, A–F), and quantitatively (see following text). The four regularities observed in the experimental data are all present in the simulated records: 1) larger spikes appear only in conjunction with larger input currents (corresponding to larger conditioned response); 2) when they do appear, larger spikes are always preceded and succeeded by smaller spikes; 3) when firing rates are compared across trials, the most prominent feature is that the firing rates of the spike amplitude classes change in step with each other; 4) the spike-rate profiles for individual amplitude classes can be fitted by Gaussian functions.
Describing these Gaussian fits in terms of the mean time of spike occurrence, profile width, and the peak firing rate leads to the three quantitative relations (Eq. 16) for the simulated data sets that were obtained previously for the experimental data sets. In particular: 1) the time of peak firing is the same for all amplitude classes (Fig. 7, G and J; first relation in Eq. 16), 2) the widths of the Gaussian fits to firing rate decrease linearly with spike amplitude to zero value at maximum spike amplitude εmax (Fig. 7, H and K; second relation in Eq. 16), and 3) the peak firing rates of the Gaussian fits decrease exponentially with spike amplitude (Fig. 7, I and L; third relation in Eq. 16).
Both simulated and experimental data show scatter around the lines of best fit. In the case of the experimental data, some of this is caused by random variability that is not present in the modeled data. However, there is also variability that is common to the two sets, which appears to arise from lumping together EMG records corresponding to different amplitudes of conditioned response. The similarity of the scatter patterns in experimental (Fig. 7, B and E) and simulated data (Fig. 7, H and K) provides further support for the view that the patterns observed in the experimental data are consistent with a common-drive mechanism for generating CRs.
RELATION TO SIMULATED TOTAL EMG SPIKE RATE.
The experimental results also showed that the firing rates of spikes of different amplitude were related to the total-spike firing rate in a systematic manner (Fig. 8).
It can be seen that the general pattern of simulated results (Fig. 8, G–L) are similar to the experimental results (Fig. 8, A–F). The firing rates for each amplitude class are approximately linearly related to total firing rate (Fig. 8, G and J). The deviation from linearity seen for the lowest-amplitude spike class arises because the firing-rate for each class in the model is linearly related to input current, whereas the relationship between input current and total spike firing-rate is not precisely linear (Fig. 3D) because of recruitment. There are suggestions of a similar curvilinear relation between the firing rate of the lowest amplitude spike class and total firing rate in the experimental data (Fig. 8, A and D), although these data show scatter not seen in the deterministic model arising from random variability in the actual EMG recordings. The relationship between spike amplitude and threshold for firing, expressed in terms of total spike firing rate, is also approximately linear (Fig. 8, H and K). Again slight deviation from linearity arises from the model neuron firing-rate being linearly related to input current, not total spike firing-rate. The corresponding plots for the experimental data also show approximately linear relations (Fig. 8, B and E), although the scatter is too great to observe any small systematic deviations from linearity. Finally, the relation between the firing-rate gradient and spike amplitude from each class is fitted well by a hyperbolic function (Fig. 8, I and L), as expected from an exponential relation between motoneuron number and threshold (Fig. 3), and as observed in the experimental data (Fig. 8, C and F).
The relations apparent in Fig. 8, G–L can be summarized by the quantitative relations of Eq. 17 found to describe the experimental data.
MODEL ROBUSTNESS.
Most of the parameters that were varied to fit data sets RB1–RB4 related to particular features of the individual data, rather than central aspects of the model itself. To assess the robustness of the model, the effects of varying the main model parameters (Eqs. 11 and 12) were therefore investigated.
1)Number of motor units n. Changing n from its baseline value of 100 altered the total- and class-firing rates equally, which rescaled hk and bk in Eqs. 16 and 17 while keeping the form of the overall relations shown in Figs. 7 and and8.8. If the number of units was reduced
![[less, similar]](/corehtml/pmc/pmcents/lsim.gif)
20 the fits to relations in Eqs. 16 and 17 were degraded by discretization errors in separating the units into spike amplitude classes. Model behavior is therefore robust to changing motor unit number provided a reasonable motor unit number is considered.2)As explained in methods, the magnitude of the input current I(t) is determined by the assignment of a maximum value to the input current (Imax) that produces the maximum NM extension. Once Imax is fixed and a value is also given to Rmax, the highest rheobase threshold in the motoneuron pool, then the motoneuron gain G is also determined (methods). The value of Rmax in fact determines how much of the force produced by the muscle is produced by recruitment and how much by rate coding. The linear relation between EMG spike amplitude and peak CR amplitude observed for the experimental data (Fig. 3) indicates that for each of the modeled data sets recruitment continues throughout the amplitude range, that is Rmax > Imax. With this restriction the overall pattern of simulated EMG results is insensitive to the value of Rmax (and corresponding value of G), provided the sampling parameter p is tuned appropriately to match a given experimental data set. This was demonstrated in simulation (results not shown), but can also be seen directly from the definition of the model because the firing rates of all amplitude classes scale by the same amount on scaling I(t) and Ri in Eq. 11.
3)The final parameter to be varied systematically was q, which determines the precise exponential relation between motoneuron number and threshold (Eq. 12). The exponent q in the rheobase distribution determines the parameters u, v in Eqs. 16 and 17 in a complicted manner. We have verified that the model results fit the general form of Eqs. 16 and 17 for q values ranging from 1 to 10.
Additional observation
Subject RB1 had an additional EMG electrode implanted in the orbicularis oculi (OO) muscle, so analysis of OO EMG was also possible for this animal. The pattern of results for spikes of different amplitudes was very similar to that seen for the retractor bulbi EMG.
DISCUSSION
The purpose of the present study was to investigate the control of conditioned NM responses, by analyzing and modeling spike amplitude in EMGs recorded from the retractor bulbi muscle during conditioning. Four main features were revealed and quantified: 1) spike amplitude increased with response amplitude; 2) smaller spikes always appeared before larger spikes; 3) subsequent firing rates covaried for spikes of different amplitude, with smaller spikes always firing at higher rates than larger ones; and 4) subsequent firing-rate profiles were approximately Gaussian for all amplitudes, with parameters related to conditioned NM response parameters so that larger conditioned responses were accompanied by larger spikes. These apparently complex features of the retractor bulbi EMG could be reproduced by a simple model of common-drive–based recruitment in the retractor bulbi muscle, consistent with previous suggestions that the control signal for conditioned NM responses codes the topography of the response in a straightforward manner.
Data analysis
RECRUITMENT AND EMG AMPLITUDE.
To our knowledge, spike amplitude in the retractor bulbi EMG has not been analyzed systematically in previous studies. The analyses here indicate a number of respects in which the retractor bulbi EMG resembles that of other muscles. One is that small spikes appear before large ones as muscle activation increases. This observation has been made in many studies (see Goldberg and Derfler 1977 for introduction to earlier literature; Akaboshi et al. 2000; Durkovic 2006; Jabre and Spellman 1996) and is generally interpreted in the context of the size principle, according to which weak (small) motor units are recruited before strong (large) ones (Henneman and Mendell 1981). If the amplitude of an EMG spike reflects the amplitude of a motor unit action potential (MUAP), which in turn is related to the number of fibers in a motor unit and/or the dimensions of individual fibers (Goldberg and Derfler 1977; Henneman and Mendell 1981), then the appearance of small spikes first is consistent with the size principle.
One caveat for this interpretation is that absolute spike amplitude is also heavily dependent on an extrinsic factor, namely electrode type and position. It is relative amplitude that exhibits the robust recruitment effect, which suggests that it is primarily influenced by factors intrinsic to the motor unit, such as the number of fibers it contains (Henneman and Mendell 1981) or the cross-sectional area of individual fibers (Goldberg and Derfler 1977). The precise effects of intrinsic and extrinsic factors in determining EMG spike amplitude have been extensively discussed and modeled (e.g., Akaboshi et al. 2000; Clamann and Robinson 1985; Ertas et al. 1995; Fuglevand et al. 1993; Keenan and Valero-Cuevas 2007; Keenan et al. 2006; Krutki et al. 2008; Nandedkar et al. 1988a,b; Stålberg and Karlsson 2001; Yao et al. 2000). In general terms, it appears that intrinsic factors are likely to be particularly dominant in small muscles, especially if they contain only one type of muscle fiber, as appears to be the case for the retractor bulbi muscle (Alvarado et al. 1967; Crandall et al. 1981; Gurahian and Goldberg 1987; Lennerstrand 1974).
This direct EMG evidence supports a previous conjecture concerning the importance of recruitment in producing conditioned NM responses (Mavritsaki et al. 2007). This study, based on a model constructed earlier by Bartha and Thompson (1992a,b), demonstrated how, in principle, recruitment of motor units in the retractor bulbi muscle could be used to make the NM plant appear more linear to a higher-level controller.
COMMON DRIVE.
A second respect in which the retractor bulbi EMG resembles that of other muscles concerns the relationship between the firing patterns of motor units. First, the firing patterns of different amplitude classes covaried and, second, small spikes always fired at a rate higher than that of large spikes (Figs. 6 and and8).8). This ensemble behavior, sometimes termed an “onion skin” pattern, has been observed in other muscles and interpreted as evidence in favor of the common-drive method of muscle-force control (De Luca et al. 1982a,b, 1996; Jabre and Spellman 1996; Masakado et al. 1995; Sauvage et al. 2006). The key feature of this method is that all motoneurons in the appropriate pool receive the same synaptic drive from a higher-level controller (De Luca and Erim 1994), so that the differences in firing rates between motor units are determined primarily by motor-unit properties (and, as a consequence, are highly regular).
The patterns of spike firing found in the retractor bulbi therefore offer a potentially important clue about the nature of the control signal sent to the motoneuron pool for conditioned NM responses. They are consistent with a single-valued command for globe retraction that is sent to all the motoneurons equally, in contrast to the multiple signals apparently sent to the extraocular muscles for the control of eye position (Dean 1997; Fuchs et al. 1988; Hazel et al. 2002; McCrea et al. 1986). However, the EMG evidence for the common-drive mechanism referred to earlier typically comes from records in which MUAPs can be assigned reliably to individual units, a process involving sophisticated recording and EMG decomposition techniques (e.g., Stashuk 2001). The data here are clearly different in that 1) only spike amplitude is used, a variable highly influenced by extrinsic factors, especially when spikes run together in interference patterns (e.g., Lewicki 1998); and that 2) spikes are assigned to classes not to individual MUs. The fact that the patterns of EMG spike activity found here resemble those obtained when individual MUAPs are identified is thus not by itself conclusive evidence for common-drive control of conditioned eyeblink responses.
There is additional indirect evidence in favor of a common-drive mechanism. The pattern of EMG spike amplitudes found for the retractor bulbi muscle during conditioned responses was also observed for the orbicularis oculi muscle and, when decomposition techniques have been applied to orbicularis EMGs during sustained contractions, they have revealed patterns of MUAPs characteristic of the common-drive mechanism (De Luca et al. 2006). The spike-amplitude patterns observed here are therefore at least consistent with a common-drive command signal. Moreover, overall EMG activity in the orbicularis muscle, which controls the external eyelid, is highly correlated with NM response amplitude during conditioning (Attwell et al. 2002; Lavond et al. 1990; McCormick et al. 1982). It appears that the orbicularis oculi and retractor bulbi pools receive very similar inputs during conditioning, a similarity that could be easily achieved using the simple, single-valued commands required by the common-drive mechanism.
Nonetheless, it was necessary to demonstrate more conclusively whether the patterns of EMG spikes observed in the present study are or are not those produced by a common-drive mechanism. These patterns were therefore modeled explicitly.
Modeling
Comparing the results of data analysis and modeling (Figs. 66–8) shows that the complex patterns of EMG spikes discerned in the data could be well reproduced by a simple common-drive model. An important issue is whether the match between model and data required values of model parameter that were realistic.
PARAMETER VALUES.
To our knowledge, the intrinsic gains and thresholds of accessory abducens motoneurons in the rabbit have not been measured. However, indirect evidence suggests that the model values for these quantities are generally reasonable. For cat abducens motoneurons, gains range from 20 to 40 Hz/nA and thresholds from 1.5 to 8.3 nA (Dean 1997; Grantyn and Grantyn 1978). The corresponding model values are 20 Hz/nA and 0 to 7.5 nA. The main problem appears to be that the model uses motoneurons with lower intrinsic thresholds than are observed experimentally. However, the low values in the model are to some extent uncertain, since the actual time course of the Gaussian motor commands is not measured directly, but instead inferred from the EMG data themselves. Thus the true value of current threshold at which the lowest threshold motoneuron starts firing is not known. Moreover, there is evidence that the red nucleus, thought to be the immediate source of motor commands for conditioned NM responses (Hesslow and Yeo 2002; Robleto and Thompson 2008; Thompson 2005), contains neurons that fire tonically (Desmond and Moore 1991). Given that inputs from the relevant (contralateral) red nucleus are excitatory, tonic drive from the red nucleus could in effect lower the thresholds of accessory abducens neurons. This possibility is considered further in the following text.
A second potential source of concern with regard to model parameter values is the wide range of EMG spike amplitudes. These can be expressed as the ratios of the largest possible spike heights εmax in the model to the threshold spike heights θ0 (Eq. 13), and take the values 27-, 14-, 50-, and 12-fold for data sets RB1–RB4 respectively. These figures appear to be consistent with the estimated range of cat RB unit strengths (i.e., 3.5; Lennerstrand 1974) if motor unit EMG spike amplitudes increase faster than their strengths. There is some evidence from skeletal muscle that suggests this is the case (Goldberg and Derfler 1977). It may be also be the case that the range of motor-unit strengths is larger in rabbit than that in cat: no direct measurements appear to have been carried out for rabbit and, although it seems clear that the cat retractor bulbi muscle is relatively homogeneous (Alvarado et al. 1967; Crandall et al. 1981; Gurahian and Goldberg 1987; Lennerstrand 1974), there is a reference to unpublished observations suggesting that there are three fiber types in the rabbit retractor bulbi (Pachter et al. 1976). For these reasons the wide range of EMG spike amplitudes used in the model appears to be consistent with such data as are currently available.
TONIC INPUTS TO ACCESSORY ABDUCENS REGULATE NM REFLEX EXCITABILITY.
The excitability of accessory abducens neurons, and as a result the threshold and gain of the NM reflex will depend on levels of tonic excitatory and inhibitory inputs. In rabbit the motor command for conditioned NM responses is thought to travel from the anterior interpositus nucleus to the accessory abducens nucleus via the red nucleus (Christian and Thompson 2003; Hesslow and Yeo 2002; Yeo and Hesslow 1998). Consistent with this view, the red nucleus is known to be essential for the production of conditioned NM responses (Bracha et al. 1993; Robleto and Thompson 2008; Rosenfield and Moore 1983), This nucleus provides a significant tonic excitatory input to the accessory abducens that is phasically modulated to create conditioned NM responses (Desmond and Moore 1991).
However, tonic excitatory input from the red nucleus to the accessory abducens nucleus should decrease the thresholds of targeted motoneurons. Consequently inactivation of the red nucleus would be expected to have a secondary effect, namely some reduction in amplitude of unconditioned responses elicited via trigeminal inputs to the accessory abducens. Such amplitude reductions have been reported following inactivations of the red nucleus (Bracha et al. 1993) and of the cerebellar interpositus nucleus, which provides excitatory drive to the red nucleus (Welsh and Harvey 1989). Consistent with this interpretation, elevating red nucleus tonic activity by applying gabazine to interpositus neurons involved in trigeminal reflex blinks increases the amplitude and duration of reflex blinks in the external eyelid system (Chen and Evinger 2006). In evaluating evidence from studies of the mechanisms and sites for generating conditioned responses, the importance of changes in these tonic modulatory influences must be considered. Excitability changes that lead to changes in the expression or performance of the reflex and conditioned responses must be dissociated to reveal mechanisms specific to learning (Hesslow and Yeo 2002; Yeo and Hardiman 1992).
Finally, it is important to note that, in addition to red nucleus inputs and trigeminal inputs, accessory abducens probably receive other modulatory inputs. There is evidence for a diversity of synaptic types in these motoneurons (Destombes et al. 1983) and facial motoneurons receive serotonergic input (Aghajanian and McCall 1980) that may modify their levels of excitability (e.g., Hounsgaard et al. 1986). So, even though conditioned responses may be produced by a common drive, there are probably other inputs to accessory abducens motoneurons that could produce effects not modeled here.
EVIDENCE FROM MOTONEURON FIRING RATES.
The common-drive mechanism for motoneuron recruitment described earlier predicts the firing patterns of motoneurons during CRs, given that the relationship between total EMG-spike firing rate and CR parameters is described by the simple first-order Eq. 15. Above its individual threshold, the firing rate for an individual motoneuron will be approximately linearly related to combined NM position and velocity (right side of Eq. 15).
At present the data required to test this prediction are not available, since the firing rates of motoneurons in the rabbit accessory abducens nucleus during the production of conditioned NM responses have not been recorded systematically (Lepora et al. 2007). Data are available for cat, but it appears that accessory abducens motoneurons in this species do not fire in relation to CRs, at least for the CSs used in the study (Trigo et al. 1999). In the case of orbicularis oculi (OO) motoneurons in the facial nucleus, we have previously shown that their general firing patterns during both URs and CRs are consistent with Eq. 15 (Lepora et al. 2007; their Fig. 13). However, details of individual motoneuron thresholds in relation to eyelid position and velocity have not been presented explicitly. Trigo (1999) showed linear fits to the firing rates of OO motoneurons as a function of eyelid position during CRs (see Fig. 16C in Lepora et al. 2007), but it appears from the fits that the motoneurons would be tonically active when the eyelid is in its resting position. Whether these data are in fact consistent with the common-drive mechanism cannot be determined.
As mentioned earlier, the first-order model derived by Lepora et al. (2007) could be used to fit firing-rate data for OO motoneurons for both CRs and URs. Analysis of the relationship between EMG and UR parameters was not possible in the present study because most of the trials used a shock US with consequent interference with the EMG record from the shock artifact. The interesting question of whether URs are also controlled by a common-drive mechanism, or in some other manner (Bartha and Thompson 1992a,b), remains to be addressed.
Control signal for conditioned eyeblink responses
The results from the present modeling study suggest that, as conjectured previously (Lepora et al. 2007; Mavritsaki et al. 2007), conditioned NM responses are generated by a relatively simple motor command that is converted into the appropriate pattern of motoneuron firing by the recruitment properties of the motoneuron pool. This idea of a simple command is consistent with evidence that very similar signals for conditioned responses are sent to both retractor bulbi and orbicularis oculi muscles (Attwell et al. 2002; Lavond et al. 1990; McCormick et al. 1982). It appears that, despite the complexities of the NM plant, conditioned NM responses are in fact relatively easy to control because the plant has been effectively simplified. The creation of a simplified “virtual plant” probably uses a combination of mechanisms (see discussion in Lepora et al. 2007), in addition to the recruitment process proposed here. The overall result is that a common-drive Gaussian signal sent to the accessory abducens nucleus is straightforwardly related to NM response parameters. In electrophysiological investigations of the mechanisms underlying classical conditioning using the NM response system, this straightforward relationship should simplify parts of the analysis, especially at control sites with relatively direct, oligosynaptic projections to the motoneuron pool.
Although the firing rate properties of eyeblink neurons resemble those of ocular motoneurons to some extent (Trigo et al. 2003), recruitment seems rather different in the two systems. The behavior of ocular motoneurons appears not to result from a common-drive mechanism, but rather from differences in synaptic drives to different motoneurons (Dean 1997; Hazel et al. 2002). The more complex commands required by the oculomotor system probably reflect both its structure, with multiple types of motor unit, and the variety and precision of the responses it controls with respect to both eye velocity and position (e.g., Büttner and Büttner-Ennever 2006). In contrast, the role of the accessory abducens and orbicularis oculi motoneurons is much more straightforward, since their only job is to produce rapid movements of the nictitating membrane and eyelid, followed by a return to their starting position at the end of the blink.
In terms of wider applicability, however, the common-drive mechanism suggested here for conditioned NM responses has been described for numerous other skeletal-muscle responses (e.g., De Luca et al. 1982a,b, 1996; Jabre and Spellman 1996; Masakado et al. 1995; Sauvage et al. 2006). This suggests that a general principle of making movement control as simple for higher-level neural structures may be followed when possible and, indeed, a variety of mechanisms for creating simplified “virtual plants” have been described (e.g., Demer 2006; Mussa-Ivaldi et al. 1994; Nichols and Houk 1976). In this context, the method of analyzing EMG records described here may be convenient for identifying putative common-drive mechanisms for movements not so far investigated.
Role of cerebellum
The role of the cerebellum in the production of CRs has been the subject of extensive investigation and dispute (Hesslow and Yeo 2002; Thompson 2005). It appears to be generally agreed that in the specific case studied here—i.e., delayed conditioning of the NM response in rabbit—the CR is controlled by signals sent to the accessory abducens nucleus from the anterior interpositus nucleus via the red nucleus (Hesslow and Yeo 2002; Thompson 2005). The question then arises of how cerebellar output could deal with the complexities of the NM plant (introduction). The present results suggest a mechanism that would allow the parameters of cerebellar output signals with a Gaussian temporal profile to be related simply to the parameters of the CR (Lepora et al. 2007). This mechanism may also be used for CRs of the eyelids, given the evidence summarized earlier that the same control signal may be applied to both NM and external eyelid CRs in delay conditioning in rabbits. A detailed model of the eyelid plant is needed to indicate how it too might be linearized.
It is important to note, however, that the nature and origin of the control signals used in trace conditioning, which depends on forebrain structures such as the hippocampus, may be significantly different. When considering conditioned eyeblink studies in cats and mice (Delgado-Garcia and Gruart 2006), it is important to remember that there are other strong drives to the orbicularis oculi motoneurons, such as those from the motor cortex. These permit voluntary blinking and are absent in the rabbit NMR system. Nonetheless, it is likely that EMG analysis of the kind used here will be helpful for elucidating whether, in accordance with the general principle of plant simplification outlined earlier, the common-drive mechanism for motoneuron recruitment is used in these other cases of learned eyeblink responses as well.
GRANTS
This research was supported by United Kingdom Biotechnology and Biological Research Council Grant BBS/B/17026 under the Integrative Analysis of Brain and Behaviour Initiative, a Wellcome Trust VIP award to N. Lepora, and National Eye Institute Grant EY-07391 to C. Evinger.
1The online version of this article contains supplemental data.
Glossary
- CS
- conditioned stimulus
- US
- unconditioned stimulus
- CR
- conditioned response
- UR
- unconditioned response
- EMG
- electromyogram
- MN
- motoneuron
- MU
- motor unit
- MUAP
- motor unit action potential
- NM
- nictitating membrane
- RB
- retractor bulbi
- k
- base total-spike rate at which fk becomes >0
- ci
- strength of motor unit i
- f(t)
- total instantaneous EMG spike firing rate at time t
- fi(t)
- instantaneous firing rate of model neuron i at time t
- k(t)
- instantaneous firing rate of spike amplitude class k at time t
- fEMG
- EMG sampling frequency
- gk
- gradient of best-fit line for spike amplitude class k
- Gi
- intrinsic gain of a motoneuron
- hk
- maximum of Gaussian-fit to firing rate for class k
- I(t)
- common drive input current
- I0
- peak input current for a single trial
- Imax
- maximum input current for data set
- k
- spike amplitude class
- l
- gradient of firing rate width σk with EMG amplitude k
- m
- gradient of base firing rate bk with EMG amplitude k
- n
- number of time-bins in fit to spike rate
- nMN
- total number of motoneurons in the pool
- N
- number of EMG spikes divided by width δ of time bin in fit to spike rate
- Nclass
- number of spike amplitude classes
- p
- proportion of sampled motoneurons
- q
- exponent defining the proportion of small to large rheobase currents Ri
- r2
- square of correlation coefficient
- Ri
- rheobase current of a motoneuron
- Rmax
- highest rheobase current in motoneuron pool
- s
- exponent defining the proportion of small to large motor unit strengths c
- t
- time
- tEMG
- times when EMG sampled
- u
- exponent of peak firing rate hk with EMG amplitude k
- v
- exponent of gradient gk with EMG amplitude k
- δ
- width of time bin in fit to spike rate
- εmax
- maximum EMG spike amplitude
- κ
- EMG amplitude in midpoint of range for class k
- θ0
- EMG threshold for spike detection
- θk
- EMG threshold for spike detectionin the spike amplitude class number k
- μ
- mean time of Gaussian-fit to total firing rate f(t)
- σ
- width of Gaussian-fit to total firing rate f(t)
- μk
- mean time of Gaussian-fit to firing rate fk(t) for class k
- σk
- width of Gaussian-fit to firing rate fk(t) for class k
REFERENCES
- Aghajanian GK, McCall RB. Serotonergic synaptic input to facial motoneurons: localization by electron-microscopic autoradiography. Neuroscience 5: 2155–2162, 1980 [Abstract] [Google Scholar]
- Akaboshi K, Masakado Y, Chino N. Quantitative EMG and motor unit recruitment threshold using a concentric needle with quadrifilar electrode. Muscle Nerve 23: 361–367, 2000 [Abstract] [Google Scholar]
- Alvarado J, Steinacker A, Bach-y-Rita P. The ultrastructure of the retractor bulbi muscle of the cat (Abstract). Invest Ophthalmol 6: 548, 1967 [Google Scholar]
- Angelaki DE, Hess BJM. Control of eye orientation: where does the brain's role end and the muscle's begin? Eur J Neurosci 19: 1–10, 2004 [Abstract] [Google Scholar]
- Attwell PJ, Ivarsson M, Millar L, Yeo CH. Cerebellar mechanisms in eyeblink conditioning. Ann NY Acad Sci 978: 79–92, 2002 [Abstract] [Google Scholar]
- Bartha GT, Thompson RF. Control of rabbit nictitating membrane movements: I. A computer model of the retractor bulbi muscle and the associated orbital mechanics. Biol Cybern 68: 135–143, 1992a [Abstract] [Google Scholar]
- Bartha GT, Thompson RF. Control of rabbit nictitating membrane movements: II. Analysis of the relation of motoneuron activity to behavior. Biol Cybern 68: 145–154, 1992b [Abstract] [Google Scholar]
- Binder MD, Heckman CJ, Powers RK. How different afferent inputs control motoneuron discharge and the output of the motoneuron pool. Curr Opin Neurobiol 3: 973–981, 1993 [Abstract] [Google Scholar]
- Bracha V, Stewart SL, Bloedel JR. The temporary inactivation of the red nucleus affects performance of both conditioned and unconditioned nictitating membrane responses in the rabbit. Exp Brain Res 94: 225–236, 1993 [Abstract] [Google Scholar]
- Büttner U, Büttner-Ennever JA. Present concepts of oculomotor organization. Prog Brain Res 151: 1–42, 2006 [Abstract] [Google Scholar]
- Chen FP, Evinger C. Cerebellar modulation of trigeminal reflex blinks: interpositus neurons. J Neurosci 26: 10569–10576, 2006 [Abstract] [Google Scholar]
- Christian KM, Thompson RF. Neural substrates of eyeblink conditioning: acquisition and retention. Learn Mem 11: 427–455, 2003 [Abstract] [Google Scholar]
- Clamann HP, Robinson AJ. A comparison of electromyographic and mechanical fatigue properties in motor units of the cat hindlimb. Brain Res 327: 203–219, 1985 [Abstract] [Google Scholar]
- Crandall WF, Goldberg SJ, Wilson JS, McClung JR. Muscle units divided among retractor bulbi muscle slips and between the lateral rectus and retractor bulbi muscles in cat. Exp Neurol 71: 251–260, 1981 [Abstract] [Google Scholar]
- Dean P. Simulated recruitment of medial rectus motoneurons by abducens internuclear neurons: synaptic specificity vs. intrinsic motoneuron properties. J Neurophysiol 78: 1531–1549, 1997 [Abstract] [Google Scholar]
- Delgado-García JM, Gruart A. Firing activities of identified posterior interpositus nucleus neurons during associative learning in behaving cats. Brain Res Rev 49: 367–376, 2005 [Abstract] [Google Scholar]
- Delgado-García JM, Gruart A. Building new motor responses: eyelid conditioning revisited. Trends Neurosci 29: 330–338, 2006 [Abstract] [Google Scholar]
- De Luca CJ, Adam A, Wotiz R, Gilmore LD, Nawab SH. Decomposition of surface EMG signals. J Neurophysiol 96: 1646–1657, 2006 [Abstract] [Google Scholar]
- De Luca CJ, Erim Z. Common drive of motor units in regulation of muscle force. Trends Neurosci 17: 299–305, 1994 [Abstract] [Google Scholar]
- De Luca CJ, Foley PJ, Erim Z. Motor unit control properties in constant-force isometric contractions. J Neurophysiol 76: 1503–1516, 1996 [Abstract] [Google Scholar]
- De Luca CJ, Lefever RS, McCue MP, Xenakis AP. Behavior of human motor units in different muscles during linearly varying contractions. J Physiol 329: 113–128, 1982a [Abstract] [Google Scholar]
- De Luca CJ, Lefever RS, McCue MP, Xenakis AP. Control scheme governing concurrently active human motor units during voluntary contractions. J Physiol 329: 129–142, 1982b [Abstract] [Google Scholar]
- Demer JL. Current concepts of mechanical and neural factors in ocular motility. Curr Opin Neurol 19: 4–13, 2006 [Europe PMC free article] [Abstract] [Google Scholar]
- Desmond JE, Moore JW. Single-unit activity in red nucleus during the classically-conditioned rabbit nictitating-membrane response. Neurosci Res 10: 260–279, 1991 [Abstract] [Google Scholar]
- Destombes J, Durand J, Gogan P, Gueritaud JP, Horcholle-Bossavit G, Tyc-Dumont S. Ultrastructural and electrophysiological properties of accessory abducens nucleus motoneurones: an intracellular horseradish peroxidase study in the cat. Neuroscience 10: 1317–1332, 1983 [Abstract] [Google Scholar]
- Disterhoft JF, Quinn KJ, Weiss C, Shipley MT. Accessory abducens nucleus and conditioned eye retraction/nictitating membrane extension in the rabbit. J Neurosci 5: 941–950, 1985 [Abstract] [Google Scholar]
- Dudeney JE, Olsen MN, Kehoe EJ. Time-specific extinction and recovery of the rabbit's (Oryctolagus cuniculus) conditioned nictitating membrane response using mixed interstimulus intervals. Behav Neurosci 121: 808–813, 2007 [Abstract] [Google Scholar]
- Durkovic RG. Functional consequences of motor unit recruitment order reversals following spinal cord transection in cat. Somatosens Mot Res 23: 25–35, 2006 [Abstract] [Google Scholar]
- Eglitis I. The glands. In: The Rabbit in Eye Research, edited by Prince JH. Springfield, IL: Charles C Thomas, 1964, p. 38–56 [Google Scholar]
- Ertas M, Stålberg E, Falck B. Can the size principle be detected in conventional EMG recordings? Muscle Nerve 18: 435–439, 1995 [Abstract] [Google Scholar]
- Fuchs AF, Scudder CA, Kaneko CRS. Discharge patterns and recruitment order of identified motoneurons and internuclear neurons in the monkey abducens nucleus. J Neurophysiol 60: 1874–1895, 1988 [Abstract] [Google Scholar]
- Fuglevand AJ, Winter DA, Patla AE. Models of recruitment and rate coding organization in motor-unit pools. J Neurophysiol 70: 2470–2488, 1993 [Abstract] [Google Scholar]
- Goldberg LJ, Derfler B. Relationship among recruitment order, spike amplitude, and twitch tension of single motor units in human masseter muscle. J Neurophysiol 40: 879–890, 1977 [Abstract] [Google Scholar]
- Gormezano I, Kehoe EK, Marshall BS. Twenty years of classical conditioning research with the rabbit. In: Progress in Psychobiology and Physiological Psychology, edited by Sprague JM, Epstein AN. New York: Academic Press, 1983, p. 197–275 [Google Scholar]
- Grantyn R, Grantyn A. Morphological and electrophysiological properties of cat abducens motoneurons. Exp Brain Res 31: 249–274, 1978 [Abstract] [Google Scholar]
- Gruart A, Yeo CH. Cerebellar cortex and eyeblink conditioning: bilateral regulation of conditioned responses. Exp Brain Res 104: 431–448, 1995 [Abstract] [Google Scholar]
- Gurahian SM, Goldberg SJ. Fatigue of lateral rectus and retractor bulbi motor units in cat. Brain Res 415: 281–292, 1987 [Abstract] [Google Scholar]
- Hazel TR, Sklavos SG, Dean P. Estimation of premotor synaptic drives to simulated abducens motoneurons for control of eye position. Exp Brain Res 146: 184–196, 2002 [Abstract] [Google Scholar]
- Heckman CJ, Binder MD. Neural mechanisms underlying the orderly recruitment of motoneurons. In: The Segmental Motor System, edited by Binder MD, Mendell LM. New York: Oxford Univ. Press, 1990, p. 182–204 [Google Scholar]
- Henneman E, Mendell LM. Functional organization of motoneuron pool and its inputs. In: Handbook of Physiology. The Nervous System. Motor Control Bethesda, MD: Am. Physiol. Soc., 1981, sect. 1, vol. II, pt. 1, p. 423–507 [Google Scholar]
- Hesslow G, Yeo CH. The functional anatomy of skeletal conditioning. In: A Neuroscientist's Guide to Classical Conditioning, edited by Moore JW. New York: Springer, 2002, p. 86–146 [Google Scholar]
- Hounsgaard J, Hultborn H, Kiehn O. Transmitter-controlled properties of alpha-motoneurones causing long-lasting motor discharge to brief excitatory inputs. In: Oculomotor and Skeletal Motor System: Differences and Similarities, edited by Freund H-J, Büttner U, Cohen B, Noth J. Amsterdam: Elsevier, 1986, p. 39–49 [Abstract] [Google Scholar]
- Jabre JF, Spellman NT. The demonstration of the size principle in humans using macro electromyography and precision decomposition. Muscle and Nerve 19: 338–341, 1996 [Abstract] [Google Scholar]
- Keenan KG, Farina D, Merletti R, Enoka RM. Influence of motor unit properties on the size of the simulated evoked surface EMG potential. Exp Brain Res 169: 37–49, 2006 [Abstract] [Google Scholar]
- Keenan KG, Valero-Cuevas FJ. Experimentally valid predictions of muscle force and EMG in models of motor-unit function are most sensitive to neural properties. J Neurophysiol 98: 1581–1590, 2007 [Abstract] [Google Scholar]
- Kehoe EJ, Dudeney JE, Ludvig EA, Neufeld J, Sutton RS. Magnitude and timing of nictitating membrane movements during classical conditioning of the rabbit (Oryctolagus cuniculus). Behav Neurosci 122: 471–476, 2008 [Abstract] [Google Scholar]
- Krutki P, Clechanowicz W, Celichowski J, Cywinska-Wasilewska G. Comparative analysis of motor unit action potentials of the medial gastrocnemius muscle in cats and rats. J Electromyogr Kinesiol 18: 732–740, 2008 [Abstract] [Google Scholar]
- Lavond DG, Logan CG, Sohn JH, Garner WDA, Kanzawa SA. Lesions of the cerebellar interpositus nucleus abolish both nictitating membrane and eyelid EMG conditioned responses. Brain Res 514: 238–248, 1990 [Abstract] [Google Scholar]
- Lennerstrand G. Mechanical studies on the retractor bulbi muscle and its motor units in the cat. J Physiol 236: 43–55, 1974 [Abstract] [Google Scholar]
- Lepora NF, Mavritsaki E, Porrill J, Yeo CH, Evinger C, Dean P. Evidence from retractor bulbi EMG for linearised motor control of conditioned nictitating membrane responses. J Neurophysiol 98: 2074–2088, 2007 [Abstract] [Google Scholar]
- Lepora NF, Porrill J, Dean P, Yeo CH, Evinger LC. Simple recruitment model describes control of conditioned nictitating membrane response by retractor bulbi muscle. Program No. 345.316. 2006 Abstract Viewer/Itinerary Planner Washington, DC: Society for Neuroscience, 2006, Online [Google Scholar]
- Lewicki MS. A review of methods for spike sorting: the detection and classification of neural action potentials. Network Comput Neural Syst 9: R53–R78, 1998 [Abstract] [Google Scholar]
- Masakado Y, Akaboshi K, Nagata M, Kimura A, Chino N. Motor unit firing behavior in slow and fast contractions of the first dorsal interosseous muscle of healthy men. Electromyogr Mot Control Electroencephalogr Clin Neurophysiol 97: 290–295, 1995 [Abstract] [Google Scholar]
- Mavritsaki E, Lepora NF, Porrill J, Yeo CH, Dean P. Response linearity determined by recruitment strategy in detailed model of nictitating membrane control. Biol Cybern 96: 39–57, 2007 [Abstract] [Google Scholar]
- McCormick DA, Lavond DG, Thompson RF. Concomitant classical conditioning of the rabbit nictitating membrane and eyelid responses: correlations and implications. Physiol Behav 28: 769–775, 1982 [Abstract] [Google Scholar]
- McCrea RA, Strassman A, Highstein SM. Morphology and physiology of abducens motoneurons and internuclear neurons intracellularly injected with horseradish peroxidase in alert squirrel monkeys. J Comp Neurol 243: 291–308, 1986 [Abstract] [Google Scholar]
- Mussa-Ivaldi FA, Giszter SF, Bizzi E. Linear combinations of primitives in vertebrate motor control. Proc Natl Acad Sci USA 91: 7534–7538, 1994 [Europe PMC free article] [Abstract] [Google Scholar]
- Nandedkar SD, Barkhaus PE, Sanders DB, Stalberg EV. Analysis of amplitude and area of concentric needle EMG motor unit action potentials. Electroencephalogr Clin Neurophysiol 69: 561–567, 1988a [Abstract] [Google Scholar]
- Nandedkar SD, Sanders DB, Stalberg EV, Andreassen S. Simulation of concentric needle EMG motor unit action potentials. Muscle Nerve 11: 151–159, 1988b [Abstract] [Google Scholar]
- Nichols TR, Houk JC. Improvement in linearity and regulation of stiffness that results from actions of stretch reflex. J Neurophysiol 39: 119–142, 1976 [Abstract] [Google Scholar]
- Pachter BR, Davidowitz J, Breinin GM. Morphological fiber types of retractor bulbi muscle in mouse and rat. Invest Ophthalmol 15: 654–657, 1976 [Abstract] [Google Scholar]
- Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical Recipes in C: The Art of Scientific Computing Cambridge, UK: Cambridge Univ. Press, 1992, p. xxvi, 994 [Google Scholar]
- Robleto K, Thompson RF. Extinction of a classically conditioned response: red nucleus and interpositus. J Neurosci 28: 2651–2658, 2008 [Abstract] [Google Scholar]
- Rosenfield ME, Moore JW. Red nucleus lesions disrupt the classically conditioned nictitating response in rabbits. Behav Brain Res 10: 393–398, 1983 [Abstract] [Google Scholar]
- Sanders DB, Stalberg EV, Nandedkar SD. Analysis of the electromyographic interference pattern. J Clin Neurophysiol 13: 385–400, 1996 [Abstract] [Google Scholar]
- Sauvage C, Manto M, Adam A, Roark R, Jissendi P, De Luca CJ. Ordered motor-unit firing behavior in acute cerebellar stroke. J Neurophysiol 96: 2769–2774, 2006 [Abstract] [Google Scholar]
- Stålberg E, Karlsson L. Simulation of the normal concentric needle electromyogram by using a muscle model. Clin Neurophysiol 112: 464–471, 2001 [Abstract] [Google Scholar]
- Stashuk D. EMG signal decomposition: how can it be accomplished and used? J Electromyogr Kinesiol 11: 151–173, 2001 [Abstract] [Google Scholar]
- Thompson RF. In search of memory traces. Ann Rev Psychol 56: 1–23, 2005 [Abstract] [Google Scholar]
- Trigo JA, Gruart A, DelgadoGarcia JM. Role of proprioception in the control of lid position during reflex and conditioned blink responses in the alert behaving cat. Neuroscience 90: 1515–1528, 1999 [Abstract] [Google Scholar]
- Trigo JA, Roa L, Gruart A, Delgado-García JM. A kinetic study of blinking responses in cats. J Physiol 549: 195–205, 2003 [Abstract] [Google Scholar]
- Welsh JP, Harvey JA. Cerebellar lesions and the nictitating membrane reflex: performance deficits of the conditioned and unconditioned response. J Neurosci 9: 299–311, 1989 [Abstract] [Google Scholar]
- Yao WX, Fuglevand AJ, Enoka RM. Motor-unit synchronization increases EMG amplitude and decreases force steadiness of simulated contractions. J Neurophysiol 83: 441–452, 2000 [Abstract] [Google Scholar]
- Yeo CH, Hardiman MJ. Cerebellar cortex and eyeblink conditioning: a reexamination. Exp Brain Res 88: 623–638, 1992 [Abstract] [Google Scholar]
- Yeo CH, Hardiman MJ, Glickstein M. Classical conditioning of the nictitating membrane response of the rabbit. I. Lesions of the cerebellar nuclei. Exp Brain Res 60: 87–98, 1985 [Abstract] [Google Scholar]
- Yeo CH, Hesslow G. Cerebellum and conditioned reflexes. Trends Cogn Sci 2: 322–330, 1998 [Abstract] [Google Scholar]
- Zhou P, Rymer WZ. Factors governing the form of the relation between muscle force and the EMG: a simulation study. J Neurophysiol 92: 2878–2886, 2004 [Abstract] [Google Scholar]
Articles from Journal of Neurophysiology are provided here courtesy of American Physiological Society
Full text links
Read article at publisher's site: https://doi.org/10.1152/jn.00204.2009
Read article for free, from open access legal sources, via Unpaywall:
https://europepmc.org/articles/pmc2775390
Free after 12 months at intl-jn.physiology.org
http://intl-jn.physiology.org/cgi/content/full/102/4/2498
Free to read at intl-jn.physiology.org
http://intl-jn.physiology.org/cgi/content/abstract/102/4/2498
Citations & impact
Impact metrics
Citations of article over time
Smart citations by scite.ai
Explore citation contexts and check if this article has been
supported or disputed.
https://scite.ai/reports/10.1152/jn.00204.2009
Article citations
Bidirectional learning in upbound and downbound microzones of the cerebellum.
Nat Rev Neurosci, 22(2):92-110, 17 Nov 2020
Cited by: 46 articles | PMID: 33203932
Review
Comparison of methods for estimating motor unit firing rate time series from firing times.
J Electromyogr Kinesiol, 31:22-31, 01 Sep 2016
Cited by: 2 articles | PMID: 27623025
Time course of the rabbit's conditioned nictitating membrane movements during acquisition, extinction, and reacquisition.
Learn Mem, 21(11):585-590, 15 Oct 2014
Cited by: 1 article | PMID: 25320350 | PMCID: PMC4201810
Adaptive timing of motor output in the mouse: the role of movement oscillations in eyelid conditioning.
Front Integr Neurosci, 5:72, 30 Nov 2011
Cited by: 32 articles | PMID: 22144951 | PMCID: PMC3226833
Cerebellum and cognition: evidence for the encoding of higher order rules.
Cereb Cortex, 23(6):1433-1443, 22 May 2012
Cited by: 68 articles | PMID: 22617850
Go to all (9) article citations
Data
Data behind the article
This data has been text mined from the article, or deposited into data resources.
BioStudies: supplemental material and supporting data
Similar Articles
To arrive at the top five similar articles we use a word-weighted algorithm to compare words from the Title and Abstract of each citation.
Evidence from retractor bulbi EMG for linearized motor control of conditioned nictitating membrane responses.
J Neurophysiol, 98(4):2074-2088, 05 Jul 2007
Cited by: 9 articles | PMID: 17615135
Discharge profiles of abducens, accessory abducens, and orbicularis oculi motoneurons during reflex and conditioned blinks in alert cats.
J Neurophysiol, 81(4):1666-1684, 01 Apr 1999
Cited by: 50 articles | PMID: 10200203
Localization of retractor bulbi motoneurons in the rabbit.
Brain Res, 226(1-2):93-106, 01 Dec 1981
Cited by: 38 articles | PMID: 7296302
Physiology of the eyelid motor system.
Ann N Y Acad Sci, 1004:1-9, 01 Oct 2003
Cited by: 6 articles | PMID: 14662442
Review
AAEM minimonograph #3: motor unit recruitment.
Muscle Nerve, 14(6):489-502, 01 Jun 1991
Cited by: 18 articles | PMID: 1852155
Review
Funding
Funders who supported this work.
Biotechnology and Biological Sciences Research Council (1)
Cerebellum as a neuronal machine: Behavioural, electrophysiological and computational analysis of classical conditioning
Dean, University of Sheffield
Grant ID: BBS/B/17026
NEI NIH HHS (2)
Grant ID: R01 EY007391
Grant ID: EY-07391