Volume 21, Issue 1 , Pages 178-183, February 2011
Piper rhythm in the activation of the gastrocnemius medialis during running
Article Outline
- Abstract
- 1. Introduction
- 2. Methods
- 3. Results
- 4. Discussion
- 5. Conclusion
- Acknowledgments
- References
- Biography
- Copyright
Abstract
The presence of temporal rhythmicity in electromyographic (EMG) signals at frequencies of 35–60
Hz was initially noted by Piper (1907). This modulation and synchronization of motor unit activity is generally accepted to represent a centrally generated coding of motor commands. The purpose of this study was to resolve and quantify the Piper rhythm in the gastrocnemius medialis (GM) muscle during running. EMG was recorded from the GM of 14 female runners during 1-h treadmill runs. The average wavelet transform was computed for EMG from series of steps taken at 2
min intervals throughout the run. The total intensity across three wavelets (center frequencies: 170, 218 and 271
Hz) was computed and a histogram indicating the incidence peaks in this signal was generated for each subject. In order to rule out effects of the analysis process, the process was repeated using simulated EMG data. Autocorrelations of the histograms were used to extract the frequency of the peaks resulting in rhythmicity at 25–55
Hz. The ability to measure superimposed rhythmicity in EMG signals during dynamic tasks allows investigation of the role of aspects of central drive during movement. In particular, the changes in central control during dynamic activities can be examined with this approach.
Keywords: EMG, Running, Time–frequency analysis, Wavelet transform
1. Introduction
The presence of temporal rhythmicity in muscle activity has long been described during isokinetic tasks as well as during slow movements. In general, one can consider the EMG recorded from movements to comprise of muscle events (ME) and pacing events (PE). Muscle events represent a large portion of the energy in the signal and have envelopes containing low frequencies. These events correspond to the central control of the muscle providing an activation burst from onset to termination. Pacing events, on the other hand, have lower signal energy content, and represent an underlying rhythmicity contained within the muscle control signal.
The characteristics of the rhythms that occur during motor tasks and mechanisms proposed to explain them vary depending on the task. For example, a physiological tremor with a frequency of 8–10
Hz can be measured during position holding. Similarly, a discontinuity in movement occurring at a similar frequency is apparent during slow movements, although it is unclear that the mechanisms are common. The role of the spinal stretch reflex in rhythmic muscle activity has long been debated with data supporting (Hagbarth and Young, 1979, Sanes, 1985) and playing down (Jacks et al., 1988, Wessberg and Vallbo, 1995) the importance of its role. Instead the 8–10
Hz rhythms were attributed to pulsatile descending signals (Wessberg and Vallbo, 1995) and the corresponding frequency modulation of active motor units (Wessberg and Kakuda, 1999).
In addition, motor activity has been seen to synchronise at frequencies in the beta (15–30
Hz), and low gamma (30–60
Hz) bands (Brown, 2000). The specific case of modulated motor unit (MU) activity in the low gamma band is termed the Piper rhythm (Piper, 1907).
Rhythmicity in the Piper band has been shown to be the product of central drive, seemingly originating from the motor cortex (Brown et al., 1998). For example, neural activity in the motor cortex has been measured and shown to be coherent with activity in various muscles (Conway et al., 1995, Salenius et al., 1996, Salenius et al., 1997, Tecchio et al., 2006). Changes in central drive relating to fatigue have been shown to manifest as a decrease in the frequency of muscle rhythmicity in the Piper band. In addition, the features of muscle rhythmicity found in different pathologies may help to reveal the pathophysiology of the disorder (Brown and Marsden, 1996, McAuley et al., 2001), or conversely, clarify the pathway of the Piper rhythm (Brown, 1997). For instance, in un-medicated Parkinson’s patients, the 40–50
Hz frequency associated with the Piper rhythm is replaced with a 10
Hz signal indicating that the Piper pathway is influenced by activity in the pallidal projections (Brown, 1997).
Piper oscillations in motor output can be extracted from groups of single MU recordings (Hagbarth et al., 1983), or surface EMG signals (McAuley et al., 2001), or by measuring the resulting muscle vibrations (Brown, 1997, Hill, 1921). In our hands, the fine structure of the rhythmic pacing events in EMG activity becomes apparent in the wavelet transformed EMG signal. Our purpose was to test the hypothesis that during running, pacing events occur with a regular timing with respect to heelstrike. In addition, the uniformity of the pacing frequency among different individuals will be investigated. A simulation of muscle activity is used to evaluate the contribution of intrinsic rhythmicity stemming from the random firing of multiple motor units and to eliminate the possibility that the observed pacing events are the result of the wavelet transform method implemented.
2. Methods
Fifteen female recreational runners (age 32.4
±
8.7
years, mass 61.3
±
6.1
kg) participated in this study. Subjects gave their written informed consent in accordance with the university’s policy on research using human subjects. The protocol was approved by the Conjoint Heath Research Ethics Board at the University of Calgary.
2.1. Measurement protocol
Speed at ventilatory threshold (sVT) was determined during a preliminary session where respiratory gas exchange was monitored (Cosmed K4b2, Rome, Italy) during an incremental running test on a motorized treadmill (Quinton Instrument Co., Seattle, WA, USA). Belt speed was increased by 0.13
m/s every 2
min until the subjects reached ventilatory threshold (Fletcher et al., 2009). Ventilatory threshold was defined as the speed at which there was both a decrease in the fraction of expired CO2 and a non-linear increase in the VE/VCO2 slope (Reinhard et al., 1979, Skinner and McLellan, 1980, Wasserman et al., 1973). Two to four weeks later, each subject returned to the lab for a 1-h endurance running session during which they ran at approximately 95% of their speed at ventilatory threshold, sVT, while electromyographic (EMG) signals were recorded from the gastrocnemius medialis (GM) muscle of the right leg.
The subjects were recreational athletes and were instructed to maintain their usual activity levels, therefore, significant training adaptation was not expected to occur during the interval between sessions (Gaesser and Poole, 1986). At the selected speed of 95% of sVT the subjects were expected to be able to complete 1-h of running while experiencing substantial levels of fatigue. In two cases 95% sVT exceeded the subjects’ perceived 1-h sustainable leg turnover (i.e., for subjects with sVT >3.75
m/s), so the belt speed was decreased to 90% of their speed at ventilatory threshold.
In preparation for the run, the skin over the belly of the GM muscle was lightly abraded, cleaned with alcohol. Bipolar silver/silver chloride (Ag/AgCl) surface electrodes (Norotrode dual electrodes, Myotronics-Noromed Inc., Kent, WA, US) were placed on the belly of the GM muscle in alignment with the direction of the muscle fibres according to SENIAM recommendations (Seniam.org). This placement ensured that the electrodes were not located over the innervation zone of the muscle. All electrodes were placed by a single experimenter to ensure consistency throughout the study. A single reference electrode was secured to the tibial tuberosity. The electrodes and preamplifiers (Biovision, Wehrheim, Germany) were secured to the skin using medical tape (Cover-Roll stretch, BSN medical GnbH, Hamburg, Germany) to minimize movement artefact and to prevent the electrodes from losing surface contact due to sweating. EMG signals were preamplified (1000×) and bandpass filtered (10–500
Hz) during acquisition. Heelstrike was detected using an accelerometer attached externally to the heel of the right running shoe and recorded simultaneously with EMG at a sample rate of 2400
Hz. Thirty seconds of EMG and accelerometer data were recorded at 2
min intervals throughout the run. Each 30
s data recording is referred to as a “lap”, during which approximately 40 steps occurred.
2.2. Data analysis
All occurrences of heelstrike were determined by detecting the onset of the sharp deceleration associated with floor contact. EMG representing individual steps comprised 600
ms of data centered on heelstrike (1440 samples). This analysis window included all potential pre- and post-heelstrike muscle activation associated with the step.
A simulated version of the EMG was generated for every step. Simulated EMG was constructed by convolving a frequency modulated pulse train with an estimated MUAP.
2.2.1.1. Generation of a modulated pulse train
The 1440 sample raw EMG for each step was lowpass filtered (15
Hz) to determine the rough envelop of the activity. Dividing all of the envelope’s values by the peak value resulted in a normalized envelope function. The normalized envelope was divided into 20 bins (72 samples each). The value of the bin provided a scaling factor between 0 and 1, which was used to determine the number of random pulses that were inserted into the corresponding bin in the pulse train. The maximum number of pulses per bin was set to 9 (corresponding to 300
pulse/s). If fewer than 3 pulses were required in a bin they were omitted from the pulse train. The resultant pulse train had the same length (1440 samples) as the original raw EMG signal.
2.2.1.2. Estimation of MUAP shape
The bin containing the maximum amplitude in the envelope function was identified for each step. The Fourier transform was used to calculate the power spectrum for the raw EMG in the corresponding bin. The power spectra from each step recorded during one lap (i.e., 30
s series of steps) were averaged. The square root of the average power spectrum was multiplied by i (square root of −1) and was submitted to the inverse Fourier transform. This yielded a symmetric waveform that mimicked the shape of the MUAP from the original data.
A single estimated MUAP was generated for each lap and convolved with the pulse train generated for each individual step. This procedure resulted in unique simulated muscle activity for a given step (referred to as an “EMG step”) corresponding to each real “EMG step”.
This method was used to generate five sets of simulated data for each set of real data. Each set of resulting simulated data differed by the random pulse train used to generate any given EMG step.
2.2.2. Quantification of pacing events2.2.2.1. Wavelet transform
The wavelet transform of each real and simulated EMG step was performed using a set of 13 non-linearly scaled wavelets (von Tscharner, 2000). The center frequencies of the wavelets were: 7, 19, 38, 62, 92, 128, 170, 218, 271, 331, 395, 466 and 542
Hz. In general, the wavelet transform produces an EMG intensity pattern where time and frequency (center frequencies of the wavelet filters) are indicated on the abscissa and ordinate, respectively, and the greyscale represents the power of the transformed EMG signal.
The EMG average wavelet pattern was found corresponding to each lap. The average patterns from real EMG were visually inspected for movement artefact and noise content and all unacceptable average patterns were rejected along with their simulated counterparts. The total power in the band of the wavelets 7–9 (center frequencies 170–271
Hz) was found by summing the corresponding wavelet intensities for each time point. Low frequency bands were omitted from the selected range due to the long time resolutions of the wavelets, and high frequency bands were omitted because of the increased amounts of power reflecting signal noise.
2.2.2.2. Characterization of pacing events
All peaks in the summed wavelet intensity signal were identified. Peaks occurring at less than 1% of the maximum value of the signal were neglected (i.e., small oscillations around zero). A raster plot was generated for each subject where each row indicated the indices of the peaks in the intensity data from a single lap. A histogram was generated for each subject representing the number of occurrences of peaks at each sample time (1.67
ms bins). Separate raster plots and histograms were made for the real and five sets of simulated EMG data. Each plot was normalized to contain a total of 100 counts distributed across all bins.
The autocorrelation was taken of the six histograms generated for each subject (one real and five simulated EMG datasets). For each subject a net autocorrelation was calculated by subtracting the mean of the autocorrelations found from the five simulated datasets, from the autocorrelation found from the real data. The net autocorrelation was smoothed (lowpass 60
Hz) and the timing of the first three extrema (T1, T2 and T3) of the net autocorrelation were identified. The amplitude (i.e., range) of the oscillations in the net autocorrelation was determined by subtracting the mean values at the minima T1 and T3 from the value of the maximum at T2. Rhythmicity in the autocorrelation was deemed to be present if the resulting amplitude was above a threshold value.
The threshold value was selected to represent the minimum random oscillations expected in the autocorrelation and was determined using the autocorrelations of the five simulations for each individual subject. One of the five simulations was selected to represent the “test” data. The remaining four simulations were used to create the average autocorrelation of the simulated data. The net autocorrelation consisted of the autocorrelation of the new “test” data minus the new average autocorrelation of the simulated data. The amplitude of the resulting net autocorrelation was found. This process was repeated five times for each subject, each iteration using a different simulation as the “test” data. The threshold amplitude was calculated for each subject using the average of the amplitudes of the five resulting net autocorrelations. The corresponding threshold value was compared to each true net autocorrelation for each subject.
The pacing frequency for each subject with a net autocorrelation amplitude that was found to be above threshold was calculated using the equation: 2/(T2
+
(T3
−
T1)). A linear regression was used to find the relationship between pacing frequency and amplitude. The R2 value was determined.
3. Results
The average maximum aerobic speed of the 14 subjects included in this study was 3.3
±
0.4
m/s. For the test sessions the subjects ran at an average speed of 3.1
m/s (±0.3
m/s), representing approximately 94% (93.6
±
3.0%) of their speed at ventilatory threshold.
The wavelet patterns from the simulated data had a strong resemblance to those of the corresponding real data (Fig. 1C and D). In general, the traces of the simulated intensity bands were characterized by smaller, and more frequent oscillations (Fig. 1E and F). This was more apparent in the raster plot and histogram representation of the peaks in the intensity trace (Fig. 2).

Fig. 1.
Example of real (left) and simulated (right) EMG data from subject 4. Data is centered around heelstrike (HS) with 300
ms before and after. Raw signals (A and B), wavelet intensity patterns (C and D) and the summed total intensities across the wavelet bands 7–9 (E and F) are illustrated. The white lines in the wavelet intensity patterns indicate the bottom and top of wavelet bands 7–9.

Fig. 2.
Example of real (left) and simulated (right) EMG data from subject 8. The raster plots (A and B) indicate the timing of peaks (x-axis) in the total intensity traces from each trial (y-axis). The histograms (C and D) indicate the number of peaks across all trials per time bin.
For real EMG data, the peaks within a subject tended to have consistent timing across all laps. This resulted in a vertical alignment of dots in the raster plots, and clusters of bins with high counts in the histograms. The peaks in the simulated data were more randomly distributed. This was characterized by less structure in the raster plots, and a corresponding broad distribution of counts across all bins in the histogram. Examples of the histograms of real (top row) and simulated (middle row) data for five subjects are shown in Fig. 3, along with the corresponding autocorrelations (bottom row).

Fig. 3.
Examples of histograms from real (top row) and simulated (middle row) EMG data from five subjects. Resultant autocorrelations are indicated in the bottom row (upper grey – simulated data; upper black – real data; lower grey – net autocorrelation; lower black – smoothed net autocorrelation; black circles – first 2 minima (T1 and T3) and first maxima (T2) in the filtered net autocorrelation).
Of the 14 subjects, all were found to have net autocorrelations with amplitudes above threshold. The mean (±std) of the threshold amplitude for significant autocorrelations was 0.157
±
0.102. The detection of local extrema in the net autocorrelations resulted in pacing frequencies that varied from subject to subject within the range from 23.1 to 54.5
Hz (Table 1). The mean (±std) of the pacing frequencies was 40.6
±
10.7
Hz. The relationship between pacing frequency and amplitude showed a trend towards larger amplitudes for lower frequencies. The linear regression had an R2 value of 0.27.
Table 1. Pacing frequencies (F) and amplitudes (A) extracted from the net autocorrelations of each subject. Pacing amplitude is listed in units of the autocorrelation exceeding the threshold for each individual.
| Subject # | F (Hz) | A |
|---|---|---|
| 12 | 23.1 | 0.13 |
| 8 | 27.9 | 0.28 |
| 2 | 28.6 | 0.40 |
| 1 | 31.6 | 0.24 |
| 10 | 32.4 | 0.13 |
| 5 | 36.4 | 0.05 |
| 9 | 37.5 | 0.13 |
| 11 | 46.2 | 0.26 |
| 4 | 46.2 | 0.09 |
| 3 | 50.0 | 0.11 |
| 7 | 50.0 | 0.07 |
| 6 | 52.2 | 0.14 |
| 14 | 52.2 | 0.11 |
| 13 | 54.5 | 0.03 |
4. Discussion
There was a detectible rhythmicity within the bursts of EMG recorded from the gastrocnemius medialis during running. This pacing rhythm overlay the primary burst of muscle activity occurring after heelstrike corresponding to the generation of force required for forward propulsion during the stance phase. The frequency of the pacing rhythm was determined to range from approximately 25 to 55
Hz. This band of frequencies corresponds with Piper rhythms previously measured in EMG (Brown, 2000, Piper, 1907). It is also consistent with currently unpublished results from our group using the same method to measure rhythmicity during isometric contractions of the abductor pollicis brevis.
Some rhythmicity was observed in the autocorrelations derived from the simulated data. In some cases this rhythmicity was consistent across all simulations for a given subject and contained peaks in the autocorrelations that were shared with those from the real data (Fig. 3, s9). This suggested that aspects of the general shape of the primary muscle burst were conserved in the envelop used to generate the simulated data. In other cases the rhythmicity in the simulated data was not consistent between simulations (Fig. 3, s1). These fluctuations in the autocorrelation represent occurrences of synchronicity arising from the random positive and negative interference of the superimposed MUAPs used to generate the simulation. In both cases, subtracting the autocorrelation of the simulated data from that of the real data resulted in a net autocorrelation containing only aspects of the signal that represented the overlying pacing rhythm. This pacing rhythm was determined to be significant in all of the measured subjects.
Rhythms within the beta and low gamma (Piper) bands have been shown to be correlated to central descending drive (Brown, 2000). This suggests that the rhythms that were measured from the gastrocnemius muscle during running reflect aspects of central control generated within the brain. The fact that the pacing events aligned with respect to heelstrike for all measurements taken throughout a given run (evident in rasterplots – Fig. 2) suggests that the heelstrike event is transmitted to higher control centers and influences the phase of the rhythmicity. Not surprisingly, the pacing varied from subject to subject both in frequency and in amplitude. This reflects varying conditions, possibly including the shoe condition, effort level or state of fatigue being experienced by the individual. Data shown in the figures of a previous study (von Tscharner et al., 2003) seemed to indicate similar pacing events occurring in the tibialis anterior. In that study, the changes seemed to be related to the footwear condition. Fatigue has specifically been shown to influence the rhythmicity at the muscle, and with respect to the brain through measurements of corticomusclular coherence (Tecchio et al., 2006, Yang et al., 2009). The effect of these factors on the EMG pacing frequency during running will be investigated in future studies.
This study has some limitations, including the selection of the cut-off frequency of the lowpass filter used to extract the envelop of the real EMG signal for use in the generation of the simulated data. Despite the use of a 15
Hz cut-off, it is possible that a portion of the low frequency component of the original pacing events were conserved within the envelop. However, all net autocorrelations were found to have a significant amplitude, suggesting that the selected cut-off did allow the greater portion of the pacing signal to be omitted from the simulated dataset.
5. Conclusion
The ability to measure superimposed rhythmicity in EMG signals during dynamic tasks allows for the investigation of the role of aspects of central drive during movement. In particular, with this method, differences in fine pacing structures overlaying the muscle events during dynamic activity can be used to investigate changes of central control in the presence external or internal influences (e.g. fatigue or footwear conditions).
Acknowledgments
Funding was provided by the Alberta Heritage Foundation for Medical Research (AHFMR), the da Vinci Foundation and adidas AG. The sponsors had no role in the study design, in the collection, analysis and interpretation of data; in the writing of the manuscript; or in the decision to submit the manuscript for publication. The authors thank Dr. Brian MacIntosh, Jared Fletcher and Shane Esau for their contributions.
References
- . Muscle sounds in Parkinson’s disease. Lancet. 1997;349(9051):533–535
- . Cortical drives to human muscle: the Piper and related rhythms. Prog Neurobiol. 2000;60(1):97–108
- . Rhythmic cortical and muscle discharge in cortical myoclonus. Brain. 1996;119(4):1307–1316
- . Cortical correlate of the Piper rhythm in humans. J Neurophysiol. 1998;80(6):2911–2917
- Synchronization between motor cortex and spinal motoneuronal pool during the performance of a maintained motor task in man. J Physiol. 1995;489(Pt 3):917–924
- . Economy of running: beyond the measurement of oxygen uptake. J Appl Physiol. 2009;107(6):1918–1922
- . Lactate and ventilatory thresholds: disparity in time course of adaptations to training. J Appl Physiol. 1986;61(3):999–1004
- . Participation of the stretch reflex in human physiological tremor. Brain. 1979;102(3):509–526
- . The Piper rhythm—a phenomenon related to muscle resonance characteristics?. Acta Physiol Scand. 1983;117(2):263–271
- . The tetanic nature of the voluntary contraction in man. J Physiol. 1921;55(Suppl.):14–16
- . Instability in human forearm movements studied with feed-back-controlled electrical stimulation of muscles. J Physiol. 1988;402:443–461
- . Levodopa reversible loss of the Piper frequency oscillation component in Parkinson’s disease. J Neurol Neurosurg Psychiatry. 2001;70(4):471–475
- . Ueber den willkuerlichen Muskeltetanus. Pflueg Gesamte Physiol. 1907;119:301–338
- . Determination of anaerobic threshold by the ventilation equivalent in normal individuals. Respiration. 1979;38(1):36–42
- . Human cortical 40
Hz rhythm is closely related to EMG rhythmicity. Neurosci Lett. 1996;213(2):75–78 - . Cortical control of human motoneuron firing during isometric contraction. J Neurophysiol. 1997;77(6):3401–3405
- . Absence of enhanced physiological tremor in patients without muscle or cutaneous afferents. J Neurol Neurosurg Psychiatry. 1985;48(7):645–649
- . The transition from aerobic to anaerobic metabolism. Res Q Exerc Sport. 1980;51(1):234–248
- . Cortical short-term fatigue effects assessed via rhythmic brain–muscle coherence. Exp Brain Res. 2006;174(1):144–151
- . Intensity analysis in time–frequency space of surface myoelectric signals by wavelets of specified resolution. J Electromyogr Kinesiol. 2000;10(6):433–445
- . Changes in EMG signals for the muscle tibialis anterior while running barefoot or with shoes resolved by non-linearly scaled wavelets. J Biomech. 2003;36(8):1169–1176
- . Anaerobic threshold and respiratory gas exchange during exercise. J Appl Physiol. 1973;35(2):236–243
- . Single motor unit activity in relation to pulsatile motor output in human finger movements. J Physiol. 1999;517(1):273–285
- . Coding of pulsatile motor output by human muscle afferents during slow finger movements. J Physiol. 1995;485(1):271–282
- Weakening of functional corticomuscular coupling during muscle fatigue. Brain Res. 2009;1250:101–112

Lisa M. Stirling received her B.A.Sc. degree (with honours) in electrical engineering from the University of Toronto, Canada, in 2002. She completed her Ph.D. in medical sciences (biomedical engineering) at the University of Alberta, Canada, in 2007. Her graduate research focused on the development of control algorithms and functional electrical stimulation protocols for restoring standing and stepping after spinal cord injury. She is currently an adjunct assistant professor at the Human Performance Laboratory, University of Calgary, where she is pursuing her interests in the neural control of movement and the application of engineering approaches to the fields of rehabilitation and movement science.

Vinzenz von Tscharner was born in Switzerland 1947. He received his diploma in applied physics and mathematics 1974 and his Ph.D. degree in biophysics at the University of Basel, Switzerland. He was a post doctorate fellow at Oxford University, Dep. Biochemistry, England in 1978 and 1979, and a post doctorate fellow at Stanford University, Dep. Biochemistry, California, USA in 1998. He returned to the Biocenter in Basel 1981. He was then research affiliate at the Theodor Kocher Institute in Bern and specialized in signal transduction studying cellular responses related to cytokin binding. He became Adj. Assistant Professor (1997) and Adj. Associate Professor (2000) at the Human Performance Laboratory, University of Calgary. His main field of research is the signal propagation controlling movement patterns of humans. This involves biophysical/biomedical measurements and the analysis of sensory systems.

Patrick Kugler was born 1983 in Germany. He received his Diploma (German M.Sc.) with honours in Computer Science at the University of Erlangen-Nuremberg in 2009. During his studies he performed a half year internship at the Human Performance Lab at the University of Calgary, where he worked on his thesis on the classification of EMG signals. Currently, he is working as a Ph.D. student at the Pattern Recognition Lab at the University of Erlangen-Nuremberg. In his research he focuses on the application of algorithms from pattern recognition and machine learning to biomechanics and human motion.

Benno M. Nigg, was born in Switzerland, and studied nuclear physics at the world renowned ETH in Zurich, Switzerland. In 1971, he switched to Biomechanics. His goal was to improve individuals’ mobility and longevity through first, the study of forces impacting the lower body, and then the development of orthotics, running shoes, and exercise prescriptions that would enhance the quality of individuals’ lives. He joined the University of Calgary as the founder and first director of the Human performance Laboratory in 1981. Since his arrival, he has built a team of about 180 co-workers that have positioned the Human Performance Laboratory with the elite biomechanics programs in the world. He has published more than 280 articles in scientific journals and authored or edited eleven books. He has received numerous international awards, including the prestigious Olympic Order for recognition of this outstanding service and accomplishments for the Olympic Movement.
PII: S1050-6411(10)00096-9
doi:10.1016/j.jelekin.2010.06.007
© 2010 Elsevier Ltd. All rights reserved.
Volume 21, Issue 1 , Pages 178-183, February 2011
