Time and frequency domain methods for quantifying common modulation of motor unit firing patterns
© Myers et al; licensee BioMed Central Ltd. 2004
Received: 30 August 2004
Accepted: 14 October 2004
Published: 14 October 2004
In investigations of the human motor system, two approaches are generally employed toward the identification of common modulating drives from motor unit recordings. One is a frequency domain method and uses the coherence function to determine the degree of linear correlation between each frequency component of the signals. The other is a time domain method that has been developed to determine the strength of low frequency common modulations between motor unit spike trains, often referred to in the literature as 'common drive'.
The relationships between these methods are systematically explored using both mathematical and experimental procedures. A mathematical derivation is presented that shows the theoretical relationship between both time and frequency domain techniques. Multiple recordings from concurrent activities of pairs of motor units are studied and linear regressions are performed between time and frequency domain estimates (for different time domain window sizes) to assess their equivalence.
Analytically, it may be demonstrated that under the theoretical condition of a narrowband point frequency, the two relations are equivalent. However practical situations deviate from this ideal condition. The correlation between the two techniques varies with time domain moving average window length and for window lengths of 200 ms, 400 ms and 800 ms, the r2 regression statistics (p < 0.05) are 0.56, 0.81 and 0.80 respectively.
Although theoretically equivalent and experimentally well correlated there are a number of minor discrepancies between the two techniques that are explored. The time domain technique is preferred for short data segments and is better able to quantify the strength of a broad band drive into a single index. The frequency domain measures are more encompassing, providing a complete description of all oscillatory inputs and are better suited to quantifying narrow ranges of descending input into a single index. In general the physiological question at hand should dictate which technique is best suited.
Common oscillations in neurophysiological activity in the human motor system have been well documented. During voluntary muscle contraction, the human central nervous system drives motor neurons at a range of frequencies which cause common modulations in the firings of these neurons. These drives are reviewed in  and  where they are summarized into four broad frequency ranges: (1) A low frequency drive at around 1–3 Hz (2) A neurogenic component of physiological tremor that occurs between 5–12 Hz and is likely to have both spinal and supraspinal components. (3) A corticospinal drive in the beta (15–30 Hz) range (4) A corticospinal drive in the low gamma (30–60 Hz) range, that increases in importance with stronger contractions and is called the Piper rhythm.
There are two distinct approaches toward the identification of these drives. The majority of the literature has examined common modulation to motor units using frequency domain methods. This methodology was first introduced by Rosenberg and colleagues  and applied by Farmer and colleagues  who used coherence analysis to identify both a significant low frequency and beta-band association between motor unit firings in the 1–12 Hz and 15–30 Hz frequency ranges respectively. Subsequently, coherence analysis has become an established technique to study bivariate motor system measurements and a number of works have used this to investigate corticomuscular interactions [5–8, 1]; tremor ; aging ; oscillatory drive in Parkinson's disease [11, 12]; dystonia ; stroke ; and cortical myoclonus .
A separate body of literature has focused specifically on the low frequency common drive. This drive was first identified by De Luca and colleagues  who demonstrated that the firing rates of concurrently active motor units (MUs) were modulated in a highly interdependent manner. They low-pass filtered the impulse trains corresponding to MU firing times to obtain the time-varying mean firing rates which they high-pass filtered at 0.75 Hz. They then performed a time domain cross-correlation analysis between pairs of zero-mean signals representing the fluctuations in mean firing rates. Peaks occurring near the zero lag location in the normalized cross correlations implied that those firings rates were essentially simultaneously modulated with virtually no time delay. This phenomenon was termed 'common drive' to indicate a common excitation to the motor neuron pool that results in concurrent fluctuations in the firing rates of motor units from the same pool. A number of subsequent studies have utilized this technique to investigate the relationship of this drive to handedness [17–19]; different proprioceptive conditions ; exercise ; task and disease ; and aging . These works have established this time domain method as an accepted means of quantifying the common low frequency modulation of MU firings.
In a recent review , it was suggested that the low frequency common drive first identified by De Luca and colleagues  using time domain methods is essentially the same low frequency drive as detected by Farmer and colleagues  using frequency domain methods. In this paper we explore the relationship between the two techniques using mathematical and experimental approaches.
Frequency domain methods: Coherence
The coherence between two zero-mean stationary random processes x1 (t) and x2 (t), at frequency f, is defined as:
where * denotes complex conjugation, N is the number of data segments employed and X1n(f) and X2n(f) are the discrete Fourier transforms of the n th data segments of x1 (t) and x2 (t). This estimate is biased and its probability density function for non-weighted and non-overlapping windows has been analytically determined . This may be used to derive the value of the estimated coherence, with a particular probability of occurrence, α, that would be obtained when the true value is zero. Any value exceeding this level is considered to be unlikely to be a false indication of coherence with (α × 100) % confidence. This confidence level is given by [24, 3]
E α = 1 - (1 - α)1/(N-1) (3)
Time domain methods: Cross correlation
The cross correlation between two zero-mean stationary random processes x1 (t) and x2 (t) is defined as:
where E [·] is the estimation operator. Assuming ergodicity, for single time-limited realizations of each random process, this is determined using the integral:
where * denotes complex conjugation and τ is the time lag between the signals. The Fourier transform of the cross correlation function, defines the cross-spectrum, (f). Cross correlation functions are unbounded measures and are typically normalized by the values of the autocorrelations at zero lag to bound the estimate between -1 and 1. The autocorrelation functions are the time domain equivalent of the auto power spectra and their value at zero lag represents the total energy in the signal. The normalized and bounded measure is known as the cross correlation coefficient, , which provides a measure of the linear association between the two signals at a given time lag and is given by:
Figure 2c depicts the normalized cross correlation function for the same two motor unit spike trains used in Figure 1 obtained using a moving average Hanning window of 400 ms duration and high pass filtering at 0.75 Hz. The function is displayed for lags up to ± 400 ms and the peak of the signal indicated at a lag of 3.5 ms. The time domain common drive coefficient is measured as 0.75.
Relationship between cross correlation and coherence
The cross correlation coefficient is related to coherence using a similar analysis to Gardner  as follows:
We begin with the real and stationary signals x1 (t) and x2 (t), where x2 (t) = α x1 (t - τ0) + n (t) is a scaled and time-delayed version of x1 (t), with additive uncorrelated, zero mean noise, n (t). The cross-correlation function is given by
since the cross correlation with the noise is zero everywhere. The cross spectrum is given by
Assume that the signals y1 (t) and y2 (t) result from passing the signals x1 (t) and x2 (t) respectively through a tunable narrowband bandpass filter with transfer function denoted by H (f):
Similarly for the autocorrelation functions we get
Thus as Δ → 0 we obtain the expression for the normalized cross correlation function as:
The peak of the cross-correlation function occurs at the time delay, τ = τ0. Thus
Thus we see that the peak of the normalized cross-correlation function between two signals after ideally bandpass filtering to contain a single frequency, is identical to the magnitude of the coherence function of the original signals at the frequency of the filter. The phase of the coherence function is the same as the phase of the cross-spectrum and provides the time delay.
For a less ideal filter that spans several frequencies the relationship is less precise and may be derived as follows:
Let W (f) be the new filter transfer function and thus the normalized cross-correlation function is:
where f1 and f2 are the cut-off frequencies of the filter. Thus when multiple frequencies are present, this may be thought of as taking the weighted summation of the cross-correlation functions at each frequency present and normalizing this by the product of the weighted summations of the autocorrelations across all frequencies present. The more narrow band the filter used, the more similar the time domain correlation and frequency domain magnitude coherence measures. As the filter encompasses a greater range of frequencies, measures from the two methods will increasingly deviate.
The low frequency time domain method employed by De Luca and colleagues  utilized a moving Hanning window as a low pass filter. The cut-off frequency of the filter is dependent on the time constant of the filter which is typically 400 ms [16, 21] but values up to 0.95 s have also been used . However different window lengths will modify the relationship between this time domain measure and the coherence function.
The effect of varying window length may be illustrated by obtaining an expression for the filter transfer function. The equation for the Hanning window is given as:
where τ is the length of the window. The discrete Fourier transform of this is given as (Kay, 1988):
In this section we demonstrate the relationship between time and frequency domain based methods to estimate the common modulating drive using empirical data. The methods are applied to data collected during isometric contractions of the First Dorsal Interosseous muscle at 20% of maximal effort. Two contractions where the activities of 4 and 5 MUs were identified yielded a total of 16 pairs of coactive MUs. The periods of concurrent activity of these MU pairs ranged between 30 s to 1 minute and were further divided into pairs of non-overlapping 10 s intervals resulting in a total of 50 pairs of 10 s long spike trains. Each method was applied to these spike train pairs and the correlation between the results yielded by the two methods were investigated as discussed below.
The time domain method was used to estimate low frequency common drive according to the method described by De Luca and colleagues . Three different time domain estimates were formed by smoothing the spike trains using Hanning windows of length 200, 400 and 800 ms respectively. These smoothed firing rate signals were then digitally high pass filtered with a low frequency cut-off of 0.75 Hz using a third order Butterworth filter to remove the mean bias discharge rates. The cross correlation coefficient function of these high pass filtered records was then obtained and the peak value of this function within ± 50 ms of the zero time lag was recorded and termed the time domain common drive coefficient.
The coherence analyses were performed in a similar manner to the procedure of Rosenberg and colleagues  for point process data. The spike trains were represented as binary pulse trains with ones corresponding to the firing times of the MU's and zeros comprising the remainder of the signal. Fourier transforms of these trains were obtained for each appropriately windowed section and then averaged according to equation (2). However where Rosenberg and colleagues  do not use overlapping or tapered data windows, we used overlapping, tapered Hanning windows of 2048 ms to optimize the variance and bias of the estimate. With any non-parametric spectral estimation technique, there is a trade-off between the variance and both the bias and resolution of the estimation. A window size of 2048 ms, gives a frequency resolution of 0.49 Hz, which is adequate to discriminate frequencies for our purposes. However, when analyzing 10 s of data using 2048 ms non-overlapping windows, only 5 different records are available and this small number of records will increase the variance of the estimate. Furthermore rectangular windows introduce an estimation bias due to the effect of their sidelobes. These concerns may be reduced by using the Welch periodogram method which uses tapered windows (to reduce spectral leakage and therefore the estimation bias) and overlapping windows (to increase the total number of windows and hence reduce the variance). The minimum variance for this method is obtained using an overlap of 62.5% . The frequency corresponding to the first zero-crossing of the Hanning filter was obtained according to equation (20) and the peak value of the coherence in the range between 0.75 Hz and this frequency was recorded.
A linear regression between the time domain common drive coefficients and corresponding frequency domain peak coherence values was performed to determine whether a linear relationship between the two indices existed. The regression r2 values are reported at a significance level of p < 0.05.
Results and discussion
As demonstrated in Figure 2, the cross-correlation peak can occur at lags slightly different than zero. A time delay or misalignment has been shown to introduce a bias into the coherence estimate that is proportional to the delay and coherence magnitude and inversely proportional to the FFT epoch duration . However for delays of the order of magnitude of ± 50 ms and FFT lengths of approximately 2 s, this type of bias is very small and unlikely to account for the observed differences between the time and frequency domain estimates.
Coherence estimates are typically formed from data records of around 1–5 minutes in length [4, 28, 29] or from pooled coherence measures of repeated trial measurements . This increases the number of non-overlapping windows in the calculation, thereby reducing the variance of the coherence estimate. Non-overlapping, rectangular windows are traditionally preferred due to the clear relationship with significance levels. Overlapping, tapered windows will allow coherence to be estimated from shorter data segments and parametric techniques, in particular multivariate autoregressive (MAR) methods are suggested for the analysis of very short duration data segments . When using short records of data (<5 s), the coherence estimates are likely to be significantly biased. However, the time domain method is more robust for such short data lengths and would therefore be preferred in these situations.
The time domain method uses a high pass filter to remove the mean bias from the smoothed signals, whereas the frequency domain coherence method simply subtracts the mean component of the signal prior to forming the estimate. Although similar, these two methods are not identical and may further explain some of the variation between the time and frequency domain techniques. A further possibility is to employ a low order polynomial detrending technique instead of high pass filtering or subtracting the mean. In general, a visual examination of the smoothed firing rate signals would indicate whether this would be necessary.
It is straight forward to quantify any time delay using the time domain technique. Although this is also possible with the frequency domain technique, this delay information is incorporated in the phase of the estimate and is therefore 2π periodic and would thus yield the same result for integer multiples of delay. For significant coherence present over a band of frequencies, Mima and colleagues  suggest a constant phase shift plus constant time delay regression model to compute time delays from coherence estimates. However for narrow band descending drives, the time domain technique provides a clearer estimate of any time delay.
It should be noted that the time domain technique may be generalized to cover any arbitrary frequency. This would necessitate bandpass filtering the signals to the frequency range of interest, removing the mean trends and then evaluating the cross-correlation function. Although this requires a priori knowledge of the drive bandwidth, this method would then be able to provide a single index to quantify a particular bandwidth.
Two separate bodies of literature offer techniques to estimate band limited common oscillatory input to motor neurons. These techniques are based in either the time or in the frequency domain. Indices derived from both measurement techniques are well correlated with each other and in the theoretical limit, the techniques are shown to be mathematically equivalent. However, for practical purposes there are a number of minor discrepancies which may favour the use of one particular method for a given application.
The time domain method offers greater resolution in time (the latency of the correlations are easily revealed) at the expense of the requirement of a priori knowledge of the bandwidth of the common modulating drive. On the other hand the frequency domain technique reveals information regarding the frequency distribution of the common modulating drives but it is more difficult to obtain accurate estimates of the coherence with short signal lengths as well as of the exact delay.
Time domain methods of estimation are preferred for short data segments and are well suited to quantifying the strength of a broad band drive into a single index. This proves useful in quantitative, comparitive analyses of the common behavior of MUs such as statistical tests that investigate the effects of aging or disease. Frequency domain measures tend to be more encompassing as they provide a complete description of all common oscillatory inputs. This facilitates qualitative analysis of distribution of coherence across frequencies and hence leads to a better understanding of the nature of the common inputs. They are well suited for estimation from large data segments, that may be assumed to be stationary, and are better able to quantify narrow ranges of descending inputs into a single index. Thus the selection of one technique over another should be dictated by the nature of the physiological question to be addressed.
The authors gratefully acknowledge the financial support of the Falk Medical Research Trust.
- Brown P: Cortical drives to human muscle: the Piper and related rhythms. Prog Neurobiol 2000, 60: 97-108. 10.1016/S0301-0082(99)00029-5View ArticlePubMedGoogle Scholar
- Grosse P, Cassidy MJ, Brown P: EEG-EMG, MEGEMG and EMG-EMG frequency analysis: physiological principles and clinical applications. Clin Neurophysiol 2002, 113: 1523-1531. 10.1016/S1388-2457(02)00223-7View ArticlePubMedGoogle Scholar
- Rosenberg JR, Amjad AM, Breeze P, Brillinger DR, Halliday DM: The Fourier approach to the identification of functional coupling between neuronal spike trains. Prog Biophys Mol Biol 1989, 53: 1-31. 10.1016/0079-6107(89)90004-7View ArticlePubMedGoogle Scholar
- Farmer SF, Bremner FD, Halliday DM, Rosenberg JR, Stephens JA: The frequency content of common synaptic inputs to motoneurones studied during voluntary isometric contraction in man. J Physiol 1993, 470: 127-155.PubMed CentralView ArticlePubMedGoogle Scholar
- Conway BA, Halliday DM, Farmer SF, Shahani U, Maas P, Weir AI, Rosenberg JR: Synchronization between motor cortex and spinal motoneuronal pool during the performance of a maintained motor task in man. J Physiol 1995, 489: 917-924.PubMed CentralView ArticlePubMedGoogle Scholar
- Salenius S, Portin K, Kajola M, Salmelin R, Hari R: Cortical control of human motoneuron firing during isometric contraction. J Neurophysiol 1997, 77: 3401-3405.PubMedGoogle Scholar
- Halliday DM, Conway BA, Farmer SF, Rosenberg JR: Using electroencephalography to study functional coupling between cortical activity and electromyograms during voluntary contractions in humans. Neurosci Lett 1998, 241: 5-8. 10.1016/S0304-3940(97)00964-6View ArticlePubMedGoogle Scholar
- Mima T, Hallett M: Corticomuscular coherence: a review. J Clin Neurophysiol 1999, 16: 501-511. 10.1097/00004691-199911000-00002View ArticlePubMedGoogle Scholar
- McAuley JH, Marsden CD: Physiological and pathological tremors and rhythmic central motor control. Brain 2000, 123: 1545-1567. 10.1093/brain/123.8.1545View ArticlePubMedGoogle Scholar
- Semmler JG, Kornatz KW, Enoka RM: Motor-unit coherence during isometric contractions is greater in a hand muscle of older adults. J Neurophysiol 2003, 90: 1346-1349.View ArticlePubMedGoogle Scholar
- Brown P, Marsden J, Defebvre L, Cassim F, Mazzone P, Oliviero A, Altibrandi MG, Di Lazzaro V, Limousin-Dowsey P, Fraix V, Odin P, Pollak P: Intermuscular coherence in Parkinson's disease: relationship to bradykinesia. Neuroreport 2001, 12: 2577-2581. 10.1097/00001756-200108080-00057View ArticlePubMedGoogle Scholar
- Salenius S, Avikainen S, Kaakkola S, Hari R, Brown P: Defective cortical drive to muscle in Parkinson's disease and its improvement with levodopa. Brain 2002, 125: 491-500. 10.1093/brain/awf042View ArticlePubMedGoogle Scholar
- Farmer SF, Sheean GL, Mayston MJ, Rothwell JC, Marsden CD, Conway BA, Halliday DM, Rosenberg JR, Stephens JA: Abnormal motor unit synchronization of antagonist muscles underlies pathological co-contraction in upper limb dystonia. Brain 1998, 121: 801-814. 10.1093/brain/121.5.801View ArticlePubMedGoogle Scholar
- Mima T, Toma K, Koshy B, Hallett M: Coherence between cortical and muscular activities after subcortical stroke. Stroke 2001, 32: 2597-2601.View ArticlePubMedGoogle Scholar
- Grosse P, Guerrini R, Parmeggiani L, Bonanni P, Pogosyan A, Brown P: Abnormal corticomuscular and intermuscular coupling in high-frequency rhythmic myoclonus. Brain 2003, 126: 326-342. 10.1093/brain/awg043View ArticlePubMedGoogle Scholar
- De Luca CJ, LeFever RS, McCue MP, Xenakis AP: Control scheme governing concurrently active human motor units during voluntary contractions. J Physiol 1982, 329: 129-142.PubMed CentralView ArticlePubMedGoogle Scholar
- Kamen G, Greenstein SS, De Luca CJ: Lateral dominance and motor unit firing behavior. Brain Res 1992, 576: 165-167. 10.1016/0006-8993(92)90625-JView ArticlePubMedGoogle Scholar
- Semmler JG, Nordstrom MA: Influence of handedness on motor unit discharge properties and force tremor. Exp RBraines 1995, 104: 115-125.Google Scholar
- Adam A, De Luca CJ, Erim Z: Hand dominance and motor unit firing behavior. J Neurophysiol 1998, 80: 1373-1382.PubMedGoogle Scholar
- Garland SJ, Miles TS: Control of motor units in human ffexor digitorum profundus under different proprioceptive conditions. J Physiol 1997, 502: 693-701. 10.1111/j.1469-7793.1997.693bj.xPubMed CentralView ArticlePubMedGoogle Scholar
- Semmler JG, Nordstrom MA, Wallace CJ: Relationship between motor unit short-term synchronization and common drive in human first dorsal interosseous muscle. Brain Res 1997, 767: 314-320. 10.1016/S0006-8993(97)00621-5View ArticlePubMedGoogle Scholar
- Patten C, Kamen G: Adaptations in motor unit discharge activity with force control training in young and older human adults. Eur J Appl Physiol 2000, 83: 128-143. 10.1007/s004210000271View ArticlePubMedGoogle Scholar
- Erim Z, Beg MF, Burke DT, de Luca CJ: Effects of aging on motor-unit control properties. J Neurophysiol 1999, 82: 2081-2091.PubMedGoogle Scholar
- Carter GC: Coherence and time delay estimation. Proc IEEE 1987, 75: 236-255.View ArticleGoogle Scholar
- Gardner WA: A unifying view of coherence in signal processing. Signal Processing 1992, 29: 113-140. 10.1016/0165-1684(92)90015-OView ArticleGoogle Scholar
- De Luca CJ, Erim Z: Common drive of motor units in regulation of muscle force. Trends Neurosci 1994, 17: 299-305. 10.1016/0166-2236(94)90064-7View ArticlePubMedGoogle Scholar
- Marple SL: Digital spectral analysis with applications. Englewood Cliffs, NJ: Prentice Hall; 1987.Google Scholar
- Kristeva-Feige R, Fritsch C, Timmer J, Lucking CH: Effects of attention and precision of exerted force on beta range EEG-EMG synchronization during a maintained motor contraction task. Clin Neurophysiol 2002, 113: 124-131. 10.1016/S1388-2457(01)00722-2View ArticlePubMedGoogle Scholar
- Gross J, Tass PA, Salenius S, Hari R, Freund HJ, Schnitzler A: Cortico-muscular synchronization during isometric muscle contraction in humans as revealed by magnetoencephalography. J Physiol 2000, 527: 623-631. 10.1111/j.1469-7793.2000.00623.xPubMed CentralView ArticlePubMedGoogle Scholar
- Amjad AM, Halliday DM, Rosenberg JR, Conway BA: An extended difference of coherence test for comparing and combining several independent coherence estimates: theory and application to the study of motor units and physiological tremor. J Neurosci Methods 1997, 73: 69-79. 10.1016/S0165-0270(96)02214-5View ArticlePubMedGoogle Scholar
- Cassidy MJ, Brown P: Hidden Markov based autoregressive analysis of stationary and non-stationary electrophysiological signals for functional coupling studies. J Neurosci Methods 2002, 116: 35-53. 10.1016/S0165-0270(02)00026-2View ArticlePubMedGoogle Scholar
- Mima T, Hallett M: Electroencephalographic analysis of cortico-muscular coherence: reference effect, volume conduction and generator mechanism. Clin Neurophysiol 1999, 110: 1892-1899. 10.1016/S1388-2457(99)00238-2View ArticlePubMedGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an open-access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.