Abstract
We assessed here true causal directionalities in cerebellar–motoneuron (MN) network associations during the classical conditioning of eyelid responses. For this, the firing activities of identified facial MNs and cerebellar interpositus (IP) nucleus neurons were recorded during the acquisition of this type of associative learning in alert behaving cats. Simultaneously, the eyelid conditioned response (CR) and the EMG activity of the orbicularis oculi (OO) muscle were recorded. Nonlinear association analysis and time-dependent causality method allowed us to determine the asymmetry, time delays, direction in coupling, and functional interdependences between neuronal recordings and learned motor responses. We concluded that the functional nonlinear association between the IP neurons and OO muscle activities was bidirectional and asymmetric, and the time delays in the two directions of coupling always lagged the start of the CR. Additionally, the strength of coupling depended inversely on the level of expression of eyeblink CRs, whereas causal inferences were significantly dependent on the phase information status. In contrast, the functional association between OO MNs and OO muscle activities was unidirectional and quasisymmetric, and the time delays in coupling were always of opposed signs. Moreover, information transfer in cerebellar–MN network associations during the learning process required a “driving common source” that induced the mere “modulating coupling” of the IP nucleus with the final common pathway for the eyelid motor system. Thus, it can be proposed that the cerebellum is always looking back and reevaluating its own function, using the information acquired in the process, to play a modulating-reinforcing role in motor learning.
Introduction
The analysis of functional relationships between motor activities (both actual movements and muscle electromyography), neural firing rates, and oscillatory rhythms at different motor and premotor levels has provided important insights regarding how neural circuits generate appropriate behaviors (Varela et al., 2001). In particular, the information collected has led to a better understanding of the physiological mechanisms underlying the role of cerebellar centers in neural control of movement (Eccles et al., 1967; Thach et al., 1992; Welsh and Llinás, 1997; Lang et al., 2006; Gerwig et al., 2007; Cheron et al., 2008). At the same time, recent analytical procedures have enabled the determination of how deep brain stimulation suppresses synchronized firing activities to improve its therapeutic efficiency in certain motor disorders (Llinás et al., 1999; Titcombe et al., 2001) by using different time-delayed feedback strategies (Rusemblum and Pikovski, 2004) and phase-resetting techniques (Tass, 2003). A widely used method for estimating functional interdependences between two physiological signals is the linear cross-correlation function. This index measures the strength (strong or weak) of association, the time delay, and the relative linearity between two processes. However, the cross-correlation function is unable to distinguish between a closed- and an open-loop system. Similar significant linear cross-correlation estimations may appear in systems both with and without feedback. Thus, these analytical procedures do not help to elucidate any causal relationship within the physiological systems involved. Therefore, to fully understand information processing at different levels of the motor system and—in particular—the functional interdependence between signals recorded during actual motor learning processes, a directional analysis to reveal causal influences is essential. The nonlinear association method (Lopes da Silva et al., 1989; Pijn et al., 1997; Kalitzin et al., 2007) and Granger time-dependent causality analysis (Granger, 1980; Geweke, 1982; Nolte et al., 2008) have proven useful in these experimental situations (Bernasconi et al., 2000; Kaminski and Liang, 2005; Witte et al., 2009).
In the present study, we have carried out a nonlinear association analysis to quantify the time-dependent coupling between simultaneously recorded neural [neuronal responses (NRs) of identified facial motoneurons (MNs) or deep cerebellar posterior IP neurons] and muscular [EMG activity of the orbicularis oculi (OO) muscle] signals. We had already preprocessed crude NRs and EMG activities, validated our hierarchical cluster analysis for averaged blocks of trials and experimental sessions and used a cumulative neuronal integration method for the phase information (Sánchez-Campusano et al., 2007). This information allowed us to develop multivariate autoregressive models of physiological time series during the motor learning processes and to determine the influence of the phase-inversion properties of interpositus (IP) neurons on the acquisition of new motor abilities (Sánchez-Campusano et al., 2007). Here, we have used time-dependent Granger causality analysis (Granger, 1980; Geweke, 1982) to reveal any possible causal relationship between time series resulting from neuronal firing rates collected from identified OO MNs or from posterior IP neurons and learned motor responses. The aim was to determine the causal involvement of the cerebellum in the generation of new motor abilities, using as an experimental model the classical conditioning of eyelid responses in alert cats.
Materials and Methods
Experimental animals.
Experiments were carried out on eight female adult cats (weighing 2.3–3.2 kg) obtained from an authorized supplier (Iffa-Credo). Experiments were carried out in accordance with the guidelines of the European Union (86/609/EU, 2003/65/EU) and Spanish regulations (BOE 252/34367-91, 2005) for the use of laboratory animals in chronic studies.
Surgery.
Animals were anesthetized with sodium pentobarbital (35 mg/kg, i.p.) following a protective injection of atropine sulfate (0.5 mg/kg, i.m.) to prevent unwanted vagal responses. A search coil (five turns, 3 mm in diameter) was implanted into the center of the left upper eyelid at ∼2 mm from the lid margin (Fig. 1A). The coil was made from Teflon-coated multistranded stainless steel wire (50 μm of external diameter). Coils weighed ∼1.5% of the cat's upper lid weight and did not impair eyelid responses. Animals were also implanted in the ipsilateral OO muscle with bipolar hook electrodes aimed for EMG recordings. These electrodes were made of the same wire as the coils and bared 1 mm at their tips.
Schematic representation of the experimental design and of recorded physiological signals. A, Diagram illustrating the stimulating (Stim.) and recording (Rec.) sites, as well as the eyelid coil and EMG electrodes implanted in the upper eyelid. Kinetic neuronal commands were computed from the firing activities of antidromically identified OO MNs located in the facial nucleus (VII n) and from neurons located in the ipsilateral cerebellar posterior IP nucleus (IP n). R n, red nucleus; V n, trigeminal nucleus. B, Diagrammatic representation of the experimental series (S1 for MN and S2 for IP neuron recordings, both obtained in simultaneity with the EMG activities of the OO muscle and eyelid position recordings) during classical eyeblink conditioning. C, For classical conditioning of eyelid responses we used a delay paradigm consisting of a tone as a CS. The CS started before but coterminated with an air puff used as an US. D–K, A set of recordings collected from the 10th conditioning session from two representative animals. D, Raster representation of action potentials generated by two identified neurons (an IP and an MN) across one learning block (n = 9 trials). Each row in the raster represents a single trial. The two single trials marked with dotted rectangles in D correspond to the global representation of collected physiological signals illustrated from E to K. E, Conditioning paradigm and firing activities of two (IP and MN) selected neurons during a single trial. Action potentials (IP spikes), marked with blue plus signs (D and E), correspond to the direct representation of the neuronal activity in the IP n (IP raw recordings, in F), and its respective instantaneous frequency (IP firing rate, in G). Action potentials (MN spikes) recorded from an OO MN are indicated with magenta plus signs (D and E) and are the direct representation of the neuronal activity in the facial nucleus (MN raw recordings, in H) and its corresponding instantaneous frequency (MN firing rate, in I). J, K, These traces illustrate the EMG activity of the OO muscle (OO EMG, in J), and the direct recording of the eyelid position by the magnetic field search-coil technique (in K). The beginning of the CR was taken as the zero reference point (in K) for subsequent analysis.
Four of the animals were prepared for the chronic recording of antidromically identified facial MNs projecting to the OO muscle. For this, two stainless steel hook electrodes were aimed at the zygomatic subdivision of the left facial nerve, 1–2 mm posterior to the external canthus. The other four animals were prepared for the chronic recording of antidromically identified left posterior IP neurons. In this case, a bipolar-stimulating electrode, made of 200 μm enamel-coated silver wire, was implanted in the magnocellular division of the right (contralateral) red nucleus following stereotaxic coordinates (Berman, 1968). A recording window (5 × 5 mm) was opened in the occipital bone of all of the animals to allow access to the facial or the IP nuclei. The dura mater was removed, and an acrylic chamber was constructed around the window. The cerebellar surface was protected with a piece of silicone sheet and sterile gauze and hermetically closed using a plastic cap. Finally, animals were provided with a head-holding system for stability and proper references of coil and recording systems. All the implanted electrodes were soldered to a socket fixed to the holding system. A detailed description of this chronic preparation can be found elsewhere (Trigo et al., 1999; Gruart et al., 2000; Jiménez-Díaz et al., 2004; Sánchez-Campusano et al., 2007).
Electrical recordings and stimulating procedures.
Eyelid movements were recorded with the magnetic field search-coil technique (Gruart et al., 1995). The gain of the recording system was set at 1 V = 10°. The EMG activity of the OO muscle was recorded with differential amplifiers at a bandwidth of 0.1 Hz to 10 kHz. Action potentials were recorded in facial and IP nuclei with glass micropipettes filled with 2 m NaCl (3–5 MΩ of resistance) using a NEX-1 preamplifier (Biomedical Engineering). For the antidromic activation of recorded neurons, we used single or double (interval of 1–2 ms) cathodal square pulses (50 μs) with current intensities <300 μA. Only antidromically identified OO MNs and IP neurons were stored and analyzed in this study (Fig. 1A,B). Site location and identification procedures have been described in detail for facial MNs (Trigo et al., 1999) and posterior IP neurons (Gruart et al., 2000; Jiménez-Díaz et al., 2004; Sánchez-Campusano et al., 2007).
Classical conditioning.
Classical eyeblink conditioning was achieved by the use of a delay-conditioning paradigm (Fig. 1C). A tone (370 ms, 600 Hz, 90 dB) was used as conditioned stimulus (CS). The tone was followed 270 ms from its onset by an air puff (100 ms, 3 kg/cm2) directed at the left cornea as a unconditioned stimulus (US). Thus, the tone and the air puff terminated simultaneously. Tones were applied from a loudspeaker located 80 cm below the animal's head. Air puffs were applied through the opening of a plastic pipette (3 mm in diameter) located 1 cm away from the left cornea.
Each animal followed a sequence of 2 habituation, 10 conditioning sessions, and 3 extinction sessions. A conditioning session always consisted of 12 blocks separated by a variable (5 ± 1 min) interval. Each block consisted of 10 trials separated by intervals of 30 ± 10 s. Within each block, the CS was presented alone during the first trial—i.e., it was not followed by the US. A complete conditioning session lasted for ∼2 h. The CS was presented alone during habituation and extinction sessions for the same number of blocks per session and trials per block and with similar random interblock and intertrial distributions (Gruart et al., 1995).
Histology.
At the end of the recording sessions, animals were deeply reanesthetized (50 mg/kg sodium pentobarbital, i.p.). Electrolytic marks were placed in selected recording sites with a tungsten electrode (1 mA for 30 s). Animals were perfused transcardially with saline and phosphate-buffered formalin. Serial sections (50 μm), including the cerebellum and the brainstem, were mounted on glass slides and stained with toluidine blue or cresyl violet for confirmation of the recording sites (Gruart et al., 2000; Jiménez-Díaz et al., 2004).
Data collection and multiparametric statistical analysis.
Eyelid position, EMG and neuronal activity, and rectangular pulses corresponding to CS and US presentations were stored digitally on a computer, using an analog–digital converter (CED 1401 Plus; Ceta Electronic Design). Commercial computer programs (spike 2 and SIGAVG; Ceta Electronic Design) were used for acquisition and on-line conventional analysis. The off-line analysis and quantification and representation programs used here were developed by one of us (R.S.-C.) with the help of MATLAB routines.
Action potentials recorded from identified OO MNs and IP neurons were computed and quantified. The quantification algorithm took into account the identification of the action potential's standard waveform and the absolute refractory (mean duration of the action potential) and relative refractory (minimum interspike time interval) periods (Sánchez-Campusano et al., 2007). The instantaneous firing rate was calculated as the inverse of the interspike intervals. Velocity and acceleration profiles were computed digitally as the first and second derivatives of eyelid position records after low-pass filtering of the data (−3 dB cutoff at 50 Hz and zero gain at ∼100 Hz) (Domingo et al., 1997).
Computed results were processed for statistical analysis using the Statistics MATLAB Toolbox. As statistical inference procedures, both ANOVA (estimate of variance both within-groups and between-groups, on the basis of one dependent measure) and MANOVA (estimate of variance in multiple-dependent parameters across groups) were used to assess the statistical significance of differences between groups. The corresponding statistical significance test [i.e., F[(m−1),(m−1)×(n−1),(l−m)] statistic] was performed, with sessions as repeated measures, coupled with contrast analysis when appropriate (Hair et al., 1998; Grafen and Hails, 2002). The orders m (number of groups), n (number of cats), and l (number of multivariate observations) were reported accompanying the F statistic values. The corresponding statistical significance tests (i.e., Student's t test and F statistic) were performed for the parameters of nonlinear correlation analysis and causal inference method (see below). These parameters were statistically significant if the resulting probability p met the requirements of the hypothesis test: p < 0.05, in which 0.05 was the predetermined statistical significance level. This hypothesis test is done by using modified Fisher's z-transformation (w) to associate each measured nonlinear association index (η) with a corresponding w-transformation. Then, each w is approximately normally distributed with values in [−inf +inf]. If the probability p is less than the predetermined significance level (0.05), then reject the null hypothesis (H0, no correlation or no difference between population means) (Belsley et al., 1980). A small value of p indicates that the two distributions (w-transformation and Gaussian) are significantly correlated (see supplemental Appendix S1, available at www.jneurosci.org as supplemental material). In all of the cases, if H0 = 0, we do not reject null hypothesis at significance level of 0.05. If H0 = 1, we reject null hypothesis at significance level of 0.05.
Collected data were analyzed using a multivariate cluster technique. Multivariate cluster analysis is an optimal tool to link the kinetic neural commands (i.e., motor and premotor neuronal activities computed from different quantitative parameters) and kinematic (i.e., motor activities computed from EMG and actual movement) parameters during motor learning; used here, it enabled us to determine the intrinsic coherence of collected data (Sánchez-Campusano et al., 2007). For this, we developed the necessary computer programs and algorithmic procedures to deal with such a huge amount of data (60 parameters quantified across 15 experimental sessions collected from eight experimental animals). Importantly, this analysis allowed neural firing properties recorded from different animals during the same conditioning paradigm to be correctly assigned to the corresponding experimental session [i.e., in agreement with the actual conditioned response (CR) evoked during the CS–US interval].
An additional advantage of this experimental approach is that once collected data were properly arranged according to the hierarchical cluster tree (and that nonsignificant data were rejected), it was possible to determine analytically the multiple and coherent evolution of kinetic and kinematic parameters, the dynamic nonlinear association functions relating MNs and IP NRs to EMG activity, the causal inferences between the firing neuronal time series and learning motor responses time series, and—finally—the relationship between the causal inference and phase-inversion properties of OO MNs and IP neurons with regard to acquired CRs (Sánchez-Campusano et al., 2007).
Multivariate analyses of physiological time series.
Multivariate analysis is extensively used with the aim of studying the relationship between simultaneously recorded signals or their equivalent time series. Multivariate time series tools [(1) nonlinear association analysis and (2) causal inferences using time-dependent causality analysis] enabled us to determine the functional relatedness, asymmetry, time delay, direction in coupling, and causal inferences between physiological recordings [i.e., neuronal activities generated in facial and cerebellar IP nuclei (Fig. 1A, VII n or IP n), and learned motor responses (OO EMG activity or conditioned eyelid responses)] collected during classical conditioning sessions. In practice, we illustrated the use of multivariate analyses for the assessment of the strength (strong or weak), functional nature (unidirectional or bidirectional), and the type (linear or nonlinear) of interdependences between these physiological time series. The analytical procedures and multivariate analyses of physiological signals developed here were designed by one of us (R.S.-C.) with the help of MATLAB routines.
(1) Nonlinear association between EMG activity and neuronal recordings.
We used nonlinear correlation analysis to investigate the functional association between the EMG activity [YEMGOO(t), for the OO muscle] and the neuronal recording [XNR(t), with XNRMN(t) for facial motoneuronal or with XNRIP(t) for IP neuronal activities]. This statistical measurement is aimed at quantifying the degree of correlation between two signals and the corresponding time delays (Lopes da Silva et al., 1989; Slaght et al., 2004; Kalitzin et al., 2007) and has also been shown to give reliable measurements of the degree and direction of functional coupling between neuronal populations (Wendling et al., 2001; Pereda et al., 2005; Ansari-Asl et al., 2006).
The basic idea is that if the amplitude of signal YEMGOO(t) is considered a function of the amplitude of signal XNR(t), the value of the first signal—given a certain value of the second—can be predicted according to a nonlinear regression curve. A fixed duration sliding window is used, in which a parameter (nonlinear correlation coefficients: η<EMG|NR
2 and η<NR|EMG
2 in the opposite direction, ranging from 0 to 1) is computed. Low values of η2 denote that analyzed signals are independent. Contrarily, high values of η2 mean that one signal can be explained by a transformation (possibly nonlinear) of the other (i.e., the two signals are dependent).
By calculating the index η2, it is also possible to estimate the delay in the coupling between the electrophysiological signals during the classical eyeblink conditioning. The index η2 can be estimated as a function of time shift (τ) between signal XNR(t) and YEMGOO(t), or vice versa. That shift for which the maximum value of η2(τ) is reached is used as an estimate of the time lag between the two signals. This procedure, in contrast to “classical” linear cross-correlation functions, can yield an asymmetric function (i.e., η<EMG|NR
2 ≠ η<NR|EMG
2), which can provide insights into possible driving–response relationships (Lopes da Silva et al., 1989; Meeren et al., 2002). Details regarding the theoretical formulation of this method can be found in supplemental Appendix S1 (available at www.jneurosci.org as supplemental material). Readers may refer to the above-mentioned reports for a detailed and practical description of this technique.
Mean and SEM values for η and τ were calculated for each animal. To determine whether time delays differed significantly from 0, the 95% confidence interval for τ was calculated for each pair of electrophysiological recordings per trial. For all analyses, the amplitude of the signals XNR(t) and YEMGOO(t) was subdivided into 10 bins (Figs. 2 and 3A in the main text and supplemental Appendix S1, available at www.jneurosci.org as supplemental material), the maximum time shift was 370 ms, and the minimum time shift was 1.0 ms (one uniform time sample for all the signals). Maximum eyelid displacements during CRs were determined in the CS–US interval, and the function corresponding to the collected data (frequency sample at 1000 Hz) in the CS–US interval was adjusted by a simple regression method. This method allowed to fix the tendency for the points near the zero level of eyelid position and to establish a standardized algorithm for all the responses across of all the blocks of trials. In this way, the typical randomness in the determination of CR onset was avoided. The onset of a CR (Fig. 1K) was determined as the latency from CS presentation to the interception of the regression function with the maximum amplitude level (supplemental Fig. S3, available at www.jneurosci.org as supplemental material). This method was applied across the successive conditioning sessions, always showing the appropriate precision and robustness.
A scatter representation of motoneuronal (MN NR, in A) and IP neuronal (IP NR, in B) responses versus the EMG activity of the OO muscle (OO EMG). A, An LPA of the nonlinear regression curves was obtained by representing the scatter plot of OO EMG amplitude (in mV) versus the amplitude of motoneuronal responses (in mV) divided by segments of straight lines. The subdivisions in amplitude bins (Nb = 10, i.e., a total of 11 vertical and horizontal magenta lines) are shown for the association between MN recordings and the OO EMG activity. Illustrated results correspond to a single trial and to data collected from Figure 1J (for OO EMG) and Figure 1H (for MN NR) of the main text. B, The same as in A but relating IP neuronal (IP NR) responses versus OO EMG. Illustrated results correspond to data collected from Figure 1J (for OO EMG) and Figure 1F (for IP NR) of the main text.
A representation of the nonlinear association method applied to the EMG activity of the OO muscle (OO EMG) and facial motoneuron (MN NR) or IP neuron (IP NR) responses. Two recorded neurons were selected from each of the experimental subjects collected from the 10th conditioning session (series S1: C10, four cats, eight motoneurons; series S2: C10, four cats, eight IP neurons). A, Two pairs of LPAs of the nonlinear regression curves were obtained by representing the scatter plot of OO EMG amplitudes (in mV) versus neuronal recording amplitudes (in mV) by segments of straight lines. The subdivision in amplitude bins (Nb = 10; i.e., a total of 11 vertical and horizontal magenta and blue dashed lines) are shown for MN and IP neuronal associations with OO EMG activities (see color code in the legends). For each mean value per bin, the ± SEM is always indicated. B, The nonlinear association indices (η) were calculated between the two pairs of electrophysiological signals (in both senses) as functions of the time shift. The start of the eyelid CR is taken as the zero reference point of the dynamic association. Four nonlinear association curves are represented. Each representation corresponds to the average of all nonlinear association curves obtained for the trials presenting conditioned eyelid responses (n < 10) and collected from all the blocks (n = 12) of the same conditioning session (C10) for all the experimental subjects (n = 4). The vertical broken lines represent the time shifts for which the maximum values for nonlinear association indices are reached (time delays). In the first case (OO EMG vs MN NR, and vice versa), the information of asymmetry and the relative time delay in coupling were positive (Δη2max > 0, Δτ > 0, respectively), whereas the index of the directional coupling between the signals was equal to unity (D± = +1). In the second case (OO EMG vs IP NR and vice versa), the degree of asymmetry of the nonlinear coupling was positive (Δη2max > 0), and the relative time delay between the signals was negative (Δτ < 0), so that the direction index was annulled (D± = 0). In this case, the time delays are located to the right of the zero reference point (i.e., always lagged the start of CR).
(2) Causal inferences according to multivariate time series analysis.
We investigated the dynamic regression models and causal inferences between neuronal firing function [SIt:S1t(MN), with fMN(t) for facial MN or S2t(IP) with fIP(t) for IP neuron instantaneous firing frequencies] and learned motor response [S0t(θ), with θ(t) for eyelid positions during conditioned eyeblink responses] using the time-dependent causality analysis as a particular case of the transfer function method (TFM), a model frequently used to measure the functional interdependence between time series (Box and Jenkins, 1976; Granger, 1980; Nolte et al., 2008).
Relationships between Ns physiological time series [corresponding to instantaneous frequencies fMN(t) and fIP(t), and conditioned eyelid responses θ(t)] can be represented by transfer function models of the form
in which
are the impulsive responses or transfer functions of the models. B is the back-shift operator such that BSIt = SIt−1 and h = 1, 2, …, Ns. The moving average ωmhi(B) and autoregressive δahi(B) operators are also polynomials in B with orders m and a, respectively. The parameters bhi are nonnegative integers representing certain periods of delay in the transmission of the effects between the input SIit and output S0ht time series. The structure of the processes of inertia or uncertainties Uht = [ϕqh(B)/ϕph(B)]nht can be represented by the univariate operators (with orders p and q) in the stochastic difference equations of the followng form: ϕph(B)SIit = ϕqh(B)nht, in which nht are Ns independent Gaussian white-noise processes with variances σh and zero means (Tiao and Box, 1981). Details of the theoretical formulation of these models and their normal GI→0 and normalized RI→02 = 1 − e−GI→0 causality indices can be found in supplemental Appendix S2 (available at www.jneurosci.org as supplemental material).
Transfer function models of this form assume that the time series, when suitably arranged, possess a triangular relationship (Geweke, 1982; Harvey, 1994), implying for example that S2t depends only on its own past (i.e., the Granger causality indices are such that G0→2 ≈ 0 and G1→2 ≈ 0); S1t depends on its own past and on the present and past of S2t (i.e., G1→2 = 0 and G2→1 < 0, for unidirectional coupling); S0t depends on its own past and on the present and past of S1t and S2t (i.e., G0→1 = 0, G1→0 > 0 and G0→2 = 0, G2→0 > 0, respectively) and so on. If S1t depends on its own past and on the present and past of S2t, and S2t depends on its own past and on the present and past of S1t, then we must have a model that allows for this feedback (i.e., high and significant values of the causality indices in both senses G2→1 > 0 and G1→2 > 0, indicating bidirectional coupling). Readers may refer to Figure 5 in the main text and supplemental Appendix S2 (available at www.jneurosci.org as supplemental material) for a detailed and practical description of this method.
Results
We recorded a total of 105 posterior IP neurons, classified as type A (Fig. 1D,F,G). Type A neurons increase their firing in the time interval between CS and US presentations across successive conditioning sessions (Gruart et al., 2000; Sánchez-Campusano et al., 2007). In addition, we recorded 102 antidromically identified OO MNs (Fig. 1D,H,I). Characteristically, OO MNs encode eyelid position during conditioned eyelid responses (Trigo et al., 1999). The two pools of neurons were recorded in separate experiments during classical conditioning of eyelid motor responses (Fig. 1J,K) using a delay paradigm (Fig. 1C,E and see Materials and Methods). The present study was centered on the analysis of data collected at CS–US intervals across the successive sessions during the motor learning process, with emphasis in the 10th conditioning session, i.e., at the asymptotic level of acquisition of the associative learning test. A more detailed description of OO MN and IP neuron firing peculiarities during classical eyeblink conditioning can be found elsewhere (Trigo et al., 1999; Gruart et al., 2000; Sánchez-Campusano et al., 2007).
Directional coupling between simultaneously recorded signals during classical eyeblink conditioning
The degree of association between the EMG activity of the OO muscle [YEMGOO(t)] and crude recordings of NRs [XNR(t)] collected from OO MNs [XNRMN(t)] and IP neurons [XNRIP(t)] was obtained by computing the nonlinear association index (η) as a function of a time shift (τ) between these muscular and neuronal signals. The shift for which the maximum of η(τ) was reached provided an estimate of the time delay between the electrophysiological signals during the associative learning process. Thus, we were able to determine whether the maximum nonlinear correlation between neuronal recordings and EMG activities was before or after the zero reference point (i.e., the moment at which the conditioned eyelid response started; arrow in Fig. 1K).
Statistically, the nonlinear association index η quantifies the reduction of variance that can be obtained by predicting the YEMG00(t) values given the values of XNR(t) (or vice versa), on the basis of a regression curve. Details regarding the rationale underlying the nonlinear correlation analysis and its dynamic association functions are described in Materials and Methods and in supplemental Appendix S1 (available at www.jneurosci.org as supplemental material). The nonlinear association functions corresponding to the trials taken from all the blocks of the same session were averaged, session by session, for each animal across the 10 successive conditioning sessions (C01–C10).
Dynamic association functions enabled describing how the mean nonlinear association indices (η) evolved during the successive sessions. Two general nonlinear correlation analyses were developed. In the first case (OO EMG vs MN NR and vice versa), the nonlinear association functions remained quasi-stationary and with high (η ≥ 0.75) mean association index values across conditioning sessions—i.e., the two electrophysiological signals were very dependent. Thus, motoneuronal activities correlate significantly (one-way ANOVA F tests, H0 = 1, F(9,27,98) = 3.06, p < 0.01, for η<OO EMG|MN NR
; and H0 = 1, F(9,27,98) = 2.51, p < 0.01, for η<MN NR|OO EMG
) with the EMG activity of the OO muscle during the performance of conditioned eyelid responses in all the conditioning sessions. Importantly, their time delays in the two directions of coupling were always of opposed signs. In the second case (OO EMG vs IP NRs and vice versa), the values of the nonlinear association indices were statistically significant (one-way ANOVA F tests, H0 = 1, F(9,27,98) = 7.26, p < 0.01, for η<OO EMG|IP NR
; and H0 = 1, F(9,27,98) = 11.02, p < 0.01, for η<IP NR|OO EMG
), although their values show only a slight coupling (0.5 ≤ η < 0.8, i.e., the two signals were slightly related) between the neuronal activity recorded at the IP nucleus and the EMG activity of the OO muscle during conditioned eyelid responses. In this particular case, the maximum association index values always lagged the zero reference point across the successive conditioning sessions.
In Figure 3 is represented the nonlinear correlation analysis carried out with OO EMG corresponding to conditioned eyelid responses and crude recordings of neuronal activity collected from either OO MNs or IP neurons during the 10th conditioning session in different animals. This method has some major advantages over other signal–analysis procedures, such as coherence and cross-correlation functions, because it can be applied independently of whether the type of relationship between the two signals is linear or nonlinear (Lopes da Silva et al., 1989; Meeren et al., 2002). In all the analyses, the amplitudes of the signals XNR(t) and YEMGOO(t) were subdivided into 10 bins (Figs. 2 and 3A in the main text and supplemental Appendix S1, available at www.jneurosci.org as supplemental material). For each mean value per bin, the SEM is always indicated. Figure 2 displays data collected from a single trial (corresponding to recordings illustrated in Fig. 1F,H,J), whereas Figure 3A illustrates mean values of the linear piecewise approximations (LPAs) curves obtained for the trials (n = 9 trials per block; e.g., Fig. 1D) collected from all the blocks of the 10th conditioning session for all subjects (see Materials and Methods).
As illustrated in Figure 3B, the nonlinear association curves ranged from 250 ms before to 250 ms after the zero reference point. This representation enabled us to determine the time of occurrence of the maximum association index (ηmax) with respect to the beginning of the conditioned eyelid response (zero reference point, in Fig. 3B)—i.e., the time of delay in coupling between EMG and neuronal signals. Two of the illustrated traces (magenta traces) correspond to the nonlinear association curves between OO EMG and motoneuronal responses [YEMGOO(t) and XNRMN(t)], whereas the other two (blue traces) represent the nonlinear association curves between OO EMG and IP NRs [YEMGOO(t) and XNRIP(t)]. For these nonlinear association curves, the null hypothesis (H0, no correlation between the two distributions) were rejected for all the w-transformations (modified Fisher's z-tests, H0 =1; p[w<OO EMG|MN NR
] < 0.01, p[w<MN NR|OO EMG
] < 0.01, and p[w<IP NR|OO EMG
] < 0.05) of the nonlinear association indices (η) (see Materials and Methods and supplemental Appendix S1, available at www.jneurosci.org as supplemental material).
In the first case (Fig. 3B, magenta curves), the maximum association indices remained with high values in the dynamic correlation range (paired modified Fisher's z-test, H0 = 1, p[w<MN NR|OO EMG
, w<OO EMG|MN NR
] < 0.05, for the significance of a difference between two measured nonlinear correlation coefficients; mean ± SEM, ηmax<OO EMG|MN NR
= 0.957 ± 0.003 and ηmax<MN NR|OO EMG
= 0.925 ± 0.006). As a result, the information of asymmetry for the maximum values of the nonlinear correlation coefficients (or nonlinear association indices) (Δηmax2 = ηmax<OO EMG|MN NR
2 − ηmax<MN NR|OO EMG
2 ≈ 0.056) showed a deviation of only 5.6% (i.e., one signal can be explained as a quasilinear transformation of the other). The time delay in the direction of preferential coupling was positive [mean ± SEM, τ<OO EMG|MN NR
= 5.22 ± 0.09 ms, the XNRMN(t) signal preceding the YEMGOO(t) signal], whereas the time delay in the opposite direction was negative [mean ± SEM, τ<MN NR|OO EMG
= −4.81 ± 0.11 ms, signal YEMGOO(t) was delayed with respect to signal XNRMN(t)]. These data are clearly indicative of a unidirectional coupling [Δτ ≈ 10 ms and D± = +1, signal XNRMN(t) was in advance; see supplemental Appendix S1, available at www.jneurosci.org as supplemental material] between the two electrophysiological signals
during the 10th conditioning session, a result that was verified across the preceding conditioning sessions (data not shown). These results strongly suggest that OO MNs certainly encode EMG activity of the OO muscle at every instant of the dynamic nonlinear correlation range. This inference is supported by the experimental fact that the total number of spikes generated by these facial MNs and their corresponding discharge rate during the CS–US interval increased progressively across the learning process in relation to the increased number of EMG action potentials and therefore with the progressive increase in the amplitude of the corresponding conditioned eyelid responses (Gruart et al., 1995; Trigo et al., 1999; Sánchez-Campusano et al., 2007).
In the second case (Fig. 3B, blue curves), the nonlinear associations were only weakly significant in the dynamic correlation range (paired modified Fisher's z-test, H0 = 1, p[w<IP NR|OO EMG
,w<OO EMG|IP NR
] < 0.05, for the significance of a difference between two measured nonlinear correlation coefficients; mean ± SEM, ηmax<OO EMG|IP NR
= 0.744 ± 0.013 and ηmax<IP NR|OO EMG
= 0.543 ± 0.017). However, moderate and low nonlinear association index values are not necessarily indicative of an absence of functional relationships between the involved signals. Thus, the statistical algorithm developed here indicates that the degree of asymmetry in coupling for the maximum values of the nonlinear correlation coefficients (Δηmax2 = ηmax<OO EMG|IP NR
2 − ηmax<IP NR|OO EMG
2 ≈ 0.26, OO EMG activity variations can be explained by the IP neuronal activity better than vice versa) presented a deviation of 26% (i.e., one signal can be explained as a transformation, possibly nonlinear, of the other). In the same way, the time delays in the two possible directions of coupling were positive (mean ± SEM, τ<OO EMG|IP NR
= 15.62 ± 0.24 ms, and τ<IP NR|OO EMG
= 29.31 ± 0.33 ms), whereas the relative time delay was negative (Δτ ≈ −13.7 ms). According to these analyses, the functional interdependence between the simultaneously recorded signals was nonlinear and bidirectional [D± = 0, XNRIP(t)⇄YEMGOO(t); see supplemental Appendix S1, available at www.jneurosci.org as supplemental material] during the 10th training session, a result that was also verified for the preceding conditioning sessions (results not illustrated). In contrast with the first case, the bidirectional coupling indicated a feedback relationship between the OO EMG and IP neuronal signals [YEMGOO(t) and XNRIP(t)] recorded during the classical conditioning of eyelid responses. These data suggest that instead of one signal exclusively driving the other (as in the first case), both YEMGOO(t) and XNRIP(t) signals influence the other.
According to these results, the bidirectional coupling between OO EMG and IP neuronal recordings [YEMGOO(t) and XNRIP(t)] could be explained if we assume that there is an indirect reinforcement effect of IP neurons on OO MNs. This inference is supported by the experimental fact that IP neurons do not directly encode eyelid kinematics (Gruart et al., 1995; Sánchez-Campusano et al., 2007)—i.e., their contribution was only weakly significant in the dynamic correlation range (according to the linear and nonlinear correlation methods; Figs. 3 and 4). The sustained increase in mean peak firing rate presented by the pool of IP neurons across conditioning (Fig. 4A) was in parallel with a decrease in the mean time of its occurrence with respect to the start of the CR (Gruart et al., 2000; Sánchez-Campusano et al., 2007). The increase in firing rate, in association with the decrease in its mean time of occurrence, caused the maximum coefficient of correlation (linear/nonlinear) between IP neuron activity (instantaneous frequency/raw NR) and CR (eyelid position/OO EMG) to decrease [Fig. 4B and Table 1 (actual values characterizing the ongoing process of motor learning) in the main text]. As a result, values of maximum linear coefficient of correlation (rmax; one-way ANOVA F test, H0 = 1, F(9,27,98) = 161.54, p < 0.01) between instantaneous IP firing and eyelid position and the maximum nonlinear association index (ηmax; one-way ANOVA F tests, H0 = 1, F(9,27,98) = 7.26, p < 0.01, for ηmax<OO EMG|IP NR
; and H0 = 1, F(9,27,98) = 11.02, p < 0.01, for ηmax<IP NR|OO EMG
) between IP NR and OO EMG during the CS–US interval always lagged the zero reference point (i.e., the start of the CR; Figs. 1K and 3B). There is a significant decrease in the maximum correlation coefficients (R = −0.97, p < 0.01, for rmax〈fIP|θ〉; R = −0.76, p < 0.05, for ηmax<OO EMG|IP NR
; and R = −0.96, p < 0.01, for ηmax<IP NR|OO EMG
; all with respect to maximum instantaneous frequency fIPmax) during conditioning sessions (see e2, e3, and e4 regression lines, respectively, in Fig. 4B).
Relationships between type A IP neuron (IPn) firing rates and the percentage of CRs across the 10 conditioning sessions. A, Learning curve (blue circles; mean (± SEM) corresponding to all the animals (n = 4) of the experimental series S2). The mean (± SEM) firing rate of 96 identified IP neurons (blue squares) during the CS–US interval is also indicated. B, The abscissa and the left ordinate illustrate the relationship between the peak firing rate of IPn (fIPmax) and the percentage of CRs (blue squares; e1 regression line). The abscissa and the right ordinate illustrate the relationship between peak firing rate (fIPmax) and the coefficient of correlation (red triangles, rmax linear coefficients, e2 regression line; green triangles, ηmax nonlinear association coefficients) with eyelid performance during CRs, across conditioning. Note that the increase in fIPmax (e1 regression line) together with the decrease in its time of occurrence (always lagged the start of CR), caused a decrease in the rmax linear coefficient (see, e2 regression line) between instantaneous frequency fIP and eyelid CRs. In turn, a similar decrease was observed in the ηmax nonlinear coefficients (see, e3 regression line for ηmax[OO EMG vs IP NR]; and e4 regression line for ηmax[IP NR versus OO EMG]) between IP NR and the EMG activity of the OO muscle (OO EMG). Thus, the strength (strong or weak) of the functional nonlinear association depends inversely on the level of expression of conditioned eyeblink responses. Also, the degree of asymmetry of the nonlinear coupling was positive (Δη2max ≈ 0.26, with δηmax ≈ 0.20). The index R is the well-known Pearson's product moment correlation coefficient, i.e., the conventional index frequently used to measure the linear correlation between two variables (e.g., fIPmax vs % CRs; or fIPmax vs rmax; or fIPmax vs ηmax).
Mean values of peak firing rates (fIPmax, see Fig. 4A,B in the main text) of interpositus neurons (IP n, n = 96) and their latencies (tIPmax, with respect to CS presentation) for all the animals (n = 4) of the experimental series S2, the percentage of conditioned responses (% CRs in Fig. 4A,B), and the maximum values of the linear (rmax, in Fig. 4B) and nonlinear (ηmax in Fig. 4B) correlation coefficients across the successive sessions during motor learning
In turn, the discharge rate of the OO MNs increased progressively across the learning process, with a relative refractory period (minimum interspike time interval) that decreased progressively in the CS–US interval. Thus, an OO MN discharge pattern that correlated significantly with eyelid position during CRs in all the conditioning sessions was obtained. In this sense, posterior IP neurons could contribute to facilitate a quick repolarization process of OO MNs, reinforcing their tonic firing during the performance of conditioned eyelid responses (Sánchez-Campusano et al., 2007). This previous experimental evidence further supports the results obtained here with regard to the bidirectional interdependence between the OO EMG [YEMGOO(t)] and IP neuronal [XNRIP(t)] signals, and the unidirectional coupling manifested by the OO EMG [YEMGOO(t)] and OO MN [XNRMN(t)] signals.
Causal inferences and temporal ordering of the neuronal influences during motor learning
A question of great interest is whether there is a “causal” relationship between two simultaneous recordings collected during associative motor learning without any specific information on the direction of the coupling. Both the linear cross-correlation function and the nonlinear correlation method are, in principle, able to indicate the delay in coupling, but inferring causality from the time delay is not always straightforward (Granger, 1980; Lopes da Silva et al., 1989). Therefore, we decided to investigate the putative functional interdependences between neuronal activity [firing rates of OO MNs, fMN(t), or IP neurons, fIP(t)] and learned motor responses [i.e., conditioned eyelid responses, θ(t)], using time-dependent causality analysis.
The physiological time series fMN(t), fIP(t), and θ(t) exhibit nonstationary behaviors. The stationary time series S1t(MN), S2t(IP), S3t(SUM), and S0t(θ) were determined here by regular differentiation of averaged relative variation functions, and the latter as a result of high-pass filtering of integrated neuronal firing activities [resulting from fMN(t) for MNs and from fIP(t) for IP neurons] and of learned motor responses [resulting from eyelid position θ(t) during conditioned eyelid responses, CRs]. In Figure 5 are represented the transfer function models for the physiological time series corresponding to data collected from the 10th conditioning session. From here on, the significant values of the transfer function (or sample impulse response) are those that are outside the confidence bounds (horizontal blue dashed lines, ∼95% confidence interval), indicating the number of SDs of the sample impulse response estimation error to compute, assuming that the input and output physiological time series are uncorrelated.
Causal inferences using TFM between the kinetic neuronal commands [neuronal activity of OO MNs and neuronal activity of cerebellar IP neurons (IP)] and kinematics (eyelid CRs). Three representative IP neurons and a representative motoneuron were selected from each experimental subject (n = 8) during the 10th conditioning session (series S1: C10, four cats, four MNs; series S2: C10, four cats, 12 IP neurons). A–D, The transfer function models assume that the stationary time series S1t(MN), S2t(IP), S3t(SUM), and S0t(θ) have a functional and dynamic relationship, implying that S0t(θ) depends on its own past and on the past of S1t(MN), S2t(IP), and S3t(SUM) (i.e., the causality indices are such that G1→0 > 0, G0→1 = 0, νk+ ≠ 0, and νk− = 0, see A for a unidirectional coupling; G2→0 > 0, G0→2 > 0, νk+ ≠ 0, and νk− ≠ 0, see B for a bidirectional coupling; G3→0 > 0, G0→3 = 0, νk+ ≠ 0, and νk− = 0, see C for a unidirectional coupling, respectively). Note that S3t(SUM) depends only on its own past (i.e., G0→3 = 0) and that the unidirectional relationship in opposite direction (as shown in C) was possible because S1t(MN) and S2t(IP) were induced toward a phase equality (i.e., the relative phase difference will be close to zero). However, S1t(MN) depends on its own past and on the past of S2t(IP), and S2t(IP) depends on its own past and on the past of S1t(MN) (i.e., significant values of the causality indices in both senses: G2→1 > 0, G1→2 > 0, νk+ ≠ 0, and νk− ≠ 0, see D for a bidirectional coupling), indicating a feedback relationship between these physiological time series, at least in the statistical sense of causality. Blue horizontal dashed lines indicate the approximate upper and lower confidence bounds (∼95% confidence interval), assuming that the input and output physiological time series are completely uncorrelated.
As illustrated in Figure 5A, the functional interdependence between S0t(θ) and S1t(MN) was unidirectional (significant values of the transfer function alone for lags > 0, vk+ ≠ 0 and vk− ≠ 0)—i.e., the Granger causality indices are such that G1→0 > 0 and G0→1 = 0, which signifies that S0t(θ) depends on its own past and on the past of S1t(MN) and, as a result, OO MNs consistently lead the OO muscle during conditioned eyelid responses.
For the curves shown in Figure 5B, the significant values of the transfer function presented a bimodal distribution for both positive (vk+ ≠ 0) and negative (vk− ≠ 0) lags, indicating that S0t(θ) depends on its own past and on the past of S2t(IP), whereas S2t(IP) depends on its own past and on the past of S0t(θ). That is, significant values of the causality indices in both senses (G2→0 > 0 and G0→2 > 0), indicate a bidirectional coupling or feedback relationship between these time series.
In contrast, the relationship between S0t(θ) and S3t(SUM) was unidirectional (G3→0 > 0, G0→3 = 0, vk+ ≠ 0, and vk− ≠ 0), which indicates that this causal inference is dependent on the phase information, as shown in Figure 5C.
Finally, and as illustrated in Figure 5D, the functional coupling (or feedback relationship, the same as in Fig. 5B) between S1t(MN) and S2t(IP) was bidirectional in the sense of Granger causality (G2→1 > 0, G1→2 > 0, vk+ ≠ 0, and vk− ≠ 0), and these causal inferences signify that the neuronal activity in the posterior IP nucleus causes both the activity present in the final common pathway (i.e., that of the MNs participating in the generation of the selected motor responses) and the CRs of the eyelid motor system and vice versa [S1t(MN) depends on its own past and on the past of S2t(IP), and S2t(IP) depends on its own past and on the past of S1t(MN)].
At the same time, the transfer function models shown in Figure 6, D and E, assume that the stationary time series [S2τ(IP) and S0τ(θ) as shown in Fig. 6D; or S2τ(IP) and S1τ(MN), as indicated in Fig. 6E] possess an unidirectional interdependence after phase synchronization. The phase corresponding to S2τ(IP) in the instant τ = t − t1 − t2 was equivalent to the phases corresponding to S0τ(θ) and S1τ(MN) in the instants τ = t − t1 and τ = t, respectively, as illustrated in Figure 6A–C. The actual values for t1 (time elapsed from activation of MN firing to the zero reference point) (Fig. 1K) and t2 (time elapsed from the zero reference point to the activation of IP neuron firing) were 5.98 ± 0.26 (mean ± SEM; range, 3.41–8.56) ms and 23.5 ± 0.31 (mean ± SEM; range, 20.41–26.59) ms, respectively. With these conditions of phase synchronization, the Granger causality indices are such that G2→0 > 0 and G0→2 = 0, and vk+ ≠ 0 and vk− ≠ 0 (Fig. 6D), implying that S0τ(θ) depends on its own past and on the past of S2τ(IP); and G2→1 > 0, G1→2 = 0, vk+ ≠ 0, and vk− ≠ 0 (Fig. 6E), which signifies that S1τ(MN) depends on its own past and on the past of S2τ(IP). This phase analysis demonstrates that causal inferences are dependent on the phase information status, as indicated in Figures 5C and 6, D and E.
Causal inferences using TFM between IP kinetic neuronal commands (activity of IP neurons) and motor activities (activity of OO MNs and eyelid CRs), during the phase synchronization. A–C, High-pass filtering (HPF) of the integrated activity. Oscillatory curves (relative variation functions) resulting from HPF (−3 dB cutoff at 5 Hz and zero gain at 15 Hz) of integrated neuronal firing activities and of eyelid position corresponding to the same set of records. (The experimental data used here are the same as those used in Fig. 5). The operator ▿d allowed obtaining the stationary time series S1t(MN), S2t(IP), and S0t(θ) after making n = d regular differentiations to the nonstationary time series (i.e., the relative variation curves, as shown in A–C). Note that in the oscillating curves shown here, components a–d are totally out of phase with components e–h. The transfer function models (D and E) assume that the stationary time series possess a direct interdependence after phase synchronization [i.e., the phases corresponding to τ = t − t1 − t2, for S2τ(IP), τ = t − t1, for S0τ(θ), and τ = t, for S1τ(MN)], implying that S0τ(θ) depends on its own past and on the past of S2τ(IP) (i.e., the indices are such that G2→0 > 0, G0→2 = 0, νk+ ≠ 0, and νk− = 0; see D for a unidirectional coupling); and that S1τ(MN) depends on its own past and on the past of S2τ(IP) (i.e., G2→1 > 0, G1→2 = 0, νk+ ≠ 0, and νk− = 0, see E for another unidirectional coupling). Blue horizontal dashed lines in D and E indicate the approximate upper and lower confidence bounds (∼95% confidence interval), assuming the input and output physiological time series are completely uncorrelated.
Revealing the functional interdependence between the kinetic neuronal commands and the performance (kinematics) of learned motor responses
Traditional analytical methods, such as cross-correlation and coherence analyses, do not reveal causality, and their application demonstrates that these indices depend strongly on the type of relationship (linear or nonlinear) between the signals involved. The genesis of generalized causal inferences between cerebellar and brainstem motor circuits involved in the acquisition of conditioned eyeblinks and in the performance of learned eyelid responses requires at least a first approach in vivo to reveal the true causal relationships between the participating centers and the activity present in OO MNs.
We used multivariate analyses (combining the information of asymmetry, time delay, direction in coupling, and causality indices) of physiological time series recorded during classical eyeblink conditioning, using a delay paradigm, to reveal the putative functional interdependence and temporal ordering of the neuronal influences (i.e., neuronal activity in the IP nucleus and firing activities of OO MNs) during the actual motor learning process. Granger causality is a time-domain measurement of functional interactions, assuming directionality and information transfer. In short, the current study aimed to evaluate a novel approach for assessing directionality in cerebellar–MN network associations during the actual performance of learned movements.
First, a relatively weak (η < 0.8) but statistically significant functional coupling [one-way ANOVA F test, H0 = 1, F(3,9,2000) = 1844.35, p < 0.01; see blue curves in Fig. 3B of the main text, and Figs. S1 and S2 in supplemental Appendix S1, available at www.jneurosci.org as supplemental material] was found between the IP nucleus and the OO muscle across conditioning (one-way ANOVA F tests, H0 =1; F(9,27,98) = 3.06, p < 0.01, for η<OO EMG|MN NR
; and H0 = 1, F(9,27,28) = 2.51, p < 0.05, for η<MN NR|OO EMG
). The flow of information in the direction “cerebellar IP neurons→OO muscle” was more intensive (η<OO EMG|IP NR
≈ 0.74) than in the backward direction (η<IP NR|OO EMG
≈ 0.54). This is the first indication of asymmetry (Fig. 3B, blue curves with a deviation of 26%, Δηmax2 ≈ 0.26 and δηmax = ηmax<OO EMG|IP NR
− ηmax<IP NR|OO EMG
≈ 0.20 in Fig. 4B) in IP neuron–OO MN interactions during classical eyeblink conditioning. We concluded that the functional nonlinear association between the IP nucleus and the OO muscle is dynamic, bidirectional, asymmetric, and dependent (inversely) on the level of expression of conditioned eyeblink responses (Figs. 3B and 4A,B; e3 regression line for ηmax[OO EMG vs IP NR], e4 regression line for ηmax[IP NR vs OO EMG]).
Second, a strong and sustained increase in “OO MNs→OO muscle” interactions was found across the 10 successive conditioning sessions (one-way ANOVA F tests, H0 =1, F(9,27,98) = 7.26, p < 0.01, for η<OO EMG|IP NR
; and H0 = 1, F(9,27,98) = 11.02, p < 0.01, for η<IP NR|OO EMG
). During the 10th conditioning session—i.e., at the asymptotic level of acquisition of the associative learning test—“OO MNs→OO muscle” coupling remained constantly high [one-way ANOVA F test, H0 = 1, F(3,9,2000) = 1844.35, p < 0.01, see magenta curves in Fig. 3B of the main text, and Figs. S1 and S2 in the supplemental Appendix S1, available at www.jneurosci.org as supplemental material] and the association between the corresponding electrophysiological time series was linear (i.e., one time series can be explained as a quasilinear transformation of the other), unidirectional, and with scarce information of asymmetry (Δηmax2 ≈ 0.056, only 5.6% of deviation).
Third, functional information transfer in IP neuron–OO MN network associations during the performance of conditioned eyeblink responses required a driving common source. This common source is represented by the causal time series S?t or S?τ (see the causal relationships r5, R6, R7, and R8, marked with an asterisk in Fig. 7B) acting before or after phase synchronization, respectively. The best candidate for this role is the cerebral motor cortex that participates in both the generation and dynamic control of learned movements (Aou et al., 1992; Troncoso et al., 2007). Thus, motoneuronal responses and cerebellar IP activities contain equal amounts of information about the common source. Finally, we studied interdependences between two indirectly connected neuronal structures that communicate via the common source: “IP neurons→common source→final common pathway of learned motor responses.”
A, Diagrammatic representation of the causal inferences among the two neuronal time series [the OO MNs (see VII n inside the magenta rectangle) and the cerebellar IP neurons (see IP n inside the blue rectangle)] and the eyelid time series in the final common pathway for performance of learned motor responses (see CR, in the black rectangle). For each of the causal inferences represented, the implied physiological series and the respective temporal dependence (t or τ) are indicated. B, Causal relationships before (ri, in the instants t) and after (Ri, in the instants τ) phase synchronization. The relationship r1 represents a unidirectional linear model with a flow of S1t(MN)→S0t(θ); r2 a bidirectional model with flows S2t(IP)⇄S0t(θ); r3 a unidirectional linear model [S1t(MN)→S0t(θ)] with additional feedback [S2t(IP)⇄S0t(θ)]; r4 a unidirectional linear model [S1t(MN)→S0t(θ)] with additional feedback [S2t(IP)⇄S1t(MN)]; and r5 a combination between a unidirectional linear relationship [S1t(MN)→S0t(θ)] and a bidirectional model with a common source S?t [S?t→S1t(MN), S?t→S2t(IP), and S1t(MN)⇄S2t(IP)]. After phase synchronization (dashed-line rectangle in A and B), R1 represents a unidirectional linear model with a flow of S1τ(MN)→. S0τ(θ); R2 a unidirectional linear model with a flow of S2τ(IP)→S0τ(θ); R3 a closed semilinear model with flows S1τ(MN)→S0τ(θ) and S2τ(IP)→S0τ(θ); and R4 a unidirectional multilinear model with flows S2τ(IP)→S1τ(MN)→S0τ(θ). R5 represents a triangular relationship [nonlinear model with flows S2τ(IP)→S1τ(MN), S1τ(MN)→S0τ(θ), and S2τ(IP)→S0τ(θ), see dashed-line triangle in A] but with spurious causality [closed semilinear (R3) and unidirectional multilinear (R4) models cannot coexist]. However, in the r5, R6, R7, and R8 causal relationships, the mathematical linear coupling of r2 and R2 vanishes and the closed semilinear (r3 and R3) and nonlinear (R5) models were rejected. Finally, the putative causal inferences depended on a common source (see S?t or S?τ, in the relationships marked with an asterisk) in open nonlinear (r5, R6, or R7) or linear (R8) models such as those represented here. The causal relationships marked with an asterisk correspond to the best candidates to explain the coupling relationships in IP neuron-facial MN circuits involved in the dynamic control of conditioned eyelid responses.
With regard to the third conclusion, within the set D1 = [S0t(θ), S1t(MN)], the definition gives S1t(MN) causes S0t(θ) (Fig. 5A). Within the set D2 = [S0t(θ), S2t(IP)], it gives S2t(IP) causes S0t(θ) (Fig. 3B). However, within the set D12 = [S0t(θ), S1t(MN), S2t(IP)], neither S1t(MN) nor S2t(IP) causes S0t(θ), although the sum of S1t(MN) and S2t(IP) would do so (Fig. 3C). How can it be decided whether S1t(MN) or S2t(IP) are causal time series for S0t(θ)? The answer is, of course, that neither is. The causal time series is a common source, S?t (see r5 in Fig. 7B), involving other neuronal centers in the generation and performance of the learned motor responses. If the set of series within which causality was discussed is expanded to include S?t, then the above apparent paradox (Figs. 5A–C and 7A) vanishes. It will often be found that constructed examples which seem to produce results contrary to common sense can be resolved by widening the set of data for which causality is defined (see S?t and S?τ in Fig. 7B).
In Figure 7A are represented the functional relationships before and after phase synchronization. Before phase synchronization, the physiological time series possess a triangular relationship with feedback (the relationship between black, magenta, and blue rectangles in Fig. 7A, according to Fig. 5A,B,D). In contrast, after phase synchronization, the triangular relationship disappeared (in our case, a combination of unidirectional multilinear and closed semilinear models without feedback; see the dashed-line triangle in Fig. 7A and R5 relationship in Fig. 7B, according to Fig. 6D,E). In this situation, there was evidently a very scarce probability to high prediction ratio and clear autoregressive structure of a closed nonlinear model without feedback (see relationship R5 in Fig. 7B: a model with a low prediction ratio and unclear autocorrelation structure of prediction error). A practical solution was that the causal relationships r2 and R2 were not considered effective, and therefore, the closed semilinear (r3 and R3) and nonlinear (R5) models were rejected. Thus, the causal relationships r5 (open nonlinear model with feedback), R6 (open nonlinear model without feedback), R7 (open nonlinear model with feedback), and R8 (linear models with feedback) (see the causal relationships marked with an asterisk in Fig. 7B) are mathematically good candidates to explain the coupling relationships in IP neuron-facial MN circuits involved in the dynamic control of conditioned eyelid responses.
On one hand, the introduction of nonlinearity into the model may be necessary to get comprehensive information about network associations in cerebellar and brainstem motor circuits during the performance of the learned response. However, application of nonlinear autoregressive models with an additional “common source” requires more careful selection of model parameters (such as dimensions and nonlinear model functions). Altogether, the use of time-dependent causality analysis provided new information regarding neuronal network connectivity during learning processes. We consider this method a good alternative to the traditional measurements of functional interactions in cerebellar–brainstem motor circuits.
In all cases (causal relationships r5, R6, R7, and R8, marked with an asterisk in Fig. 7B), IP neurons were lagging OO MNs in a temporal window of τ − t = −(t1 + t2) ≈ −(29.48 ± 0.57) [mean ± SEM, range, −(23.68–35.29)] ms before their maximum neuronal activation. Thus, the IP nucleus detects the onset of the significant information from the MNs at t = τ (Figs. 6C and 7A) and uses that functional information to reinforce to a greater or lesser extent (adjusting their relative refractory periods in the CS–US interval) the performance of conditioned eyeblink responses at t = τ + t1 (Figs. 6B and 7A), modulating the phase of MN responses at t = τ + t1 + t2 (Figs. 6A and 7A). Finally, the cerebellum always lags eyelid activation, and cerebellar activities remain “silent” at the actual and current status during the start and preperformance of the learned motor response. Thus, it could be said that the cerebellum is looking back and reevaluating its own function with the information acquired in the process to play a modulating-reinforcing role in motor learning (Fig. 8).
Schematic representation of the reinforcing-modulating role of cerebellar posterior IP neurons during the acquisition of an associative learning task. Neuronal inputs arriving at the facial nucleus (OO MNs) and carrying conditioned signals p(t) need the reinforcing-modulating role of deep cerebellar nuclei signals m(t). To be efficient, IP neuronal signals need to go through a learning process to become 180° out-of-phase OO MN firing. Thus, IP neuronal activities (following a relay in the red nucleus) reach OO MNs right at the moment of maximum motoneuronal hyperpolarization (Trigo et al., 1999), and IP neurons facilitate a quick repolarization of OO MNs, reinforcing their tonic firing during the performance of classically conditioned eyelid responses. VII n, facial nucleus; IP n, cerebellar posterior IP nucleus.
Discussion
The present results allow proposing a definite role of IP neurons on the final common pathway for the eyelid motor system (i.e., the MN pool) that, by progressively inverting phase information, modulates and/or reinforces eyelid motor responses inversely (but not oppositely) to the initial contribution of OO MNs. Moreover, in accordance with previous reports (Gruart et al., 2000; Sánchez-Campusano et al., 2007; Poulet and Petersen, 2008), we have shown here that the neuronal activity in the IP nucleus does not lead the performance of learned motor responses but follows neural motor commands originated in different neuronal sources. As a whole, IP neurons could be considered to behave as a neuronal phase-modulating device supporting OO MN firing during learned eyelid movements.
To interpret our results, we have put together various theoretical considerations: theories regarding the distributed nature of the learning process (Llinás and Welsh, 1993; Bloedel and Bracha, 1995; Delgado-García and Gruart, 2006), the modulatory and reinforcing role of IP neurons on the final common pathway for the eyelid motor system (Gruart et al., 2000; Jiménez-Díaz et al., 2004; Sánchez-Campusano et al., 2007), and suggestions regarding the presence of “driving” and “modulating” connections in neuronal networks (Crick and Koch, 1998). We propose here that the “common source” is not merely a trigger but also acts as a distributor of significant functional information about the dynamic control of motor responses (see causal relationships r5, R6, R7, and R8, marked with an asterisk in Fig. 7B). In particular, the significant information in relation with the beginning and performance of conditioned eyelid responses may easily propagate to those neuronal centers (i.e., cerebellar IP neurons and facial MNs) that have multiple connections with the common source in cerebellar-brainstem circuits.
Although highly speculative at the present moment, we can suggest that the “driving common source” mentioned in the present study is located at the motor cortex and/or in related cortical areas. Indeed, many studies on cortical learning-related plasticity have shown changes in the size of cortical representations as a result of long-term changes in the strength of intracortical synaptic connections (Sanes and Donoghue, 2000; Hayashi et al., 2005). In this regard, it has been described recently a progressive increase of CS-evoked field potentials in the vibrissal motor cortex and a significant increase in the firing rate of layer V pyramidal neurons during the acquisition of a pavlovian associative learning task in behaving mice (Troncoso et al., 2007). Similar findings were reported in a seminal study on the firing of cortical neurons during classical eyeblink conditioning in alert cats (Aou et al., 1992). In the latter two cases, recorded cortical neurons fired in advance to the onset of the respective vibrissal or eyelid CRs. From a hodological point of view, selective motor-related cortical areas seem to project both to the cerebellum and to the facial nucleus as shown with different neuroanatomical and electrophysiological procedures (Fanardjian and Manvelyan, 1984, 1987; Morcuende et al., 2002; Ramnani et al., 2006; Glickstein and Doron, 2008).
To explain the quantitative differences in coupling strength between cerebellar IP neurons and facial nucleus MNs (Figs. 5D and 6E), we used two conceptually different kinds of network association—i.e., “driving” (Fig. 5A,C) and “modulating” (Fig. 5B,D before and Fig. 6D,E following phase synchronization) connections. We hypothesized that the not very efficient “IP neuron-MN” coupling via a common source (r5, R6, R7, and R8, marked with an asterisk in Fig. 7B) in cerebellar-brainstem motor circuits might prompt the sequential and out-of-phase release of the activity of the IP neurons (Sánchez-Campusano et al., 2007) and the propagation of the modulated information to the final common pathway of the eyelid system, which would certainly be helpful for a proper performance of consolidated conditioned eyeblink responses. These are “modulating connections,” which—although weak—are effective because they induce an optimized increase of firing activity in corresponding neurons (i.e., IP neurons and MNs) across the successive conditioning sessions (Trigo et al., 1999; Gruart et al., 2000; Sánchez-Campusano et al., 2007).
Thus, the IP nucleus could modulate the pathway of transition of functional information via bidirectional (see r5 in Fig. 7B) or unidirectional (see R6, and R7 in Fig. 7B) coupling with OO MNs or across a putative common source to the final common pathway of the eyelid motor system (see R8 in Fig. 7B). This optimized modulation may be a “modulating-reinforcing” indirect force for learned movements. Strong “common source activity→IP NRs” and “common source activity→OO MN responses” coupling may correspond to physically powerful “driving connections” that permit spreading of an appropriate rate of discharge of IP neurons and OO MNs during the CS–US interval and are, therefore, involved in the initiation and performance of learned motor responses.
In the present study, we identified type A neurons by their antidromic activation from the red nucleus, assuming that those neurons are the one involved in the control of OO MNs, because of their activation of red nucleus neurons projecting to the facial nucleus (Gruart et al., 2000; Morcuende et al., 2002; Jiménez-Díaz et al., 2004). But, according to a previous study (Gruart and Delgado-García, 1994), type A neurons can also be activated antidromically from the pontine nuclei, the medial longitudinal fascicle, and the restiform body. It is still possible that some of the type A neurons included in this study are GABAergic neurons projecting to the medial accessory subdivision of the inferior olive, because their axons pass very close to the red nucleus stimulating site. This latter possibility is suggestive of a putative involvement of inferior olive intrinsic circuits in the reinforcing role play by the cerebellum in this type of associative learning. For example, it has been recently shown that the electrotonic coupling among olivary neurons by gap junctions is fundamental for a proper timing of conditioned eyelid responses in wild-type and Cx36-deficient mice (Van Der Giessen et al., 2008 for details). Indeed, as shown in a previous study from our group, the reinforcing role of cerebellar circuits of ongoing conditioned eyelid responses is highly dependent on its adequate phase modulation with respect to intrinsic facial MN oscillatory properties (Sánchez-Campusano et al., 2007).
In summary, the reinforcing-modulating role of posterior IP neurons could help to explain cerebellar contribution to the proper performance of ongoing motor responses but without excluding prediction or a feedforward dynamics model. The IP neurons are able to exert a phase-modulating action with indirect predictive function. Thus, for a reevaluating process, the cerebellum must use their previous activity during the preceding quantum of movement to provide a phase-modulating action on the following quanta across the CR (Fig. 8). Something is thus in advance in the loop and the cerebellum is looking back for prediction, but their actual action occurs not before but during the learned action. Also, the contribution of the cerebellar cortex could play an important function on the actual phase-modulating mechanism of IP neurons during performance of learned movements. Purkinje cells, acting via their collateral and interneuronal networks in addition with the participation of the olivary system (Velarde et al., 2004; Bracha et al., 2009), may contribute for maintaining the information dynamics during and between trials (Hong and Optican, 2008). The recent demonstration of the existence of the “brain states” (Poulet and Petersen, 2008) determined by oscillation of the membrane potential in synchronized pyramidal cells may be compared with our experimental data showing similar oscillation in the OO MNs (Trigo et al., 1999; Sánchez-Campusano et al., 2007).
Footnotes
-
This work was supported by grants from the Spanish Ministry of Research and Innovation (BFU2008-00899/BFI) and Junta de Andalucía (BIO-122 and CVI-2487). We thank Drs. Fabrice Wendling and Guy Cheron for their comments and suggestions on a previous version of this manuscript. We thank Yulaima Favier and Roger Churchill for their editorial help.
- Correspondence should be addressed to Prof. José M. Delgado-García, División de Neurociencias, Universidad Pablo de Olavide, Ctra. de Utrera, Km. 1, Sevilla 41013, Spain. jmdelgar{at}upo.es
