active fibres. To ensure that muscles remained relaxed in the first two experiments in which spontaneously active afferents were investigated, the RMS filtered EMG signal was inspected visually and total surface area under the EMG curve was calculated and compared with baseline. Physiological conditions Apnoea. Subjects were asked to perform an inspiratory- capacity apnoea, in which subjects hold their breath at maximal lung volume against a closed glottis for 30–40 s [19]. Nociceptive stimulation. Two 23 G cannulae were inserted; one ,1 cm into the belly of the tibialis anterior muscle of the same leg form which the spindle afferent was being recorded; the other subdermally into the overlying skin ,5 mm away. Intramuscular bolus injection of 5% hypertonic saline (0.5 ml) caused sensations of deep, dull and diffuse pain, whereas subdermal injections of 5% hypertonic saline (0.2 ml) caused sensations of localised, sharp burning pain. After an appropriate baseline period was recorded, injections were made at unexpected times and in quasi-random order. Subjective pain level was evaluated on a visual-analog scale from 0 to 10, where 0 was described as ‘‘no pain’’ and 10 as ‘‘the worst pain ever experienced by the subject’’. Subjects indicated instantaneous pain level by a labelled potentiometer dial. Voluntary contraction. Subjects were asked to dorsiflex their ankle and gently apply force against force transducer and keep it at a constant level set to activate the muscle spindle at the rate of about 10 imp/s. Muscle force production was displayed on a computer monitor during voluntary muscle contractions for visual feedback to the subject. Analyses of pulse-wave effects on muscle spindle discharge rate The pulse wave effect was estimated both in the time and frequency domain [20]. In the time domain the effect was determined by R-wave triggered averaging. Instantaneous dis- charge rate traces were cut into segments between each of two consecutive R-wave peaks (RR intervals). Segments containing discharge rates that were more than three standard deviations (SD) different from the mean were removed to exclude artefacts and errors in automatic spike recognition. The remaining data were Fisher z-transformed to have zero mean and unit variance. If n is the number of data points in all artefact-free segments, this is xi{x achieved with xi~ , where x is the mean and s the standard s deviation over all data points x1, …, xn. The normalized data segments were then averaged with respect to the R-wave triggers. Note that due to normalization, the variance of the average segment must lie between 0 and 1. If the modulation of discharge rate is uncorrelated across segments, then the variance will be 0. If all segments are identical (fully correlated), then the variance of the average segment will be exactly 1. Otherwise, the variance will lie between 0 and 1. Hence, the variance of the R-triggered average quantifies the amount of the total variance explained by the effect of arterial pulsations (pulse wave modulation). To statistically test whether this modulation is caused by the pulse wave or by random events not related to cardiac rhythmicity, surrogate data were generated by circular permutation, i.e. each data segment was shifted a random number of samples before re- computing the average. A set of 1000 surrogates was generated. Because only the alignment of segments was changed, the variance of the surrogate averages represents the explained variance under the hypothesis. A straightforward statistical comparison between this distribution and the value derived from the Modulation of Human Muscle Spindle Discharge experimental data then permits formal testing of the hypothesis. This nonparametric technique of generating a distribution has a well-established role in statistics cf [21,22]. The pulse wave modulation was considered statistically significant if the explained variance was larger than 95% (p,0.05) of the surrogate distribution. In the frequency domain the pulse wave effect was assessed by means of the power spectral density of the discharge rate traces. The power spectral density is the Fourier transform of the autocorrelation function and can be used to test for modulations at the heart rate frequency. Power spectral density was estimated using Welch’s periodogram method using overlapping Hanning windows (window length, 7 s; overlap 3.5 s) [23]. To test for significant periodicity in the discharge rate, a set of 1000 surrogates were generated. To this end, the order of the spikes was locally permuted while leaving the inter-spike intervals the same, i.e. destroying possible temporal dependencies while leaving the mean and variability of discharge rate unchanged [24]. For each surrogate signal the power spectrum was computed and this surrogate distribution was then used to define the 95% confidence interval. Principal component analysis Principal component analysis (PCA) was employed to compare R-wave averages and power spectra across muscle spindles and to test for pulse wave modulation at the population level. PCA has commonly been used for unbiased statistical comparison of multivariate data across conditions [25–27]. Separate PCAs were conducted for R-wave averages and power spectra. The time and frequency axes were renormalized to compare discharge modu- lation across recordings with varying heart rates. The time axes of the averages were resampled using linear interpolation such that each RR interval contained an equal number of samples. The frequency axes of the power spectra were transformed from a linear scale to one relative to the heart rate frequency spanning 0 to 6 times the heart rate frequency. The broken-stick method was used to determine the number of principal components that were considered relevant for interpre- tation [28]. Of those principal components, the data were projected onto the eigenvector to determine the dominant features in the R-wave averages and power spectra. To determine the stability of those projections, the standard errors of each data point in the projection were estimated through 1000 bootstrap samples [29]. Bootstrap samples were generated using sampling-with- replacement and PCA was recalculated for each bootstrap sample. The ratio of the original projection to the bootstrap standard error is approximately equivalent to a z-score if the bootstrap distribution is normal [29]. Data points in the projection were considered stable if their z-scores were larger than 2.58 (p,0.01). Finally, eigenvector coefficients of significant principal compo- nents were compared across pain conditions and regressed against other variables to investigate correlations with the amount of discharge modulation [26,27]. Statistical analysis A combination of parametric and nonparametric statistics was used to describe the data and test for statistical significance [30]. To describe the prevalent effect characteristic for single muscle spindle afferents during normal conditions we primarily report median values. By contrast, we report mean values in the context of possible physiological effects because, regardless of skewed distributions, means better reflect functional consequences of converging inputs from a population of muscle spindles. The Wilcoxon matched-paired signed-ranks test was used to detect the PLoS ONE | www.plosone.org 3 April 2012 | Volume 7 | Issue 4 | e35091