Early monkey Corticomuscular coherence between motor cortex , somatosensory areas and forearm muscles in the monkey

Corticomuscular coherence has previously been reported between primary motor cortex (M1) and contralateral muscles. We examined whether such coherence could also be seen from somatosensory areas. Local field potentials (LFPs) were recorded from primary somatosensory cortex (S1; areas 3a and 2) and posterior parietal cortex (PPC; area 5) simultaneously with M1 LFP and forearm EMG activity in two monkeys during an index finger flexion task. Significant beta-band (^20 Hz) corticomuscular coherence was found in all areas investigated. Directed coherence (Granger causality) analysis was used to investigate the direction of effects. Surprisingly, the strongest beta-band directed coherence was in the direction from S1/PPC to muscle; it was much weaker in the ascending direction. Examination of the phase of directed coherence provided estimates of the time delay from cortex to muscle. Delays were longer from M1 (^62 ms for the first dorsal interosseous muscle) than from S1/PPC (^36 ms). We then looked at coherence and directed coherence between M1 and S1 for clues to this discrepancy. Directed coherence showed large beta-band effects from S1/PPC to M1, with smaller directed coherence in the reverse direction. The directed coherence phase suggested a delay of ^40 ms from M1 to S1. Corticomuscular coherence from S1/PPC could involve multiple pathways; the most important is probably common input from M1 to S1/PPC and muscles. If correct, this implies that somatosensory cortex receives oscillatory efference copy information from M1 about the motor command. This could allow sensory inflow to be interpreted in the light of its motor context.


INTRODUCTION
Oscillations in the beta band (15 30 Hz) are found throughout the motor system and behave in a task-dependent manner, but their function is as yet unknown. These oscillations appear as rhythmical fluctuations in both EEG and field potential recordings, as well as synchronous cell activity at the same frequency. They occur in sensorimotor cortex of both monkeys and humans (Tiihonen et al., 1989;Murthy and Fetz, 1992;Salmelin and Hari, 1994). Beta-band oscillations have been shown to be altered in various diseases including Parkinson's Disease and dystonia (Hutchison et al., 2004).
Coherence is the standard technique for measuring the strength of correlations between two signals in the frequency domain (Rosenberg et al., 1989). It provides amplitude and phase information about the relationship between two signals at a particular frequency. The phase lag between two signals, and its dependence on frequency, can also be used to infer a conduction delay if the relationship is purely feedforward (Rosenberg et al., 1989). More recently Granger causality has been used (Brovelli et al., 2004). This attempts to predict oscillations in one area from the past history of oscillations in another, and can help determine the direction in which oscillations are propagated.
Oscillations in primary motor cortex (M1) are coherent with EMG and cerebellar activity (Conway et al., 1995;Baker et al., 1997;Salenius et al., 1997;Halliday et al., 1998;Courtemanche and Lamarre, 2005;Soteropoulos and Baker, 2006). Early monkey ECOG recordings over S1 and EMG in humans suffering from intractable epilepsy. The slopes of the phase-frequency relationship for this coherence suggested that S1 led the EMG with a short delay. However this does not rule out an afferent contribution to corticomuscular coherence. Gerloff et al. (2006) studied patients suffering from congenital hemiparesis in whom M1 had relocated to the hemisphere ipsilateral to the paretic hand, whilst S1 remained on the contralateral side. In these subjects it was possible to separate signals from M1 and S1 in MEG recordings; corticomuscular coherence was seen from M1 but not S1. However, these patients have undergone substantial reorganization of their nervous system and it is unclear how much they reflect the relative contributions of M1 and S1 in normal humans. Tsujimoto et al. (2009) used field potential recordings in monkey to investigate corticomuscular coupling using the directed transfer function (DTF). They showed coherence between both pre-and post-central sites and EMG; the DTF suggested that the major causal influence was from cortex to muscle in all cases. However, this study was not able to separate the different somatosensory recordings by cortical area. Baker et al. (2006) reported that coherence between EMG and afferent discharge was considerably smaller for putative cutaneous fibers than those likely to be muscle spindle afferents. S1 areas 3b and 1 receive mainly cutaneous information, compared to the predominantly deep receptive fields of areas 3a and 2. An afferent contribution to corticomuscular coherence would therefore be most likely therefore in areas 3a and 2.
In this paper we use field potential recordings from identified motor, somatosensory and posterior parietal areas of the monkey cortex together with EMG to probe the contributions to corticomuscular coherence of the different sensorimotor areas. We apply both standard coherence analysis, and directed coherence (Granger causality). By targeting our recordings on cortical areas known to receive proprioceptive information, and by estimating the delays in the causal interactions between sites, we are able to address specific hypotheses on the pathways which might underlie the observed effects. Our results support the idea that beta-band oscillations are important in sensorimotor integration between motor and somatosensory/parietal areas.

BEHAVIORAL TASK
Two female rhesus macaques (M. mulatta) were trained to perform a finger flexion task for food reward. The index finger was inserted into a narrow tube, which splinted the finger and constrained movement to the metacarpophalangeal (MCP) joint. The tube was mounted on a lever which rotated on an axis aligned to the MCP joint. Lever movement was sensed by an optical encoder, and a motor exerted torque in a direction to oppose flexion. This was programmed to act like a spring (initial torque 48 mNm). The task required movement into target (6 , 12 , 18 or 24 flexion) and holding for 2 s (torque required at target either 64 mNm or 128 mNm). The four positions and two torques gave eight task conditions. The different conditions were presented in a pseudorandom order. Motor torque then rose, and the animal released the lever to obtain its reward (food rewards were given randomly every one to five trials). Visual feedback of target and lever position was provided via a computer screen. The analysis reported here focuses on the hold period of this task, as periods of steady holding have previously been shown to contain the strongest betaband activity (Gastaut et al., 1952;Baker et al., 1997Baker et al., , 2001.

SURGICAL PREPARATION
Following behavioral training, each monkey was implanted with EMG patch electrodes over 10 forearm and hand muscles, using the technique developed by Miller and Houk (1995). Wires from the EMG electrodes were tunneled subcutaneously to a connector on the back. Only recordings from first dorsal interosseous (1DI) and flexor digitorum profundus (FDP) are reported here, as these muscles displayed consistent activity during the task in both animals. We chose to include an intrinsic hand muscle as well as a forearm muscle as these should have different conduction delays from the cortex. The flexor digitorum superficialis (FDS) electrode gave a poor recording in monkey M. In monkey L the results for FDS (not shown) were similar to those for FDP. Cross-talk between these two EMGs was assessed to be negligible (magnitude of correlation between third-order differentiated unrectified EMG less than 0.02; Kilner et al., 2002). In a subsequent surgery the monkey received a stainless-steel headpiece (to allow head fixation) and a recording chamber placed over the central sulcus (Lemon, 1984;Baker et al., 2001). Surgeries were carried out under aseptic conditions, and general anesthesia (3.0-5.0% sevoflurane inhalation in 100% O 2 supplemented with a continuous intravenous infusion of alfentanil, 0.025 mg/kg/h). A full program of post-operative analgesia (10 µg/kg buprenorphine; Vetergesic; Reckitt and Colman Products, 5 mg/kg carprofen; Rimadyl; Pfizer) and antibiotic care (10 mg/kg cefalexin; Ceporex; Schering-Plough Animal Health or 15 mg/kg amoxycillin; Clamoxyl LA; Pfizer) followed surgery. All procedures were carried out under the authority of licenses issued by the UK Home Office under the Animals (Scientific Procedures) Act 1986, and were approved by Newcastle University's Local Ethical Review Panel.

RECORDING
In daily experiments, a 16 channel microdrive (Eckhorn and Thomas, 1993), loaded with glass-insulated platinum electrodes (M1) or tetrodes (area 3a, area 2, and area 5), was used to record LFPs, together with single unit activity as part of a related study. After completion of recording in M1, three microwires electrodes (50-µm diameter stainless steel wire insulated with Teflon, AM790500) were implanted transdurally and fixed in place to permit recording of M1 LFP simultaneously with other areas. LFPs (bandpass 1-100 Hz) were sampled continuously at 500 Hz, and saved to hard disk together with lever position, EMGs (bandpass 30 Hz-2 kHz, sampling rate 5000 Hz) and task behavioral markers. LFPs were recorded with negativity upward (inverting amplifier) and referenced to the stainless-steel headpiece. This large isopotential annulus probably yields a derivation similar to the linked-ear reference of human EEG. Similar results to those reported here were obtained using a bipolar derivation produced by differencing LFPs from the same area.
The different cortical areas were identified by electrode location, by a clinical examination of unit receptive fields and by noting the motor responses to intra-cortical microstimulation (ICMS). Cells and L is the total number of non-overlapping sections (Rosenberg et al., 1989). Several bins would be expected to be above the significance limit by chance, due to the multiple comparisons implicit in examining a coherence spectrum with many frequency bins. The number of points above the significance level defined by Eq. 1 which would be expected by chance at the 95% level was estimated using the binomial distribution. The coherence spectrum (between 0 and 45 Hz) was only considered significant if the total number of points above the significance line was greater than this value ( 4 points for LFP-LFP coherence and 5 points for LFP-EMG coherence). A similar approach to significance assessment was used in our previous work . Phase was calculated by taking the argument of the cross-spectrum. Confidence limits on the phase were calculated as (Rosenberg et al., 1989):

Directed coherence
Directed coherence was calculated using the method outlined in Baker et al. (2006), adjusted to use only data from the hold phase of the task. Rectified EMG and LFP recordings were down-sampled to a 200 Hz sampling rate. An autoregressive (AR) model of order 100 (see paragraph below for details on choice of model order) was fitted to the 400 points from each hold period of the task (using a publicly available program ArFit, Schneider and Neumaier, 2001). The AR model coefficients were then averaged together and used to calculate directed coherence and phase as in Baker et al. (2006). Although coherence can also be estimated from the AR model, we found this gave more variable estimates than the approach described above, which we used to match previous work in the field. In this paper, we used the normalization of directed coherence suggested by Geweke (1982).
where H ij is the directional transfer function representing the causal influence of signal j on signal i, H ii is the directional transfer function representing the causal influence of signal i on itself, C ii or C jj are the covariances of the noise innovations of each signal in the AR model, and complex conjugation is denoted by . This is similar to the approach we have used previously , but is more generally applicable to signals with non-white power spectra. Using this normalization, the directed coherence can be interpreted as the proportion of the variance in signal i which is explained by the past history of signal j (a coefficient of determination, Pierce, 1982). This contrasts to the alternative normalization used in the DTF, where the magnitude of a causal effect is expressed as a fraction of all causal effects on that signal (Tsujimoto et al., 2009). In the mathematical literature on AR models, much is made of the choice of model order. Several methods exist which seek to maximize the model's fit to the available data, whilst avoiding in M1 showed clear ICMS (13-18 biphasic pulses, 300 Hz, 0.2 ms per pulse, currents up to 60 µA) effects at low threshold ( 20 µA), usually digit or wrist movement. Cells in area 3a typically responded to single digit movement or wrist movement in one direction. ICMS effects were rarely observed in this area even at higher intensities (up to 60 µA, Widener and Cheney, 1997). Area 2 cells were located caudal to area 1 (cutaneous receptive fields) and tended to respond to single or multiple digit extension/flexion or wrist movement. The border between area 2 and area 5 was the most difficult to define. Most area 5 cells were located in the rostral bank of the intraparietal sulcus whereas area 2 cells were mainly on the post-central gyrus (Lewis et al., 1999). The area 5 cells had less well-defined proprioceptive receptive fields that were more often related to proximal joints (especially elbow) than cells in the other areas. Other more complex responses were seen in accordance with Mountcastle et al. (1975). Many area 5 cells were quiescent, unlike S1 cells, and only fired when the monkey began the task. Although we believe that most sites were correctly assigned, it is possible that the area 3a recordings included a few recordings from area 3b and the area 2 recordings included a small number of area 1 recordings.

ANALYSIS
Coherence and directed coherence analysis techniques were used to examine functional linkages between cortical LFPs and EMGs and also between M1 and S1/area 5 LFPs. Whilst coherence is a measure of synchronization between two areas at a given frequency, directed coherence assesses the ability of the past history of one signal to predict another. Significant directed coherence may indicate a causal influence of one signal on another.
In all cases, only data from the 2 s hold phase of the task were used. LFP recordings from all available electrodes (between 2 and 12) in each recording session were averaged to yield a single, lownoise estimate of oscillatory activity in that area (see . EMGs were full-wave rectified before further processing. Coherence and directed coherence were calculated by combining together all available recordings within a given monkey and area and treating them as one long recording. Before combination, LFP and rectified EMG from the different sessions were normalized by dividing them by their standard deviation. This procedure is similar to the "pooled coherence" method of Amjad et al. (1997) but avoids possible artifactually significant results if signal amplitudes covary across recording sessions (Baker, 2000). Data from all the different task conditions were used, as in a separate analysis no consistent differences were found between the coherence spectra for the different task conditions.

Coherence
To calculate LFP-EMG coherence, LFP was first up-sampled to 5 kHz to match the sampling rate of the EMG, and two nonoverlapping 4096 point long sections were used from each hold period for the Fourier analysis (for details of coherence calculation see Baker et al., 1997). For LFP-LFP coherence, the original waveforms sampled at 500 Hz were used, with a 256 point Fast Fourier Transform and four non-overlapping sections per hold period. This gave a frequency resolution of 1.22 Hz for LFP-EMG coherence and 1.95 Hz for LFP-LFP coherence. Coherence was considered significant (P < 0.05) if it was greater than Z where the binomial distribution was used to determine the number of significant frequency bins required to accept the spectrum as showing a significant effect (in this case, more than four bins in the 0-45 Hz range).

Phase delays
For both coherence and directed coherence phase delays were calculated by fitting a line to the phase -frequency plot using linear regression and calculating the slope; the slope divided by 2 gave the delay in seconds. Since phase-frequency relationships do not in general have zero intercept (Mima et al., 2000), both slope and intercept were unconstrained in the regression. Delays are presented as the maximum likelihood value returned by the regression fit, and the 95% confidence interval. All analysis routines were implemented in the MATLAB package (The MathWorks Inc, Natick, MA, USA).

Example of analysis
Figure 1 illustrates application of coherence and directed coherence to simulated data, with the aim of demonstrating the relative abilities of the two methods. In Figure 1A, two signals have been simulated as Gaussian white noise. Signal B receives a delayed input from signal A, with a delay of 20 ms. The coherence spectrum is flat in this situation (Figure 1Aa), correctly identifying the lack of any frequency specificity in the coupling. The coherence phase (Figure 1Ab) shows a linear phase-frequency dependence; the negative slope indicates that A leads B. The value of the slope implies a delay of 19.9 ms, close to the 20 ms delay used to simulate the data. The directed coherence for this situation has a magnitude similar to coherence for the A B direction, but is below the significance limit for the B A direction (Figure 1Ac). The phase of the A B directed coherence is similar to that for coherence, and linear regression reveals a delay of 19.9 ms from the slope of the phase-frequency relationship. For this simple situation, therefore, directed coherence adds little over-and-above that which can be learnt from coherence analysis. Figure 1B illustrates a more complex situation, where each of the two simulated signals has inputs to the other, but with different delays. This kind of reciprocally connected loop is often encountered in neural circuits. The coherence spectrum (Figure 1Ba) now shows oscillatory features, as the white noise inputs to the simulation interact with the feedback loop delay. Although the coherence phase (Figure 1Bb) still shows a linear relationship with frequency, the slope of this relationship does not relate in any simple way to the delays used to simulate the data.
By contrast, the directed coherence spectra (Figure 1Bc) are both flat, correctly identifying the lack of frequency specificity in the coupling between signals A and B. Both A B and B A directed coherence are above significance, corresponding to the reciprocal nature of the interaction. Finally, the slopes of the directed coherence phase-frequency relationships (Figure 1Bd) estimate delays of 19.9 ms and 29.6 ms, in close agreement to the known interaction delays. For reciprocally coupled systems, therefore, directed coherence appears to yield results which are easier to interpret than standard coherence analysis. excessive numbers of free parameters and consequent overfitting to a limited dataset. If such overfitting occurs, the model parameters begin to represent noise fluctuations in the data, rather than the genuine spatio-temporal relationships between the variables which are sought. This is important if the goal is to use the model to extrapolate from the recorded data to make predictions of future values (e.g., when attempting to predict financial time series). In this case, overfitting will lead to spurious predictions which extrapolate noise into the future. By contrast, in this work we use the model parameters to assess the strength of inter-relationships between the recorded signals. We therefore begin by choosing an arbitrarily high model order. We then carry out statistical tests on the resulting directed coherence spectra to assess which features are likely to result purely from noise, and which reflect genuine effects. This approach is comparable to the use of standard coherence spectra, in which we measure coherence at many frequencies, and then test which frequency bins have coherence significantly different from zero. The choice of model order is important mainly in so far as it alters the frequency resolution of the directed coherence. Using a model order of 100 for data sampled at 200 Hz will result in a frequency resolution of 2 Hz. Although apparently smoother spectra can be obtained using interpolation methods and lower model orders (Ding et al., 2000), the number of free parameters defining the spectrum remains equal to the model order. We chose an order of 100, because it provided sufficient frequency resolution to examine how the phase of the directed coherence varies with frequency. In cases where phase was linearly related to frequency, measurement of the slope of this relationship allowed estimation of the time delay, a key parameter in constraining hypotheses of what pathways might underlie the effects.
Importantly, the analysis reported here used very large datasets: between 2668 and 27460 trials of the task for each cortical area (mean 12319). AR models were thus fitted to between 1 and 11 million sample points, depending on the cortical area. In such circumstances, overfitting of the data is not a concern. A similar approach to model order was taken by Baker et al. (2006). Baker et al. (2006) found that the analytical significance limits normally applied to coherence were also appropriate for directed coherence. In the present work, we found that this was not the case, probably due to the use of data from discrete trials rather than a continuous recording. This was tested by simulating two uncorrelated white noise signals and fitting the AR models to small sections of the signals (400 points), corresponding to single "trials". We found that the analytical significance levels were too low and gave an excessive number of false positives ( 5%). Significance limits were therefore estimated by numerical Monte Carlo simulation. Two signals were generated as independent Gaussian random vectors with the same number of trials as the original data, and the directed coherence calculated. This was repeated 50 times with different random numbers. Directed coherence estimates at all frequencies and for all simulations were rank ordered, and the 95th percentile used as an approximate P 0.05 significance level when using that number of trials. Additionally a permutation test (not shown) was used to calculate significance; this estimated the significance limits by shuffling the original data, and produced similar results. As with the standard coherence spectra, underlying this coherence in more detail, directed coherence was also calculated ( Figure 2E). This showed a significant peak in the LFP EMG direction at 18 Hz, whilst for the EMG LFP direction it was very close to the significance limit, although there was a small peak at 40 Hz. In a linear system which incorporates a delay, phase should be linearly related to frequency. The slope of this relationship provides an estimate of the delay (Rosenberg et al., 1989). A regression line was fitted to the points of Figure 2D over the frequency band 14.6-28.1 Hz. The slopes implied a delay of 15.1 3 ms (regression slope 95% confidence limits) with cortex leading muscle. In these data from a single penetration, there were not enough significant directed coherence points (for either direction) to use in a regression analysis (Figure 2F).

LFP-EMG COHERENCE
Before combining together different recordings from one cortical area, it was important to check that they gave consistent results when analyzed individually. Figure 3A accordingly overlays coherence phase-frequency plots taken from single recording sites in monkey M. There was clearly a high level of consistency between sessions. Single site phase-frequency plots for monkey L showed similar internal consistency. The remainder of the analysis therefore combined all recordings from one monkey and cortical area, treating the available data as one continuous long recording. The very high signal-to-noise ratio which this produced allowed detailed analysis of both the magnitude and phase of the coherence and directed coherence.

HISTOLOGY
At the end of experiments, monkeys were deeply anesthetized (pentobarbitone, 60 mg/kg IP) and perfused through the heart with phosphate buffered saline (pH 7.2) followed by 4% formal saline fixative. For both monkeys 50 µm parasagital sections of the sensorimotor cortex were cut and stained with Cresyl Violet. These were used to confirm the location of the different cortical areas, and together with details of penetration sites and electrode depths to reconstruct the approximate locations of the recording sites. The electrodes rarely left visible tracks because of their small diameter (80 µm; see Mountcastle et al., 1991). Figure 2A shows estimated locations of the recording sites together with the estimated area boundaries for a parasagittal slice taken from monkey L (at 18 mm lateral to midline).

RESULTS
Recordings were made from a total of 28 sites in M1 (12 in monkey M and 16 in monkey L), 34 sites in area 3a (18 in M and 16 in L), 12 sites in area 2 (11 in M and 1 in L) and 19 sites in area 5 (6 in M and 13 in L). All recordings were within the hand representation, as assessed by the responses to ICMS (M1) or by single unit receptive fields (S1/area 5). Figure 2B shows example recordings of LFP, EMG (from the FDP muscle) and task lever position during a penetration into M1 in monkey L. The coherence between LFP and rectified EMG is shown in Figure 2C, with the corresponding phase-frequency relationship in Figure 2D. There was a clear 20 Hz peak in the coherence spectrum, which rose above the theoretical significance limit (dotted line, P 0.05). In order to investigate the pathways Figures 3B,C show the coherence spectra between the four cortical areas and two EMGs (FDP and 1DI) for monkeys M and L respectively. All plots have a clear peak at 20 Hz, indicating that not just M1, but also S1 and area 5 activity are coherent with EMG activity in the beta range. In monkey L both M1 and area 3a had small peaks at 40 Hz. The source of this coherence is unclear although corticomuscular coherence in the gamma range has previously been reported for a different type of movement (Omlor et al., 2007). Given its presence in only one of the monkeys, this peak was not investigated further.
In addition to the 20 Hz peaks, coherence was often significant at low frequencies (0-5 Hz). This coherence probably resulted from small residual movements present during the hold period, which would be represented both in the EMG and in the LFP (movement related potentials). Such movements tend to occur more often at the start of the task hold period than the end; in a separate analysis (not shown), we calculated coherence separately for the first and second halves of the hold phase. Significant coherence at low frequencies was largely confined to the first part of the hold phase, supporting the interpretation that this is the consequence of small movements immediately after entering the target zone.
Figures 3D,E present the phase-frequency relationships for the corticomuscular coherence, in different muscles and areas (only frequency points where there was significant coherence were plotted). The results from the two monkeys have been overlaid. In these plots a positive slope indicates that cortex leads EMG whilst a negative slope indicates that EMG leads cortex. The plots from areas 3a, 2 and 5 showed consistency with a positive slope in the beta band. By contrast, the coherence phase for M1 exhibited a more complex relationship; there was a negative slope for points outside the beta band, but a steep positive slope in the 13-25 Hz range which encompassed the peaks in the coherence spectra. These general traits were consistent across monkeys and for the two muscles studied. Regression lines were fitted to the phasefrequency plots to estimate the delays which these linear relationships imply (solid lines Figures 3D,E). The delays for the beta range (15-24 Hz; dotted lines in Figures 3D,E) are shown in Table 1. All slopes were significantly different from zero except for those for area 2 in monkey L.
Although coherence phase can provide useful information, its interpretation is complex in situations of reciprocal connectivity (Riddle and Baker, 2005;Baker et al., 2006). We therefore also calculated directed coherence to provide further insight into the mechanisms underlying corticomuscular coherence.  direction by contrast showed no consistent peak, although values were significant across a broad frequency range (Figures 4C,D). These results suggest that the beta oscillations in the cortex cause those in EMG; feedback influences from EMG to the cortex cover a much broader frequency range with no clear peak in the beta frequency band.

LFP-EMG DIRECTED COHERENCE
The results of the directed coherence analysis for the 1DI and FDP muscles are shown in Figure 4. In all areas investigated, both monkeys showed a peak close to 20 Hz for the LFP EMG direction (Figure 4A for monkey M and Figure 4B for monkey L; plots for 1DI and FDP are overlain). Spectra for the EMG LFP were 62.7 21.8 ms and 62.1 16.0 ms respectively. The most obvious route for somatosensory corticomuscular coherence is via connections to M1; this should give longer delays for somatosensory compared with motor corticomuscular coherence. Surprisingly, the phase delays for the somatosensory areas were actually shorter than M1 for both 1DI (area 3a: 33.0 14.2 ms; area 2: 41.1 16.5 ms; area 5: 38.1 18.9 ms) and FDP (area 3a: 43.4 8.5 ms; area 2: 21.6 17.5 ms). The phase-frequency plots showed no consistent trends for EMG LFP directed coherence (not shown).
The phase-frequency plots for LFP EMG directed coherence are shown in Figures 4E,F. There appears to be a linear relationship over the beta range for all areas, and results were consistent between the two monkeys. Linear regressions were fitted to points between 16 and 26 Hz, combining data from the two monkeys. For most areamuscle combinations the slopes were significantly different from zero (P 0.05). The only exception was the area 5 -FDP combination where there were not enough points for the regression analysis ( Figure 4F) the results were less consistent (Figures 6B,F). In monkey M there were clear peaks for all three areas ( Figure 6B); however these were at a lower frequency ( 15 Hz) than the peaks in the coherence spectra. In monkey L M1 had a clear causal relationship to area 2 at 20 Hz ( Figure 6F; red line), whilst the causal influences from M1 to areas 3a and 5 were smaller and at a higher frequency (20-30 Hz). The phase relationships were consistent across the three parietal areas (Figures 6,C,D,G,H), but varied slightly between the two monkeys, with longer delays for monkey L than monkey M in both directions. Lines were fitted to the data between 16 and 30 Hz for monkey M and 12 and 26 Hz for monkey L. For M1 to parietal cortex directed coherence in monkey L there was a phase shift between 16 and 18 Hz (Figure 6H), the cause of which is unknown. For this plot the line was fitted between 18 to 26 Hz. The phase delays are given in Table 2.

DISCUSSION
In this paper, we have shown that there is significant coherence between somatosensory/parietal areas and contralateral muscles. It has recently been demonstrated that beta-band oscillations in monkeys are actually larger in somatosensory/parietal areas than in M1 . Since there is coherence between muscles and the muscle spindle afferents , somatosensory/ parietal corticomuscular coherence might be expected as a simple consequence of oscillatory feedback from the periphery. However, surprisingly, directed coherence analysis indicated that oscillatory signals COHERENCE BETWEEN S1/AREA 5 AND M1 Figure 5 shows the coherence spectra between S1/area 5 and M1. As with the LFP-EMG coherence, all recording sessions within a given area and monkey were combined and treated as one long recording before coherence calculation. All areas in monkey L showed a 20 Hz peak in their spectra ( Figure 5B). In monkey M there were no clear peaks ( Figure 5A) although the coherence spectra around 20 Hz were above significance for all three areas. The lack of well defined peaks and lower coherence in general may be due to a lower signal:noise ratio in the microwire recordings from this animal. This has been discussed in a previous paper using data from the same monkeys . In both monkeys, there was no linear relationship between phase and frequency ( Figure 5C, regression slopes not significantly different from zero, P 0.05), and phases clustered around zero. This zero phase synchronization probably reflects the strong reciprocal connectivity between M1 and S1/area 5 areas; there may also be a contribution from electrical cross-talk between the recordings from these closely-situated regions. This is a problem for coherence analysis, however directed coherence ignores correlations at zero lag including those generated by electrical cross-talk.

DIRECTED COHERENCE BETWEEN S1/AREA 5 AND M1
In agreement with a previous study (Brovelli et al., 2004), we found that the S1 and area 5 have a clear Granger-causal relationship in the beta band to activity in M1 (Figures 6A,E). In the reverse direction is often non-zero, implying the existence of a constant phase offset as well as a constant delay (Mima et al., 2000). A phase intercept of -/2 might be expected based on the known relationship of LFP to spiking activity . Riddle and Baker (2005) observed sloped phase-frequency relationships with implied delays of 10 ms in some human subjects (compared with corticomuscular conduction delays estimated using transcranial magnetic brain stimulation, TMS, of 22 ms), whilst others had flat phase-frequency plots. In this work, data from one monkey yielded M1 corticomuscular delays from coherence analysis of 8-11 ms, which is consistent with the known conduction times between M1 and hand/forearm muscles (Palmer and Fetz, 1985;Lemon et al., 1987). For the second monkey the delays were longer than expected (20-30 ms). Riddle and Baker (2005) attributed the different phase delays across their pool of human subjects as reflecting different, subjectspecific contributions from feedforward and feedback pathways. Similar factors may explain the difference between coherence phase delays measured from the two monkeys used here. This could be a genuine difference in the neural circuitry of the two animals, or alternatively could reflect differences in the way the monkeys performed the behavioral task. For example, our task required the monkeys to place the index finger in a tube. Detailed differences in hand posture could markedly alter the afferent feedback during the task, and change the balance between feedback and feedforward pathways.
In contrast to the delays estimated from coherence, those found from directed coherence analysis in M1 were substantially longer ( 60 ms) than direct estimates of conduction delays in the macaque in somatosensory/parietal regions flowed mainly from cortex to the periphery rather than vise versa -the same direction as that observed for M1. These results agree with the findings of Ohara et al. (2000), who found that there was coherence between S1 (ECoG recordings) and EMG, and that the slope of the phase-frequency relationship indicated that S1 led EMG. They also agree with the work of Tsujimoto et al. (2009), who showed that causal influences from both S1 and M1 to muscle were primarily in the descending direction.

ESTIMATED CONDUCTION DELAYS
A powerful feature of coherence analysis is the ability to extract the time delay associated with communication between two areas, by examining the slope of the phase-frequency relationship. However, there is controversy in the literature over the delays indicated by corticomuscular coherence from M1. Some authors report delays broadly consistent with conduction over fast corticospinal pathways (Brown et al., 1998;Gross et al., 2000;Mima et al., 2000), although the actual values are often slightly smaller than expected. Others report flat phase-frequency relationships . In addition, even if linear relationships are observed, the phase intercept  (A), but for M1 3a/2/5 direction. The highest significance level out of area 3a, area 2 and area 5 is shown but the difference in 95% significance levels between the different areas is small compared to the size of the directed coherence. (C,D) as (A,B) for monkey L. (E,G) phase spectra for areas 3a/2/5 M1 directed coherence for monkeys M and L respectively. (F,H) phase spectra for M1 areas 3a/2/5 directed coherence for monkeys M and L respectively. For (E-H) regression analysis was used to fit phases at frequencies between dotted lines; black line represents best fit. The measured delays up to 49 ms may therefore not be an unreasonable estimate of the time for oscillatory propagation from one area to another. Finally, it should be noted that both coherence and directed coherence are linear analysis methods. Neural systems are known to contain non-linear elements (e.g., the spiking threshold), and the numerical magnitude of coherence can be substantially affected by these non-linearities Soteropoulos and Baker, 2006). It is possible that non-linearities could also affect estimates of delay. However, in most cases we found significant linear relationships between phase and frequency. This indicates that a major part of the transformation from one signal to another can be modeled as a simple delay; any residuals from a straight-line relationship may result from non-linearities.

POSSIBLE PATHWAYS
Several pathways could provide the substrate for descending corticomuscular coherence from S1. Figure 7 shows four schematics of possible pathways using the main known connections between M1, S1 and the periphery. One possibility is that M1 and S1 activity is synchronized by common input from another region (Figure 7, Pathway 1); this is a reasonable option, as -for example -both SMA (Darian-Smith et al., 1993) and SII (Jones, 1986) are known to provide input to primary motor and somatosensory cortices. However, given the close spatial proximity of the pre-and post-central regions which we investigated, we would expect that conduction delays from the distant area to each would be similar. If synchronization with muscle was then mediated via fast corticospinal axons from M1, this would lead to similar delay estimates from all of the studied cortical areas to muscle. In fact, the delays differed between parietal and motor cortex.
Alternatively, oscillations from S1 could propagate to M1, and thence to the periphery (Figure 7, Pathway 2). There are strong bidirectional anatomical links between the somatosensory areas and M1 (Jones et al., 1978), and directed coherence shows strong influences from S1 to M1 (Brovelli et al., 2004;present Figures 6A,C). However, if this pathway is correct, the delay from S1 to muscle should be longer than from M1. In fact, corticomuscular delays estimated from directed coherence were shorter for parietal regions than for M1.
The corticospinal tract does not originate solely in motor cortex; around 40% of fibers come from S1/area 5 (Philips and Porter, 1977). These axons mainly terminate in the dorsal horn, where they are assumed to influence spinal sensory relays (Galea and Darian-Smith, 1994;Seki et al., 2003). An exception is some connections from area 3a, which have recently been shown to make direct monosynaptic connections to motoneurons (Rathelot and Strick, 2006). Because microstimulation of area 3a does not usually produce overt movements (Widener and Cheney, 1997), Rathelot and Strick (2006) suggested that area 3a corticomotoneuronal cells may target gamma motoneurons. In that case, it is conceivable that oscillations could pass from area 3a to gamma motoneurons (connections E F, Figure 7, Pathway 3). Muscle spindle nuclear chain fibers exhibit brief twitch times (Bessou et al., 1968), such that beta-band activity in static gamma motoneurone discharge should be passed reasonably faithfully to the discharge of spindle afferents. Group Ia spindle afferents have monosynaptic connec-from stimulation. They were also highly consistent between the two animals ( Figures 4E,F). Using computational modeling Williams and Baker (2009) showed that several factors can lead to increased delay estimates from coherence measures, including the time course of unitary EPSPs in motoneurons, and the width of the motor unit action potentials recorded in EMG. It is possible that these considerations are partly responsible for the long delays found from M1 using directed coherence analysis in the present work. However, not all factors influencing coherence lead to longer delays. For example, motoneuron firing properties can introduce a phase advance, which shortens the estimated delay (Williams and Baker, 2009). In all, no satisfactory explanation for the long measured delays can be found at the present time.
Coherence analysis between M1 and S1/area 5 revealed zerophase synchronization of oscillatory activity. This is a common finding in studies of distant brain regions, and has been reported previously between both visual and motor areas across hemispheres (Engel et al., 1990;Murthy and Fetz, 1996;Roelfsema et al., 1997), as well as between M1 and cerebellum (Soteropoulos and Baker, 2006). Studies using computational models have shown that reciprocally connected areas can synchronize at zero lag so long as part of the coupling is inhibitory (Van Vreeswijk et al., 1994;Ernst et al., 1995;Gerstner et al., 1996;Sturm and Konig, 2001). Rajagovindan and Ding (2008) examined empirically the factors which could produce synchrony at zero-lag, concluding that it can result from positively correlated common input, reciprocal interaction or a combination of these two mechanisms. However, it is possible that electrical cross-talk between M1 and S1/area 5 recordings could have produced zero-phase synchronization artifactually. This is unlikely: available evidence suggests that LFP signals spread no more than a few millimeters (Steriade et al., 1996;Kreiman et al., 2006). It may represent even more localized activity: a recent report suggests that LFP in primary visual cortex represents the activity of cells within just 250 µm of the recording site (Xing et al., 2009). Because it does not depend on analysis of simultaneous sample points (the past history of one signal is used to predict the present value of another), electrical cross-talk cannot affect directed coherence.
Directed coherence in the M1 S1/area 5 direction yielded delay estimates which varied between the two animals (31.8 ms and 48.9 ms), as did the estimates for the opposite direction (10.1 ms compared with 36.4 ms). A delay of up to 49 ms for conduction between M1 and the very close somatosensory areas investigated appears excessively long. Following stimulation of somatosensory cortex in the awake rabbit, units in M1 respond with latencies as short as 1 ms (Swadlow, 1994). However, although the fastest cortico-cortical axons may conduct at velocities of 10 m/s, they are greatly outnumbered by fibers with lower conduction velocity (Swadlow, 1990(Swadlow, , 1991(Swadlow, , 1994; the mean conduction velocity for axons projecting from M1 to S1 has been estimated as 0.9 m/s. The actual mean conduction velocity is likely to be lower due to preferential sampling of large cells both for stimulation and recording (see Swadlow, 1998). Although the fastest fibers are responsible for the onset latency of stimulation effects, measures of oscillatory coupling are more likely to be sensitive to the mean conduction time across the whole population of fibers. Additionally, just as for corticomuscular coherence, the time course of the post-synaptic EPSP is also likely to add extra delays (Williams and Baker, 2009 Whilst it is tempting to interpret significant directed coherence as indicating a causal influence, other network topologies may produce this result. We have demonstrated that there is directed coherence in the S1/area 5 M1 direction (connection D, Pathway 2), however, it is also present in the other direction (connection G, Pathway 4), although rather weaker (Figures 6B,D). Somatosensory/ parietal areas and the motoneurons therefore receive common input from M1, which could lead to the observed corticomuscular coherence. Because the delay from M1 to S1/area 5 (connection H, mean delay 40 ms across two animals; Figures 6F,H) is shorter than from M1 to EMG (connections B C, delay 61 ms for 1DI; Figure 4E), directed coherence would detect an influence in the S1/area 5 EMG direction. The observed delay from S1/ area 5 to EMG based on the directed coherence phase is broadly consistent with that expected if this interpretation is correct (connection G 21 ms versus observed delay of 35 13 ms for 1DI). It is possible that the measured delay between S1/area 5 and EMG reflects a mixture of two pathways, the first a causal influence via the motor cortex (Pathway 2, Figure 7), and the second common input from motor cortex to S1/area 5 areas and EMG (Pathway 4, Figure 7). This might explain why the measured delay is larger than that expected for the common input pathway alone. Support for this explanation of parieto-muscular coherence comes from the recent study of corticomuscular coherence in hemiparetic patients by Gerloff et al. (2006). These patients have an abnormal separation across the midline of M1 and S1 relating to one hand. If this resulted in reduced strength of interaction between somatosensory and motor areas in these patients, M1 would provide reduced common input to S1 and muscle, and corticomuscular coherence from S1 would be reduced.
In this paper, we have chosen to use bivariate methods (coherence and directed coherence). Further extensions to both coherence and directed coherence are possible, in which multivariate analysis is used to reveal the underlying network Blinowska et al., 2004;Barrett et al., 2010). However, we are concerned that such methods may prove especially sensitive to the non-linearities of neural data. We already know that the magnitude of coherence can be markedly altered by these non-linearities Soteropoulos and Baker, 2006). However, if coherence remains above significance, our conclusion that there is an interaction between two signals is unaffected. Consider by contrast one important multivariate method, namely partial coherence. This seeks effectively to predict coherence between two signals B-C, based on the measured coherence A-B and A-C. The measured coherence B-C is then compared with the prediction, to discover whether there is any additional B-C interaction other than that expected due to the interaction of B and C with A. In this case, the impact of non-linearities on coherence magnitude estimation will severely distort the calculation, and may lead to erroneous inferences. It is likely that other multivariate methods will suffer from similar problems. Until the impact of this has been fully assessed, we feel that our bivariate analyses are at the limit of what may currently be considered reliable.

FUNCTIONAL SIGNIFICANCE
As discussed above, interpretation of directed coherence as indicating causal connections is difficult. However, significant directed coherence does unambiguously show that activity in one area tions to alpha motoneurons, and are capable of passing beta-band activity to EMG . A route involving area 3a, gamma motoneurons and muscle spindle afferents could therefore conceivably generate coherence between parietal areas and muscle (Figure 7, Pathway 3).
Several pieces of evidence argue that corticomotoneuronal connections from area 3a are not the primary pathway underlying our observations. Firstly, using both coherence and directed coherence measures, results for area 3a were very similar to those from areas 2 and 5 (Figures 4 and 5). However, only area 3a has corticomotoneuronal axons. Secondly, the passage around the peripheral stretch reflex loop would introduce substantial additional time delays, compared with the pathway from M1 involving a direct corticomotoneuronal connection to alpha motoneurones. Yet the estimated delays from parietal areas to EMG were actually shorter than from M1. FIGURE 7 | Schematic showing possible pathways underlying coherence between S1/area 5 and EMG. Lettered labels are referred to in text.

www.frontiersin.org
July 2010 | Volume 4 | Article 38 | 13 predicts that in another. Directed corticomuscular coherence from S1 was often as large or larger than that from M1, meaning that somatosensory activity could accurately predict the peripheral motor output. Common input to S1 and muscle from M1 could thus be viewed as the supply to S1 of efference copy information. This would place S1 in a good position to predict not just the evolving motor output, but also its sensory consequences. Interpretation of sensory input critically depends on placing that input into the context of the movements that elicited it (Blakemore et al., 1998;Bays et al., 2005Bays et al., , 2006. The present observations may therefore further support previous suggestions that beta-band oscillations have a role in sensorimotor integration (MacKay, 1997;Brovelli et al., 2004;Riddle and Baker, 2005;Baker et al., 2006).