| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
1Department of Psychology, University of Sheffield, Sheffield, United Kingdom; 2School of Psychology, University of Birmingham, Edgbaston, Birmingham; and 3Department of Anatomy and Developmental Biology, University College London, London, United Kingdom; and 4Department of Neurobiology and Behavior, State University of New York, Stony Brook, New York
Submitted 27 February 2007; accepted in final form 3 July 2007
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSREFERENCES |
|---|
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSREFERENCES |
|---|
; Thompson 2005). Nonetheless, the nature of the control signal sent from the cerebellum to the relevant motoneurons is still unclear. Initial data from multiunit recordings in the interpositus nucleus (McCormick and Thompson 1984) indicated a firing rate profile that simply mimicked the profile of the conditioned NMR. Consequently, many models of the cerebellum in eyeblink conditioning essentially stop at the interpositus nucleus, with no representation either of efferent target neural populations such as the red nucleus and motor nuclei, or of the plant itself (Balkenius and Morén 1999; Fiala et al. 1996; Garenne and Chauvet 2004; Gluck et al. 2001; Hofstötter et al. 2002; Medina and Mauk 2000; Moore and Choi 1997). (The term plant refers to "that which is controlled," in this case by the signal sent from the motoneurons, and corresponds to the relevant muscles together with orbital or eyelid tissue). Nevertheless, the assumption that the plant has very simple dynamics is not obviously true and has been questioned explicitly: "it is somewhat unlikely that the actual discharge rate of interpositus neurons could be correlated with the profile of CRs determined from nictitating membrane displacement, as this is a passive, highly damped movement" (Delgado-García and Gruart 2005; p.375). The question is to what extent the system peripheral to the interpositus nucleus presents the cerebellum with complex nonlinear control problems that must be taken into account to produce appropriately timed and shaped conditioned eyeblink responses.
Part of the answer to this question has been provided by a detailed model of a portion of the peripheral system (Bartha and Thompson 1992a,b). This model relates the firing patterns of accessory abducens motoneurons to the NMR, describing the detailed mechanism by which contraction of the retractor bulbi (RB) muscle retracts the globe, displacing Harder's gland and thus forcing the nictitating membrane across the cornea (Eglitis 1964). The model is also consistent with the Lennerstrand (1974) recordings of RB isometric force as a sigmoidal function of input frequency with only a narrow linear range. A subsequent reimplementation of this model (Mavritsaki et al. 2007) has identified two main sources of nonlinearity. The first nonlinearity originates from the sigmoidal pattern in muscle force, which is reproduced in the NMR itself if the system uses rate coding as a control strategy (Fig. 8 of Mavritsaki et al. 2007). In rate coding, all units fire with the same frequency, and force increases occur by increasing frequency. The second major nonlinearity was observed when the system used a simple recruitment control strategy in which units (all of equal strength) either did not fire or fired with a fixed frequency. In this strategy, force is increased by increasing the number of units firing. However, increasing the number of active motor units also increases muscle stiffness, so that the equilibrium position of the muscle in relation to an elastic load increases sublinearly with the number of units recruited (Fig. 9 of Mavritsaki et al. 2007). The presence of these clear nonlinearities suggests that concerns about plant complexity (Delgado-García and Gruart 2005) are entirely justified and seem to preclude any simple relation between control signal and NMR.
|
|
|
|
|
|
|
|
|
); the other was recruitment of motor units of appropriately increasing strength (Fig. 11 of Mavritsaki et al. 2007). These strategies in effect produce a simplified "virtual plant" that is easier to control, because its input signals are restricted to a set over which the plant appears linear, while still maintaining the full response range of the system. Thus, for example, doubling output amplitude can be achieved merely by doubling the amplitude of the input signal, without concern for nonlinear effects such as response saturation. The existence of appropriate recruitment strategies means that the presence of substantive nonlinearities in the output systems peripheral to the interpositus nucleus does not necessarily invalidate the simple models of cerebellar control referred to above. Furthermore, the presentation to the cerebellum of a simple virtual plant in the case of NMR conditioning would be suggestive of a more general operating principle whereby physiological output systems are structured to simplify the task of neural control (Angelaki and Hess 2004; Mussa-Ivaldi et al. 1994; Nichols and Houk 1976).
|
|
; Disterhoft and Weiss 1985; Disterhoft et al. 1985), and one study suggested that such motoneurons fire only in relation to unconditioned NMRs in the cat (Trigo et al. 1999). An early study of multiunit motoneuron firing in rabbit showed activity profiles that strongly resembled, and were highly correlated with, the time-course of the NMR (Fig. 8 of Cegavske et al. 1979), suggesting that "[w]hatever the abducens units do, so will the NM do" (p. 605). Nevertheless, the recordings were from the abducens rather than accessory abducens, only four examples of relevant data were shown, and no explicit demonstration of linearity was undertaken. This study therefore sought to obtain more data on the relation between firing in the accessory abducens motoneuron pool and NMR dynamics by recording the EMG from the RB muscle during the production of conditioned NMRs, and analyzing the EMG records to estimate the occurrence of motor unit action potentials (Gruart et al. 1995; Sanders et al. 1996). Because of the technical difficulties of electrode implantation in this muscle in combination with the recording of NM position during conditioning, these results constitute, to the best of our knowledge, the first systematic description of RB EMG during conditioning. Portions of the findings have been presented previously in abstracts (Lepora et al. 2005; Mavritsaki et al. 2001). |
METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSREFERENCES |
|---|
Data were obtained from two different laboratories: one at University College London (UCL; C. H. Yeo) and the other at the State University of New York at Stony Brook (SUNY; C. Evinger).
UCL procedures
EMG electrodes were implanted in the RB muscle of three male Dutch belted rabbits (2.0–2.2 kg), subsequently referred to as RB1-3. Under gaseous general anesthesia (halothane 2% in O2-N2O) (see Yeo et al. 1985 for details of all surgical and conditioning procedures), a midline incision of the scalp was made, and the soft tissues were reflected laterally to reveal bone. Ophthalmic anesthetic, proxymetacaine hydrochloride 0.5% wt/vol (Opthaine, Squibb), was applied to the corneal surface and to the conjunctiva. After 5–10 min, the scleral margin was held with miniature forceps, and the globe was pulled gently downward and forward. Two to four microwires (Teflon-coated, stainless steel, 50 µm) were inserted into the muscle using a 30-gauge hypodermic needle passing from the scalp opening into the muscle. The wires were separated by 
0.25–1.0 mm. The ends of the wires were stripped of insulation for varying lengths (0–0.5 mm) to enable recording of various numbers of motor units. The final 0.5 mm of each wire was bent back to form a "fish-hook" that retained placement in the muscle as the hypodermic needle was withdrawn. The microwires were passed dorsally through the soft tissues above the orbit toward the scalp incision, where they were terminated with gold-plated crimp connectors that were inserted into a six-channel, nylon pedestal (Plastics 1, Roanoke, VA). Miniature, stainless steel skull screws were placed lateral to the pedestal, and dental acrylic was used to bond the pedestal to the skull and the screws. The scalp was sutured and drawn close to the implanted pedestal. A removable screw cap was fitted to cover the implant, and a nylon monofilament suture loop was placed in the right NM. After 1 wk for postoperative recovery, all subjects underwent one session of adaptation. 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 410 (RB2 and RB3) or 560 ms (RB1), 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. 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 NMR data were recorded for a 2-s period starting 1 s before CS onset. EMG data were sampled at 5–40 kHz and NMR data at 1 kHz. After acquisition, RB1 and RB2 were tested with triplets, in which one trial was CS alone, one was CS + US, and one was US alone. In addition, RB1 was tested with varying CS intensities (65, 75, 85, and 95 dBA).
SUNY procedures
Surgical procedures were identical at SUNY and UCL. To insert EMG electrodes, a pair of Teflon-coated stainless steel wires (76 µm diam bare, 140 µm diam coated, A-M Systems 791000) were inserted into one slip of the RB with a 30-gauge hypodermic needle at SUNY. 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. One rabbit (RB4) was a New Zealand white and the other was a Dutch belted rabbit (RB5).
After initial acquisition of the conditioned response, RB4 was tested under a variety of conditions, including variations in CS and US intensity and extinction. RB5 was tested with variations in CS intensity only.
Data analysis
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 because there was spontaneous movement of the NM just before CS presentation, there was interference in the EMG or NMR signals, or there was evidence from the EMG recording suggesting that the electrode moved during the CR.
EMG analysis
As explained in the INTRODUCTION, the method of analysis adopted here was to estimate the occurrence of motor-unit action potentials from spikes in the EMG records from the intramuscular electrodes (Sanders et al. 1996). The procedure chosen was the simple one of threshold crossing; i.e., whenever the EMG record increased past a threshold, 
, a spike was registered (Fig. 1A). Given an EMG signal x = (x0, x1, x2,...) sampled at frequency fs Hz, [i.e., at discrete times t = (0, 1/fs, 2/fs,...)], this procedure gives a spike sequence, s, that is a binary signal of zeros and ones
 | (1) |
One potential problem with the above procedure is that, at high frequencies, spikes begin to overlap, producing what is termed an interference pattern (Sanders et al. 1996). To estimate the possible effect of such overlap on estimates of spike frequency, we compared the results of 1) summing different EMG records and extracting spikes with 2) extracting spikes from the individual records and summing them. At frequencies where spike interference is significant, the number of extracted spikes from the summed EMGs will be less than the sum of the spikes from the individual EMGs. Our results indicated that the spike frequency at which interference became appreciable was dataset dependent (RB1 = 400 spikes/s, RB2 = 350 spikes/s, RB3 = 400 spikes/s, and RB4 = 550 spikes/s). The dataset dependence appeared to arise from differences between data sets in the width of the EMG spikes. Our observed frequency ranges almost always lay below the interference frequency (peak frequency of RB1 = 450 spikes/s, RB2 = 225 spikes/s, RB3 = 425 spikes/s, and RB4 = 350 spikes/s). There may, however, have been some interference effects for a small number of RB1 and RB3 recordings. However, these records comprise a small proportion of the available data, and no appreciable biasing was observed in our analysis.
A second potential problem with our method of analysis concerns choice of threshold for each subject, which was complicated by variability in the EMG records between trials, possibly the result of electrode movement produced by muscle contraction. The main kinds of variability are fluctuations in general sensitivity (both spike size and background noise changed) or the apparently random appearance of tonic spikes clearly unrelated to the CR, possibly originating in adjacent extraocular muscles. This is a classic signal to noise separation problem, which can only be solved by explicit choice of a criterion for separating signal and noise. Our criterion was an EMG threshold for each subject such that, for all of the trials used in the analysis, there was neither a bias toward trials where EMG spikes occurred without a corresponding CR (threshold too low) nor toward trials where no EMG spikes occurred but a CR did (threshold too high). Presence of a CR depended on whether the NMR passed 0.1-mm amplitude. Further details are given below in Response linearity.
Once spikes had been extracted, their frequencies were calculated from the number of spikes in a 10-ms time interval. The noise levels of individual records usually meant that data had to be pooled across trials based on CR amplitude. Because the initial part of the resulting frequency profiles appeared to resemble a normal distribution (see RESULTS), they were fitted with a Gaussian bell-shaped curve with the equation
 | (2) |
is a measure of its width. These parameters were estimated from the frequency dataset as follows.
1) The mean is the mean of all the times at which spikes occur, given by
 | (3) |
2) Height h is the peak value of the frequency dataset.
3) The width is found by equating the spike total of the frequency dataset to that for the Gaussian curve (i.e., matching areas under the curves). A well-known mathematical identity gives
 | (4) |
= 50 ms is the width of the frequency bin. Fitting was restricted to the interval from CS onset to one ISI after US onset (ISI = interstimulus interval, that is, time between CS and US onsets), to exclude a tail in the frequency profile of some records. This range actually covers most of the frequency variation of each record, because the peak frequency typically occurs significantly before US onset (see Fig. 4 for example).
Response linearity
Response linearity was assessed in two ways. First, the relation between peak CR amplitude and peak EMG frequency was examined. For a given subject, the records from individual trials were sorted on the basis of CR amplitude and averaged over a suitable trial range to reduce the effects of noise (a range of 3 trials gave a good balance between reducing noise and choosing similar records). Each group of records contributed a point to a scatterplot of peak CR amplitude against peak EMG frequency. The best linear and quadratic fits to these points were determined for each subject, using standard statistical techniques. The fits were determined for a range of thresholds for EMG spike extraction. This threshold choice does not depend on fit linearity, because linearity and intersecting a given point are independent measures. The threshold chosen was the one that gave a line of best fit that passed through the origin. This choice was based on the criterion that too high a threshold would give examples of CRs without any EMG spikes and too low a threshold would give EMG spikes without any CRs. In fact the goodness of the linear fit was typically similar over a wide range of threshold values (Fig. 5E), reinforcing the point that this method of choosing the threshold did not bias the analysis in favor of linearity.
The second method of studying response linearity was to determine the best-fit linear filter that related instantaneous EMG frequency to CR amplitude. (The signal-processing term "filter" describes the transformation of a time-varying input signal to a time-varying output signal. In biological or electrical systems, for example, these transformations often involve processes that can be represented by linear differential equations.) Because initial analysis with filters of increasing complexity indicated that, for each subject, the data could be well fitted with a first-order filter, the fitting procedure is described for the first-order case. A first-order filter has a gain (g) and time constant (t0). A unit impulse input to the filter produces an instantaneous response of height g, which decays exponentially with time constant t0. The filter output yi at the ith time step is related both to its output yi-1 at the preceding time step and to the current spike input s by the difference equation
 | (5) |
t where t is the time step). The form of Eq. 5 is well known from signal processing theory, and it can be seen heuristically that it gives the impulse response appropriate for a first-order filter because the bs term drives a response that is a multiplicative gain b of the spike input, while the ayi-1 term makes this response exponentially decay over the following time steps. A system latency between the spike inputs and NM response can also be included in the model by delaying the input signal by an appropriate number of time steps. For a latency 
seconds, the signal s is delayed by d = fs time steps, giving the difference equation
 | (6) |
, g, and t0.
The latency was estimated directly from the data. Examining the start points of the EMG and NM response data showed the latency was typically between 10 and 30 ms for each trial of all datasets with a mean of 20 ms, values that probably reflect the highly damped nature of the plant (Mavritsaki et al. 2007). A sensitivity analysis showed our results did not depend appreciably on the particular value in this observed range, and so the value = 20 ms was used. The parameters g and t0 were estimated by regression analysis, using as an estimate for the fit error the standard RMS difference between the predicted output y from Eq. 6 and the measured output ydata
 | (7) |
, where ej was the fit error of the jth trial. The best fit is found by minimizing the dataset error ê over the parameters g and t0. This optimization is achieved with a standard gradient descent algorithm. The overall quality of the fit can be described with the nonlinear regression coefficient
 | (8) |
Modeling
DETAILED MODEL.
The detailed model of the relation between motoneuron firing and conditioned NMR used in this study was that described in Mavritsaki et al. (2007), which is based closely on a previous model of Bartha and Thompson (1992a,b). The inputs to the model are the firing rates of 100 simulated motoneurons, and its output is instantaneous NM position (Fig. 7). Three processes are represented explicitly: 1) the production of isometric force by a motor unit in the retractor bulbi in response to its input train of motoneuron action potentials; 2) the conversion of this notional force to actual force taking into account the muscle's position and velocity; and 3) the action of this force on the mechanics of globe retraction and Harder's gland compression. In the linearized model used here, motoneurons are recruited as the input signal to the motoneuron pool increases, and in addition, increase their firing rates (details in Mavritsaki et al. 2007). The model was run with a time step of 1 ms, and the output smoothed by taking a 50-ms running average of the simulated isometric force, plotted every 1 ms.
SIMPLE MODEL.
This model is derived from the observation that the relation between EMG spikes and conditioned NMRs can be approximated by a first-order filter. The discrete time version of such a filter is given by Eq. 5 that relates amplitude yi to the spike input si at discrete time ti, given here in more convenient form
 | (9) |
, with
 | (10) |
t = 1/fs 0 gives the differential equation for this continuous time model
 | (11) |
Equation 11 can be interpreted as a notional muscle force F acting on a first-order viscoelastic plant
 | (12) |
 | (13) |
 | (14) |
and 
are averaged over the same time interval T. The relation (Eq. 14) gives an alternative method for qualifying whether the data are consistent with a first-order linear system. Estimates of the time constant t0 can be obtained from the filter fit described above in Response linearity. This specifies F/c from the NM response amplitude y(t). To summarize, the above equations describe a simplified plant model that can be used 1) to predict instantaneous NM displacement from motoneuron firing rates and 2) to estimate notional muscle forces and hence input frequencies from NMR profiles.
|
RESULTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSREFERENCES |
|---|
Data for RB 1–4
RB EMGS AND CONDITIONED NMRS.
Because there is little previous information available on RB EMGs and conditioned NMRs in rabbits, we first describe some basic features of the data. Figure 1A shows an example of spike extraction from the RB EMGs from subject RB4. The full EMG and spike records for that trial, together with the corresponding conditioned NMR, are shown in Fig. 1B. It can be seen that the peak of the CR occurs just before US onset (dotted vertical line, 350 ms after CS onset), whereas the EMG signal, which starts 170 ms after CS onset, is almost finished by time of US onset. This record also shows why selecting an appropriate spike threshold is not straightforward (see METHODS) if some of the spikes that are part of the noise (i.e., unrelated to the CR) are larger than some of the spikes that are part of the signal (i.e., related to the CR). Finally, it is apparent from both panels in Fig. 1 that an extracted spike can correspond to excursions of different amplitude in the original EMG. These amplitude differences are ignored in these analyses, where all spikes are treated as equal, but will be the subject of a subsequent report on possible recruitment of motor units in the RB muscle.
To examine the frequency profiles of the extracted spike records, it proved necessary to pool data across records to reduce the often substantial variability between individual trials. For the examples shown in Fig. 2, data were averaged for each subject over the largest 20% of its CRs and corresponding EMG spike frequency profiles. In all subjects, the CR peak was close to US onset (range, 28 ms before to 23 ms after), with peak EMG frequency 45–75 ms earlier than CR peak. The main differences between the subjects were 1) the values of CR amplitude, which ranged from 2 to 6 mm; 2) the values of EMG peak frequency, which ranged from 200 to 400 spikes/s, possibly reflecting the number of motor units that were sampled by the electrode; and 3) how quickly the NM returned to its resting position. It can be seen that this return, which also appears in the recordings of conditioned NMRs from other laboratories (cf. Fig. 10A), was fastest in RB2 and RB4 and much slower in RB3. The relationship of this response "tail" to EMG activity is not straightforward: there is prolonged EMG activity in the case of subject RB3 but not for RB1, although both subjects have substantial NMR "tails." This lack of relationship may arise because of EMG sampling problems when only few units are active, a view consistent with the irregular appearance of prolonged EMG discharge from trial to trial in subject RB3 (cf. Figs. 2 and 4). Tonic activity in the RB EMGs continuing after the conditioned response has also been observed in a study of eyelid conditioning (Leal-Campanario et al. 2004).
To characterize the signals that control conditioned NMRs, it is important to know how both motoneuron activity and response profile vary with response amplitude. The variation for CRs is shown in Fig. 3. Although in certain cases the traces were still quite noisy even after pooling across sets of three responses, it can be seen that the change in shape with amplitude is relatively slight and much less pronounced than the differences in shape between subjects.
This was also true for EMG spike frequency profiles (Fig. 4), although in this case, the differences between subjects were not as marked as for CR profiles. An important feature of the EMG frequency profiles was that they could be well fitted by Gaussian curves (see METHODS) over the time period starting at CS onset and continuing for twice the duration of the IS1 (time period equals 1 s for RB1 and 0.7 s for RB2-4). Although the timing of the peak of the fitted Gaussians varied relatively little within subjects, there was a slight tendency for larger EMG signals to have earlier peaks. It should be noted that these Gaussian fits are entirely empirical; that is, they do not depend on any prior assumptions of system linearity, and they are significant not only for indicating the actual temporal profile of motoneuronal activity during CRs, but allowing it to be approximated by just three parameters (time and height of peak and width; cf. Figs. 10 and 11).
EMG spikes linearly related to conditioned NMR
Because the inputs to the system as indicated by EMG spike frequency profiles have roughly similar shapes independent of their amplitude, it is possible to exploit the principle of superposition to test for linearity without knowing the exact relation between input and output. If input and output are linearly related, and the temporal profile of the input is invariant, peak output should be directly proportional to peak input. Figure 5 shows this to be a reasonable approximation for these data, with straight line fits giving r2 values of 0.79–0.94. For each subject, the data points were scattered around the straight line in a manner suggesting a noisy rather than a systematically nonlinear relation between the two variables. In fact, fitting a quadratic curve to the data gave very little improvement in the proportion of variance accounted for (range, 0–3% improvement) when both quadratic and linear fits were constrained to pass through the origin. In a further test of robustness, the datasets were fitted with straight lines using different values for the threshold for spike extraction (see METHODS). There was typically a broad range of thresholds that gave good linear fits (Fig. 5E). This is important in indicating that the method for selecting the threshold value, i.e., that value that gave a line of best fit that passed through the origin, did not bias the subsequent analysis toward linearity.
Given the above evidence for overall linearity between input and output, the next step was to search for the linear filter that would give the best fits to the temporal profiles of the responses. As indicated in METHODS, it turned out that a first-order filter was adequate. Filter parameters were estimated for each subject using that subject's entire data set of CR and EMG spike profiles. Typical fits obtained with these parameters are shown in Fig. 6. The example records shown in this figure are averaged over three trials to reduce noise (pooled on the basis of CR amplitude) and show median fit values as determined by R2, the proportion of variance accounted for. Given the noisiness of the data and the fact that the EMG spikes were a only a sample of total motoneuron activity, the fits shown in these panels (R2 = 0.90–0.97) suggest that a first-order filter is a reasonable approximation of the relationship between EMG spike input and conditioned NMR output.
To check for overfitting in our analysis, we repeated the above procedure with separate fitting and validation sets. Filter fits were determined on the central 50% of trials (ordered by NMR peak amplitude) and goodness of fit determined on the outer 50% of trials. Parameter and fit values (R2 = 0.90–0.98) were consistent with the above method.
Modeling
A detailed model of the relation between the firing rates of motoneurons in the accessory abducens nucleus and the displacement of the NM (Fig. 7) has been described by Bartha and Thompson (1992a,b). This model includes a number of features that result in that relation being nonlinear. The first issue addressed in this section therefore is how these results could be at all compatible with the detailed model. The second issue is the extent to which a much simpler linear model, based on the linear filter described above, can account both for the present data and for other published data on the temporal profiles of conditioned NMRs.
DETAILED MODEL.
Bartha and Thompson's detailed model is organized into three parts (Fig. 7), with the first describing the production of isometric force by an individual motor unit, the second deriving the actual force from isometric force using the muscle's length and velocity, and the final part calculating the instantaneous position of the NM as a function of the actual force. A subsequent implementation of this model (Mavritsaki et al. 2007) showed that its behavior depended critically on the nature of the control signals sent by the motoneuron pool. For simple rate-coding, in which muscle force is determined only by variations in the frequency with which all the units fire, there is a nonlinear, sigmoidal, relationship between firing frequency and NMR amplitude (Fig. 8 of Mavritsaki et al. 2007 ), which reflects the relationship between the firing frequency and isometric force of a single motor unit (Lennerstrand 1974).
However, it was possible to show that if motoneuron recruitment were combined in a suitable manner with frequency modulation, the model behaved in a much more linear fashion (Fig. 10 of Mavritsaki et al. 2007). The response of this linearized version of the model to Gaussian-shaped frequency profiles is shown in Fig. 8A. It can be seen that peak NMR amplitude is linearly related to peak input frequency. The best-fit linear filter for this version of the model can be estimated from input-output data with the same techniques as those used for the actual data (METHODS), and an example of filter performance is shown in Fig. 8B. The CR output of the detailed model (labeled "biophysical model" in Fig. 8B) is approximated by the output of the linear filter (with time constant of 0.1 s and gain of 0.12). Thus these data showing a simple relation between EMG spikes and CR profiles can be consistent with a detailed model of the RB muscle and associated plant. It is important to note that, although recruitment is one way of improving the model's linearity, it is not necessarily the only way, so the results shown in Fig. 8 do not prove that the real system does use recruitment. This point is addressed further in the DISCUSSION.
LINEAR MODEL.
As indicated in Fig. 6, the fact that the relationship between EMG spikes and conditioned NMRs can be approximated by a linear filter suggests a simple linear model of CR production, in which the summed motoneuron frequency profiles signal desired muscle force, and this force acts on a first-order plant (Fig. 7, bottom row). A simple test of such a model is to run the temporal profiles of a CR back through the first-order plant to estimate the force. In practice, this was achieved by using the time constant t0 of the fitted filter (values given in the legend to Fig. 6) to construct the force F 
+ y/t0 from Eq. 12. This estimated force can be compared with the frequency profile of the associated EMG to verify they are linearly related, as assumed by Eq. 14.
For these datasets, peak EMG frequencies and estimated peak forces were linearly related (data not shown: r2 range, 0.80–0.94). Figure 9 shows records with median goodness of fits for data pooled over three trials. The reasonable fits between the actual EMG and inferred force profiles are consistent with the simplified model in which muscle input and force are related by a gain term only. It is of interest that this simple model gives a relationship between the firing rate FR of the motoneuron pool and NM position y and velocity dy/dt
 | (15) |
shown in A have been run back through a first-order model to generate the input profiles shown in B. Because only NMR data are available, the actual time constant of the best-fit first-order filter cannot be determined. We instead used a value of 0.16 s, which is the mean of the values found for RB1-4 (see legend to Fig. 6). The resulting input profiles could be well fitted with Gaussian curves, as could the EMG spike-frequency profiles in this study (Fig. 4). Both the peaks and widths of the Gaussian curves shown in Fig. 10 vary with ISI (interval between CS and US onset in paired trials: actual data in Fig. 10 from CS alone trials). The latency of the peak increases with ISI, as does the width of the curve. In addition, the height of the peak diminishes as the ISI increases.
More details concerning the first two of these relationships are shown in Fig. 11, for both the data from Smith (1968) and data from four additional studies (Coleman and Gormezano 1971; Millenson et al. 1977; Welsh 1992; Yeo et al. 1997). For each study, the conditioned NMR profiles were run backward through a range of simple models with time constants from 0.07 to 0.23 s (as in RB1-4) and the predicted motoneuron signals fitted with Gaussian curves. The time of the peak of each Gaussian (relative to CS onset) is plotted against ISI in Fig. 11A, with error bars from the spread of time constants. The scatter plot is reasonably well fitted by a straight line, whose equation indicates that the time to peak of summed motoneuron activity is approximately the same as the ISI. This would be consistent with a wealth of data concerning CR timing. The time to peak of the inferred activity is also linearly related to its width (Fig. 11B). The two equations for time to peak (Fig. 11) can be combined to yield an expression that relates the onset latency of the activity to the ISI
 | (16) |
denotes the number of SDs before the mean at which the activity is deemed to start, so that µ = /2+ . A starting criterion of 5% of peak amplitude corresponds to 2.5, which gives
 | (17) |
15 ms and a second related to the ISI (about one half its value). As such, it is broadly consistent with measurements summarized by Gormezano et al. (1983), indicating that "the frequency distribution of CR onset latencies is centered at about the midpoint of the CS–US interval" (p. 216). It is of interest that, although some data shown in Fig. 11 were obtained with a counterweighted transducer (Coleman and Gormezano 1971; Millenson et al. 1977; Smith 1968) and some with a noncounterweighted transducer (Welsh 1992; Yeo et al. 1997), there seems to be no obvious difference between them. Additional observations
DATA ANALYSIS FOR SUBJECT RB5.
Analysis of EMG spikes for subject RB5 indicated that their peak frequency measured over 50 ms (60 spikes/s; Fig. 12E) was substantially lower than the peak frequencies found for RB1-4 (200–400 spikes/s shown in Fig. 2) and that individual spikes all had similar temporal profiles as indicated by principal component analysis (Stashuk 2001). It was therefore possible that, in this particular subject, the EMG electrode was recording primarily from a single motor unit rather than from a number of units that would allow a reasonable estimate of motoneuron pool activity. Moreover, the CRs in RB5 tended to be either small (2 mm) or large (8 mm), making it difficult to relate EMG parameters to CR parameters.
|
50 ms). Inspection of the data suggests that these short ISIs correspond to brief bursts (2–4 spikes) that occur 20% of the times that a single spike would be expected from a Poisson process. Whether these bursts arise from noisy inputs to a single motoneuron or from multiple motoneurons remains to be determined. Spectral analysis of NMR
Although this study is concerned with conditioned responses of the NM, the claim that external eyelid responses are driven by a neural oscillator, running at 8 Hz in rabbit (Gruart et al. 2000) and 20 Hz in cat (Domingo et al. 1997), is of relevance if both conditioned NM and external eyelid responses are driven by a common cerebellar command. We therefore studied whether there is evidence for similar input oscillations in these data. Because the spike trains extracted from the RB EMG records were noisy, we used instead the inferred force records generated by Eq. 12 and shown in Fig. 9. Power spectra were determined from the fast Fourier transforms of the inferred force records of subjects RB1-5 and were found to contain no consistent peaks across subjects apart from that at 0 Hz to be expected from Gaussian input (data not shown). Above 10 Hz, the spectra indicated the presence of high-frequency noise in amounts that varied from subject to subject.
Analyses of eyelid records had assumed force was simply proportional to acceleration (Domingo et al. 1997; Gruart et al. 2000) and so had examined the power spectra of acceleration profiles. We therefore also checked acceleration profiles from these data (even though for a first-order model, acceleration cannot be simply related to force input) to see whether evidence for oscillations could be found. As might be expected, differentiating the position traces greatly increased the effects of high-frequency noise (cf. Fig. 1 of Welsh 1992) >10 Hz, but no consistent peaks were found below 10 Hz.
It appears therefore that these data provide no evidence for a neural oscillator in the control of the rabbit's NMR. There are a number of possible reasons for this apparent discrepancy with the eyelid data. One is that the oscillations actually arise from particular mechanical properties of the eyelid (plus recording coil) itself, which are not shared with the NMR. Another is that the actual NM mechanics, as represented by first-order plant with a time constant of 0.1 s, would attenuate an 8-Hz component in the input by a factor of 6 in comparison with low-frequency (<0.5 Hz) components, so that its presence would be hard to detect. Third, an analysis of conditioned eyelid responses in rabbits after section of the facial nerve suggested that "the oscillatory component recorded from eyelid responses comes mainly from the facial motor system" (Leal-Campanario et al. 2004; p. 1552), i.e., not from the RB component.
|
DISCUSSION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSREFERENCES |
|---|
). Application of the first-order model to data from a range of previous studies suggested that these CRs also were produced by summed motoneuron firing with a Gaussian profile, whose width and peak latency varied with interstimulus interval. These results strongly suggest that the cerebellum deals with a virtual plant that substantially simplifies control of CR amplitude and timing in NMR conditioning. Methodological issues
The difficulties of extracting motor unit action potential (MUAP) trains from the EMG are well known (Christodoulou and Pattichis 1999; Farina et al. 2001; Lewicki 1998; Loudon et al. 1992; McGill and Dorfman 1985; Stashuk 2001). Nevertheless, we avoided one of the major problems, namely the variations in MUAP shape and amplitude that make it difficult to classify MUAPs according to the motor unit of origin, because our objective was to estimate total MUAP activity, not the firing characteristics of individual motor units. The main methodological problem for this study was whether simple spike extraction from the EMG provided a good estimate of the number of summed MUAPs. Experimental studies have indicated that the soleus EMG can indeed behave as the simple algebraic sum of MUAPs (Day and Hulliger 2001) at least up until a certain level of muscle activity. We attempted to assess whether this was also true for the EMG records produced by our RB electrodes and concluded that the maximum spike rates obtained were still within the linear range. Moreover, insofar as synchronization would be a source of bias in estimating total MUAP activity (Day and Hulliger 2001), the evidence of substantial trial to trial variability in spike timing from a putative single motor unit (Fig. 12) suggests it is unlikely to be a major concern when the motor response of interest is the conditioned NMR.
Relation to previous data
The conditioned NMRs recorded here are similar to those obtained in previous studies (Gormezano et al. 1983), with respect to both temporal profile (including the long "tail" denoting a slow return to resting position lasting 1 s) and maximum amplitude of 2–8 mm (Garcia et al. 2003; Millenson et al. 1977; Welsh 1992; Yeo et al. 1997). The relatively modest change in the shape of the CR profile with increases in CR amplitude is also consistent with previous reports of conditioned NMR and eyelid profiles. These indicate that during acquisition there are 1) slight reductions in onset latency (Gormezano et al. 1983; Gruart et al. 1995; Schneiderman 1966; Schneiderman and Gormezano 1964; Smith 1968; Smith et al. 1969), although there seems to be great variability between the latencies of individual CRs early in acquisition (Frey and Ross 1968; Garcia et al. 2003; Schneiderman 1966; Schneiderman and Gormezano 1964; Smith et al. 1969) and 2) no consistent change in peak latency (Gruart et al. 1995, 2000; Leal-Campanario et al. 2004; Smith 1968; Welsh 1992). Thus in general terms, the conditioned NMRs recorded here are representative, and this is confirmed by the finding that a model derived from them can also be applied to data from other laboratories (Fig. 11).
As far as relations between CR profiles and EMG activity in the RB muscle are concerned, there seem to be very few previous data, possibly because of the technical difficulties of electrode implantation in this muscle in combination with the recording of NM position during conditioning. Thus Berthier (1992) shows only a single record of NM position and RB EMGs. These results constitute, to the best of our knowledge, the first systematic description of RB EMG during conditioning.
There have, however, been previous reports that indirectly support there being a close and possibly linear relation between summed RB motor unit firing and CR temporal profiles. First, activity in the EMG of the orbicularis oculi (OO) muscle is highly correlated with conditioned NMR amplitude (Lavond et al. 1990; McCormick et al. 1982) to such an extent that some studies use OO EMG activity as a measure of conditioning in its own right (Ivarsson and Svensson 2000; Mauk and Thompson 1987; Skelton 1988). Second, there have been a number of reports suggesting linear relationships between aspects of OO EMG activity and the size and velocity of unconditioned blink responses (Evinger et al. 1991; Gruart et al. 1995; Manning and Evinger 1986; VanderWerf et al. 2003). Finally, single unit recordings of activity in the motoneurons that innervate the OO muscle have suggested a linear relationship between their firing rate and the position of the eyelid for CRs (Fig. 16 of Trigo et al. 1999), and between firing rate and lid velocity for URs (Fig. 5 of Trigo et al. 1999). It may be recalled that for a first-order linear plant, input is related to a combination of position and velocity of output (Eq. 15). However, for slow movements the position term will dominate the relationship, and for fast movements the velocity term will predominate. In the data reported by Trigo et al. (1999) the URs were rapid (peak velocities 
700°/s), whereas the CRs, produced by trace conditioning with a brief airpuff CS, were almost an order of magnitude slower (peak velocities, 80°/s). Simulations show that our first-order model can produce the relationships between input signal and lid movement for both the URs and CRs shown in Figs. 5 and 16 of Trigo et al. (1999) (results shown in Fig. 13).
|
