Recruitment in retractor bulbi muscle during eyeblink conditioning: EMG analysis and common-drive model

Abstract

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.

Fig. 1.

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.)

Fig. 2.

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.

Fig. 3.

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.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.

Fig. 4.

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).

Fig. 5.

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 θ₀ 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).

Fig. 6.

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 I₀ 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.

Fig. 7.

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 r² value (where r is the correlation coefficient) are shown on B and E. The best-fit exponential and its r² value are shown on C and F. The ratio of peak firing rates h_k/h₁ 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).

Fig. 8.

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 θ₀. 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).