# Imaging human cortical responses to intraneural microstimulation using magnetoencephalography

George C O'Neill, Roger H Watkins, Rochelle Ackerley, Eleanor L Barratt, Ayan Sengupta, Michael Asghar, Rosa M Sanchez Panchuelo, Matthew J Brookes, Paul M Glover, Johan W Wessberg, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

George C O'Neill, Roger H Watkins, Rochelle Ackerley, Eleanor L Barratt, Ayan Sengupta, et al. Imaging human cortical responses to intraneural microstimulation using magnetoencephalography. NeuroImage, 2019, 189, pp.329-340. 10.1016/j.neuroimage.2019.01.017. hal-02048371

HAL Id: hal-02048371

https://hal.science/hal-02048371

Submitted on 25 Feb 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Imaging human cortical responses to intraneural microstimulation using magnetoencephalography

<figure><img src="image_4.png" /></figure>

Check for
updates

George C. O'Neill<sup>a, *</sup>, Roger H. Watkins<sup>b</sup>, Rochelle Ackerley<sup>b, c</sup>, Eleanor L. Barratt<sup>a</sup>, Ayan Sengupta<sup>a</sup>, Michael Asghar<sup>a</sup>, Rosa Maria Sanchez Panchuelo<sup>a</sup>, Matthew J. Brookes<sup>a</sup>, Paul M. Glover<sup>a</sup>, Johan Wessberg<sup>b</sup>, Susan T. Francis<sup>a</sup>

## ARTICLE INFO

### Keywords:

Magnetoencephalography

Somatosensory cortex

Intraneural microstimulation

Single-unit

## ABSTRACT

The sensation of touch in the glabrous skin of the human hand is conveyed by thousands of fast-conducting mechanoreceptive afferents, which can be categorised into four distinct types. The spiking properties of these afferents in the periphery in response to varied tactile stimuli are well-characterised, but relatively little is known about the spatiotemporal properties of the neural representations of these different receptor types in the human cortex. Here, we use the novel methodological combination of single-unit intraneural microstimulation (INMS) with magnetoencephalography (MEG) to localise cortical representations of individual touch afferents in humans, by measuring the extracranial magnetic fields from neural currents. We found that by assessing the modulation of the beta (13–30 Hz) rhythm during single-unit INMS, significant changes in oscillatory amplitude occur in the contralateral primary somatosensory cortex within and across a group of fast adapting type I mechanoreceptive afferents, which corresponded well to the induced response from matched vibrotactile stimulation. Combining the spatiotemporal specificity of MEG with the selective single-unit stimulation of INMS enables the interrogation of the central representations of different aspects of tactile afferent signalling within the human cortices. The fundamental finding that single-unit INMS ERD responses are robust and consistent with natural somatosensory stimuli will permit us to more dynamically probe the central nervous system responses in humans, to address questions about the processing of touch from the different classes of mechanoreceptive afferents and the effects of varying the stimulus frequency and patterning.

## 1. Introduction

In humans, our sense of touch is conveyed by an array of discrete low-threshold mechanoreceptive afferents distributed throughout the skin, and haptic sensations are derived from a combination of activity across stimulated receptors (Johnson, 2001; Saal and Bensmaia, 2014). The sensory afferents from these receptors project to the somatosensory cortices, via the thalamus, enabling the decoding of the nature and location of the stimulus evoking the sensation. Our somatosensory abilities support a wide range of systems in the body, from integrating with the motor cortex for force feedback and fine-motor control (Ackerley and Kavounoudias, 2015), to emotional impact, providing feedback from social situations (McGlone et al., 2014).

Much of our knowledge about human mechanoreceptive afferents and

their electrophysiological properties stems from microneurography studies (Vallbo and Hagbarth, 1968), where an electrode is placed in a peripheral nerve to register axonal impulses. Microneurography allows the characterisation of *single* mechanoreceptive afferent responses when the associated receptive field in the skin is excited by touch. Four distinct low-threshold mechanoreceptor classes exist in glabrous (non-hairy) skin, each with their own distinct spiking pattern in response to somatosensory stimuli and spatial representations within the skin (Vallbo et al., 1984a). In addition to registering impulses, it is possible to pass a small current (of approximately a few microamperes) down the electrode to elicit artificial impulses in the same mechanoreceptive afferent, in the absence of contact to the skin. Here, a small, corresponding simulated tactile sensation can be felt, i.e. the projected receptive field (at the same site as the physiological receptive field), which is thought to represent the sensation conveyed by

<footnote><sup>a</sup> Sir Peter Mansfield Imaging Centre, School of Physics and Astronomy, University of Nottingham, Nottingham, UK</footnote>

<footnote><sup>b</sup> Department of Physiology, University of Gothenburg, Gothenburg, Sweden</footnote>

<footnote><sup>c</sup> Aix Marseille Univ, CNRS, LNSC (Laboratoire de Neurosciences Sensorielles et Cognitives – UMR 7260), Marseille, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: george.oneill@nottingham.ac.uk (G.C. O'Neill).</footnote>
these individual mechanoreceptors. This process has been termed single-unit intraneural microstimulation (INMS; Torebjörk and Ochoa, 1980; Vallbo, 1981) and allows precise impulse patterns to be evoked in individual nerve fibres, something traditional transdermal or intraneural stimulation cannot provide. The use of low current in INMS allows single afferents to be stimulated and avoids both the recruitment of numerous afferents and involuntary movements that may be induced by efferent activation when using higher currents, such as with median nerve stimulation. Single-unit INMS has shown that physiologically-defined mechanoreceptor types have distinct perceptual correlates in humans (Ochoa and Torebjörk, 1983; Vallbo et al., 1984b). Fast adapting type I (FA1) and type II (FA2) afferents elicit a vibration or fluttering sensation when stimulated, whilst slowly adapting (SA) type I (SA1) give rise to a pressure or pulling sensation. Slowly-adapting type II (SA2) afferents do not typically produce a clear, identifiable sensation. The type I units, which are more numerous (Johansson and Vallbo, 1979), provide small, point-like sensations. The use of single-unit INMS has the advantage over traditional tactile stimuli of being able to produce a defined and controllable pattern of firing within a single afferent and single afferent class. This allows the precise examination of how input patterns from peripheral mechanoreceptive afferents are perceived and represented in the central nervous system.

Despite knowledge of the peripheral and perceptive properties of single mechanoreceptive afferents, little is known about how their responses to stimulation manifest in the human brain. Using single-unit INMS, it is possible to probe the sense of touch at the "quantal" level of individual afferents, something which cannot be achieved for any of the other senses. There have been relatively few studies in which neuroimaging has been combined with single-unit INMS to investigate the neural effects of selective mechanoreceptive afferent stimulation. A study by Kelly et al. (1997) used high density electroencephalography (EEG) over the primary somatosensory cortex (S1) and showed that stimulating single FA1 or SA1 units projecting to the hand gave characteristic peaks in power at the same frequency at which the units were stimulated. Only two studies have been performed using functional magnetic resonance imaging (fMRI) to map the spatial profile of the cortical responses to single-unit INMS (Sanchez Panchuelo et al., 2016; Trulsson et al., 2001). The earlier study, performed at 3.0 T, showed that individual type I units have corresponding neural activations in S1 and the secondary somatosensory cortex (S2), and the responses were in good spatial agreement with activations from vibrotactile stimulation (Trulsson et al., 2001). The more recent study was performed at ultra-high field (7.0 T) using high spatial resolution to show that functional activation patterns of units located on the hand digits were spatially well-localised within the expected digit regions, ascertained from maps of digit somatotopy obtained from vibrotactile stimulation to the skin of the fingertips (Sanchez Panchuelo et al., 2016).

Magnetoencephalography (MEG) measures the extracranial magnetic field associated with intracellular current flow (Hämäläinen et al., 1993). MEG provides unique spatiotemporal resolution, offering vastly improved temporal resolution over fMRI (which is limited by the hæmodynamic lag of the BOLD effect, and it being an indirect measure of neural activity) and an improved spatial resolution compared to EEG (as magnetic fields are not distorted by the skull to the same extent as electric potentials). In this paper, we performed single-unit INMS, using a recently developed system designed for compatibility with MEG (Glover et al., 2017). We aimed to demonstrate that this approach can localise the cortical response to single-unit INMS and capture the spatiotemporal dynamics of responses to simulation of individual mechanoreceptors, providing the feasibility of performing future studies to probe how the dynamics of different patterns of activity in these mechanoreceptors are represented cortically.

## 2. Methods

### 2.1. *Experimental procedure*

A total of 9 healthy participants (age 24–68 years, 5 female)

volunteered to undergo microneurography/single-unit INMS and simultaneous MEG recording. Each participant was given detailed information about the procedure and provided their written informed consent. Ethical approval was granted by the University of Nottingham Medical School Ethics Committee, with all procedures conducted in line with the Declaration of Helsinki.

Experimental sessions consisted of 3 phases: 1) microneurography for the characterisation of a single mechanoreceptive afferent (Vallbo and Hagbarth, 1968); 2) assessment of the sensation to single-unit INMS; 3) MEG recordings during INMS. These steps were repeated and ~2–3 individual mechanoreceptive afferents were typically stimulated during MEG data collection in a single experimental session. The participant was seated comfortably with their head just below the MEG helmet, and were instructed to sit as still as possible with their arm/hand held still using a vacuum cushion to minimise any slight movement of the arm/hand, which would result in slight movements of the electrode and loss of the unit. A high-impedance, insulated tungsten recording/stimulating electrode (15 mm length, shaft diameter 200 µm, tip diameter ~5 µm; FHC, Bowdoin, ME) was then inserted percutaneously into the left median nerve, approximately 30 mm proximal from the wrist crease. A second uninsulated reference electrode was placed subcutaneously 50 mm away from the first and both were connected to an INMS system specifically designed to operate without generating measurable magnetic interference (Glover et al., 2017). A diagram of the electrode setup is shown in Fig. 1.

For experimental phase 1, the INMS system was in 'record mode', and a real-time visual and audio readout of the potential difference between the two electrodes was output to the experimenter within the magnetically shielded room. The location of mechanoreceptive afferents was determined by gentle touch across the hand, delivered by the experimenter, between minor adjustments of electrode position, until a unitary recording was detected. The precise location of the unit receptive field on

<figure><img src="image_5.png" /><figcaption>**Fig. 1.** Experimental setup for microneurography and single-unit micro-stimulation. A) A simplified diagram showing how a single mechanoreceptive afferent is interrogated. A recording/stimulating electrode is placed in the median nerve to register activity from a single nerve fibre. If mechanical stimulation is applied to the afferent's receptive field, the electrode in 'record mode' will register the spiking pattern being sent to the central nervous system. In 'stimulate mode', passing a small current down the electrode will stimulate the mechanoreceptive afferent, eliciting a tactile sensation at the location of the receptive field in the skin. B) The experimental setup on a participant in the study.</figcaption></figure>
the skin and force activation threshold of the unit was then identified using von Frey monofilaments and marked with ink. Once the mechanoreceptive afferent was identified, the response pattern and spatial extent on the skin was examined to identify the receptor type as SA1, FA1, SA2 or FA2 (Vallbo et al., 1984a). The INMS system was then switched to 'stimulate mode', to assess the sensation to INMS (experimental phase 2). A train of constant current pulses (pulse width 200 µs, amplitude <10 µA) was delivered at 60 Hz for 1 s. The trains of pulses were delivered repetitively and the current was increased gradually (resolution 0.1 µA), until the participant felt a sensation in their hand, which was time locked to the stimulation. The participant was asked to describe the nature and location of the sensation; if there was an exact correspondence with the physiological properties of the unit (e.g. location, quality, size) with only a single unit stimulated at that current, the investigation proceeded to the MEG data collection (experimental phase 3). The paradigm during the MEG recording consisted of 80 trials of a 1 s long train of pulses (as described above), with a base inter-trial-interval of 10 ± 1 s to jitter the stimulus. After every 10 trials, the quality and intensity of the sensation was confirmed verbally with the participant, as small movements of the electrode can diminish the intensity of the sensation, presumably reflecting failure to generate impulses with each stimulation pulse. In the case where a participant no longer felt the sensation, the pulse amplitude was increased until a similar quality of sensation was again felt again in the same location (Torebjörk et al., 1987). If a comparable sensation was achieved by increasing the stimulation intensity, the MEG protocol was resumed, but if comparable sensation intensity was not recovered, or there was evidence of a change in the unit stimulated (e.g. its quality/location), the MEG recording was terminated.

Immediately after INMS, the participant underwent a short experiment, where mechanical vibrotactile stimulation was applied to the skin at each microstimulated unit's receptive field location to determine the amplitude at which the vibrotactile stimulation matched the INMS sensation. A small piezoelectric vibrotactile stimulator (Dancer Design, St. Helens, UK) controlled by an Arduino microcontroller (http:// arduino.cc) was used to deliver a vibrotactile stimulation at 60 Hz. The amplitude was adjusted until the sensation felt similar in intensity to that from the INMS. These amplitudes were noted and used in a subsequent MEG experiment in which vibrotactile stimulation at each micro-stimulated unit's receptive field was carried out with identical timings to the INMS paradigm. This vibrotactile/MEG follow-up session was performed at a later date to allow for the maximum number of units to be stimulated within each INMS experimental session.

## 2.2. Data acquisition

MEG data were recorded using a 275 channel MEG system (CTF; Coquitlam, BC) in synthetic 3rd order gradiometer configuration, with a sampling rate of 1200 Hz, and a hardware anti-aliasing low-pass filter of 300 Hz applied prior to signal digitisation. An analogue trace of the INMS stimulation pulses was recorded concurrently to allow for offline synchronisation of the neuromagnetic recordings. All participants wore a bespoke flexible headcast (details of the manufacturing procedure can be found in Liuzzi et al., 2017; Meyer et al., 2017), designed to minimise head movement within the MEG system. The headcasts include cavities designed to fit three head position indicator (HPI) coils, which were periodically energised during the experiment, to determine their position within the MEG sensor array. The HPI coil locations relative to the brain were defined during headcast manufacture, meaning potential source localisation errors typically introduced when misaligning the sensors and the anatomy by other means were minimised (Troebinger et al., 2014). For a single participant, the headcast failed to fit and so an alternative method to coregister the brain to the sensors was employed. In this case the HPI indicators were attached to the head at the nasion and bilateral preauricular locations. A measurement of coil locations relative to the scalp was recorded using a 3D digitiser (Polhemus; Colchester, VT).

Coregistration was achieved by matching the digitised surface with a surface extracted for the anatomical image. Anatomical images were acquired using a 3.0 T MR System (Philips Healthcare; Best, Netherlands) using an MPRAGE sequence (1 mm³ isotropic resolution).

## 2.3. Data analysis

MEG data were visually inspected for artefacts. If an epoch contained excessive artefacts (for example generated by muscle or eye movement, eye blinks, or from SQUID resets), it was omitted from further analysis. The sensor data were then band-pass filtered to the beta band (13–30 Hz), as this is one of the oscillatory components modulated by somatosensory stimulation (Hari and Salmelin, 1997; Cheyne, 2013). MR images were segmented and cortical surfaces were reconstructed using FreeSurfer 5.3 (Dale et al., 1999; Fischl et al., 1999). Further segmentation of the scalp, as well as inner and outer skull surfaces, were generated using modified routines from FieldTrip (Oostenveld et al., 2011). The surface which defined the grey/white matter boundary was decimated to 12500 vertices per hemisphere and each vertex was used as a location for source reconstruction. Source reconstruction at each vertex was achieved with a Linearly Constrained Minimum Variance (LCMV) beamformer, with the data covariance matrix generated from the entire duration of the experiment (Brookes et al., 2008). Dipole approximation for forward modelling was calculated using a 3-shell boundary element model (BEM; Stenroos and Nummenmaa, 2016), with the optimal dipole orientation for maximal signal variance derived using an eigenvalue decomposition approach (Sekihara et al., 2004). With beamformer weights for a given vertex, w, generated, they are employed to quantify changes in power across the cortex. For every location, an activation index, A, was generated to quantify changes in source power during the stimulation (0–1 s after stimulus onset) compared to a baseline period (8.9–9.9 s after stimulus onset). For the i'th trial, A was calculated at each location using Equation (1),

$$A(i) = \frac{\mathbf{w}C_s(i)\mathbf{w}^T - \mathbf{w}C_b(i)\mathbf{w}^T}{\mathbf{w}C_s(i)\mathbf{w}^T + \mathbf{w}C_b(i)\mathbf{w}^T}, \quad (1)$$

where $C_s(i)$ and $C_b(i)$ are the single-trial sensor-level covariance matrices which correspond to the stimulation and baseline epochs respectively. Activation indices were averaged across trials to generate a final activation image for each experiment. The activation index is scaled between -1 and 1, with 0 denoting no modulation in power between the two epochs. To allow for a group level analysis, all images were transformed to a standard cortical surface derived from the MNI-305 template brain. In addition, to interpret source localizations multimodal parcellations published in the Glasser atlas (Glasser et al., 2016) were used. The Glasser atlas is a multi-modal parcellation of the brain derived from Human Connectome Project (HCP-MMP1.0) and was converted using the HCP Workbench software to fit the MNI-305 brain. In addition, we also compared our results to previously derived functional markers. Here we used probabilistic atlas of the hand area generated from fMRI responses to vibrotactile somatosensory stimulation of the left hand (Sengupta et al., 2018).

## 2.4. Statistical inference of significant power changes

To identify which locations corresponded to a significant change in power during stimulation (INMS or vibrotactile), we employed non-parametric permutation testing (Maris and Oostenveld, 2007; Nichols and Holmes, 2002; Pantazis et al., 2005; Singh et al., 2003) at both the individual subject and group levels. At the individual subject level, our null hypothesis was that the oscillatory power measured in the stimulation and the baseline epochs within a trial were identical. Therefore, if we were to randomly shuffle the labelling of the covariance matrices within trials, prior to generating A in Equation (1), then averaging across trials to form a null averaged activation index, A, would not give a significantly
different result to the genuine averaged activation index, $\bar{A}$. 2000 permutaions of $\bar{A}$ were generated to create a null distribution for each vertex location. To correct for multiple comparisons, the threshold free cluster enhancement score (TFCE; Smith and Nichols, 2009) was calculated for every vertex on each iteration of $\bar{A}$. For a given location $r$, the TFCE score was calculated as

$$TFCE(\bar{A}(\mathbf{r})) = \int_{0}^{h_{max}} e(\bar{A}(\mathbf{r}), h)^{0.5} h^2 dh \quad (2)$$

where $e(\bar{A}(\mathbf{r}), h)$ is a function which returns the spatial extent of a cluster at $\mathbf{r}$, should $\bar{A}(\mathbf{r})$ exceed a threshold value $h$, and 0 otherwise. Positive and negative clusters were treated as non-contiguous, and the maximum absolute value of TFCE on each iteration was used to generate a global null. Any vertices where the absolute TFCE of $\bar{A}$ exceeded the 95th percentile of this null were considered significant.

At the group level, our null hypothesis was that the activation maps simply represented zero-mean noise, so randomly multiplying any of the experiments' trial averaged activation images by $-1$ and then group averaging to generate a null activation map, $\bar{A}$, would not give a significantly different result to the group-averaged image $\bar{A}$. All $n!$ possible permutations (where $n$ represents the number of images in the group) were calculated and similarly to the group level analysis, multiple comparisons correction was handled by finding the critical TFCE threshold. By performing statistical analyses at both the group and individual level, we can determine whether changes in oscillatory power which are observed at the group level changes are potentially a) applicable markers for studying changes in individual experiments and b) generate a probabilistic map of where the cortical source of stimulation to the hand area localises to.

# 3. Results

## 3.1. Electrophysiology

Across the cohort of participants, 39 single-unit recordings were identified. 19 of these provided percepts that matched the unit physiology, indicating successful single-unit INMS (Torebjörk et al., 1987) and cortical responses to stimulation were recorded during MEG. An excerpt from a single-unit recording from an FA1 can be seen in Fig. 2A. Here, we observed high spiking rates (occasionally exceeding 200 impulses $s^{-1}$) when a von Frey monofilament was pressed into the receptive field location in the skin. Data from 3 units were discarded due to too few MEG trials being recorded before losing the unit sensation. The properties of each of the 16 units with full MEG recording protocols can be found in Table 1, and their locations on the hand in Fig. 2B. 10 units were identified as mechanoreceptor type FA1, 5 as SA1, and 1 as FA2. In the results presented here, we primarily focused on the 10 FA1 units, given that these formed the largest group, and their perceived sensation — vibration/ flutter — is similar to that of vibrotactile stimulation. For completeness, the localisation results of the other types, 5 SA1 units and a single FA2 unit, can be found in the Supplementary material.

## 3.2. MEG

Fig. 3 identifies where in the brain changes in beta power during single unit INMS across the 10 FA1 units were located. Fig. 3A shows a conjunction map of the number of individual units that had a significant ($p_{corrected} < 0.05$) change in power. We observed a large degree of overlap across the units in the contralateral primary somatosensory and motor (M1) cortices, with 8/10 units in spatial agreement; the other 2 units showed no significant changes in power at any spatial location in the brain (see Supplementary material). In the ipsilateral S1, we observed overlap of activations in 2/10 units. Fig. 3B shows the map of group-averaged activation index ($\bar{A}_{INMS}$), thresholded to show only the

<figure><img src="image_6.png" /><figcaption>Fig. 2. A) Excerpt from a single-unit peripheral nerve recording from a fast adapting type 1 (FA1) afferent (unit 6a), showing the unitary responses to mechanical stimulation of the receptive field in the skin whilst in record mode. Upon single-unit INMS, the subject felt a 3 mm spot defined to have a 'prickle' sensation. The location of the receptive field can be found in panel B. (B) Receptive field locations in all 16 units in the right hand which were micro-stimulated whilst in the MEG scanner, with adjacent numbers corresponding to the Unit ID in Table 1. The shapes correspond to the receptor type; dark blue circles represent FA1 afferents (vibration sensation in a very small area); blue squares show SA1 afferents (pressure or pulling in a small area), and the cyan triangle represents an FA2 afferent (diffuse vibration sensation).</figcaption></figure>

significant changes in beta power in response to INMS at the group level. This map shows that there was a strong decrease in oscillatory power, which was (much like the high overlap region at the individual level) situated primarily over the contralateral S1 and extended into M1. However, it is important to note that these three areas with the largest reductions in beta power could be in fact a single diffuse reconstructed source that is represented in multiple locations when constrained to a cortical surface. To help clarify where the origin of this source may be, the location of maximal reduction is represented as a magenta dot and is situated in S1. Next, we compared our results to previous findings; Fig. 3C is a zoomed representation of the contralateral sensorimotor area, with the Brodmann areas from the Glasser atlas (Glasser et al., 2016) and fMRI derived hand area (Sengupta et al., 2018) in grey and black, respectively. Here, we see that the area of with the largest spatial overlap (yellow) is primarily situated within the hand area. The magenta dot from Fig. 3B has also been overlaid and is situated within the high spatial overlap area, implying the largest reduction in beta power is a common feature across units, rather than being driven by a small subset. Fig. 3D is a zoomed representation of the findings in Fig. 3B. Again, the Glasser atlas is overlaid, but in green for better contrast, and the hand area remains in black. The location of the peak reduction in beta power, represented by the magenta dot, is situated on the borders of Brodmann areas 1 and 2, confirming that the peak modulation of power occurred in S1, rather than in M1. Given that the magenta marker was situated at the largest group effect, and the highest amount of overlap across individual subjects, this makes for a reasonable candidate location to place a virtual
<caption>**Table 1**</caption>

Metadata on the 16 microstimulated units recorded in the MEG scanner. Unit identifiers are based on the order of collection. Note all locations are on the palmar side of the hand unless otherwise noted (D denotes digit number).

<table><thead><tr><th>Participant ID</th><th>Unit ID</th><th>Type</th><th>Location</th><th>Sensation</th><th>Threshold current (uA)</th><th>Perceptive Field Size (mm)</th></tr></thead><tbody><tr><td rowspan="2">1</td><td>1a</td><td>FA1</td><td>Medial D3</td><td>Vibration</td><td>2.2</td><td>Not recorded</td></tr><tr><td>1b</td><td>SA1</td><td>Proximal D4</td><td>Pulling</td><td>0.9</td><td>0.1</td></tr><tr><td rowspan="2">2</td><td rowspan="2">2a<sup>a</sup></td><td>FA1</td><td>Distal D3</td><td>Tingle</td><td>2.2</td><td>2</td></tr><tr><td>3a</td><td>Palm</td><td>Pulling</td><td>1.4</td><td>0.1</td></tr><tr><td rowspan="5">3</td><td>3b</td><td>SA1</td><td>Palm</td><td>Pulling</td><td>2.8</td><td>1</td></tr><tr><td>3c</td><td>FA1</td><td>Palm</td><td>Vibration</td><td>3.7</td><td>1</td></tr><tr><td>3d</td><td>FA1</td><td>Distal D3</td><td>Flutter</td><td>3.3</td><td>1</td></tr><tr><td>3e</td><td>FA1</td><td>Dorsal Medial D3</td><td>Flutter</td><td>0.9</td><td>1</td></tr><tr><td>4a</td><td>FA1</td><td>Distal D4</td><td>Vibration</td><td>3.4</td><td>1</td></tr><tr><td rowspan="5">4</td><td>4b</td><td>FA1</td><td>Distal D4</td><td>Vibration</td><td>4.2</td><td>2</td></tr><tr><td>4c</td><td>FA2</td><td>Distal D4</td><td>Tapping</td><td>4.4</td><td>7</td></tr><tr><td>4d</td><td>SA1</td><td>Medial D4</td><td>Point-like pressure</td><td>2.4</td><td>0.1</td></tr><tr><td>4e</td><td>SA1</td><td>Medial D4</td><td>Pressure</td><td>3.2</td><td>3</td></tr><tr><td>4f</td><td>FA1</td><td>Distal D4</td><td>Vibration</td><td>2.0</td><td>3</td></tr><tr><td>5</td><td>5a</td><td>FA1</td><td>Distal D1</td><td>Tingle</td><td>4.2</td><td>2</td></tr><tr><td>6</td><td>6a</td><td>FA1</td><td>Palm</td><td>Prickle</td><td>1.1</td><td>3</td></tr></tbody></table>

INMS
Beta: 0-1 after stimulaiton onset

<figure><img src="image_7.png" /><figcaption>**Fig. 3.** Results of the INMS from the group of 10 FA1 units studied overlaid onto inflated cortical surfaces, gyri are represented as light grey surfaces sulci are dark grey. A) Results of the individual unit analysis in beta band (13–30 Hz) power for a given location, represented as a conjunction map, showing the number of FA1 units that exhibit significant changes. B) Group-averaged activation image of the 10 FA1 units at a statistical threshold (p < 0.05), showing microsimulation elicits a reduction in beta oscillatory power over the contralateral sensorimotor area, with the peak reduction in beta power represented by the magenta dot) occurring in the postcentral gyrus. C) Zoomed portion of panel A, with the Brodmann areas as defined by the Glasser atlas (Glasser et al., 2016) overlaid in grey and hand area derived from a functional atlas (Sengupta et al., 2018) in black. Here the region representing the largest overlap across units (yellow; N = 8) is contained primarily in the postcentral gyrus and within the predefined hand area. D) Zoomed portion of panel B, with both the Glasser atlas (this time in green) and the hand area overlaid. The peak reduction of beta power represented as the magenta dot) is located in postcentral gyrus. E) A time-frequency spectrogram of the peak location from panels B, C and D, showing a distinct event-related behaviour in the beta band which extends beyond the initial 1 s of stimulation. Interactive versions of the cortical plots can be found at http://georgeoneill.github.io.</figcaption></figure>

electrode to explore the other spectral and temporal properties during the experiment. The results of such are shown in Fig. 3E, which shows a group-averaged time-frequency spectrogram of INMS at this location. The spectrogram shows that it is the beta oscillatory band that showed

the largest fractional change in power during stimulation, where we observed event-related desynchronisation (ERD). After the stimulation ceased (t = 1), post-event rebound (PER) occurred, with an overshoot in beta power at t = 2–3 s, before returning to baseline. Performing the

<footnote><sup>a</sup> Unit acquired without headcast.</footnote>
same localisation analysis on beta power in the PER epoch showed this to be significant at the group level, but individual testing showed that this was only significant for three units, all from the same participant (see Supplementary material). We also found similar temporal properties in the t = 0–1 s and t = 2–3 s epochs in the alpha band, but few units showed individual level significance and the locations of group level significance had a reduced activation index magnitude compared to beta band activity (see Supplementary material).

The equivalent statistical maps from the vibrotactile experiments are shown in Fig. 4. In Fig. 4A, we observed similar results at the individual level, as depicted in Fig. 3A, with spatial overlap between the vibrotactile stimuli applied at the receptive field location for each unit and INMS in contralateral S1 in 9/10 corresponding vibrotactile stimuli. In the ipsilateral hemisphere we saw a greater number of overlapping responses between vibrotactile stimulation and INMS (3/10) with the area of overlap more parietal than contralateral S1. Fig. 4B shows the group-averaged activation index ($\bar{A}_{Vibrotactile}$), again with a statistical threshold. A significant decrease in beta band oscillatory power was found over the contralateral postcentral gyrus during vibrotactile stimulation, agreeing with the findings of INMS. Furthermore, the magenta dot, which represents the location of the maximum reduction of beta power during stimulation, was located in the exact same vertex as the INMS experiments (i.e. 0.0 mm distance between INMS and vibrotactile

peaks). Fig. 4C is a zoomed-in representation of the plot in Fig. 4A, which again overlays the Glasser Atlas and the fMRI derived hand area. We observed that our regions exhibiting the highest level of significance across the subjects occurred within the hand area, and within the Brodmann Areas that corresponded to the postcentral gyrus. Fig. 4D is a zoomed-in version of the group level image. We see that the peak location also occurred within the hand area. We note that this location coincided with one of the regions where $N = 9$ in Fig. 4A and C. The time-frequency vibrotactile data is shown in Fig. 4E, where we observed a similar behaviour in the beta band, whilst the alpha band showed a rebound starting later than for the beta band, which may be related to slight movement of the finger in response the vibrotactile stimulation (Caetano et al., 2007), a finding not observed with INMS. Analyses for the alpha band during the stimulation and resynchronisation epochs, and the beta in the resynchronisation epoch, can be found in the Supplementary material, but none of these other conditions showed a high likelihood of significance in individual subjects.

Fig. 5 shows the spatiotemporal similarities of responses to single-unit INMS and vibrotactile stimulation. Fig. 5A provides a conjunction map representing the number of unit locations for which there was a significant change in beta power in both the INMS and corresponding vibrotactile experiments. Here, we saw that for any given location, there was a maximum overlap in 8 subjects over the contralateral S1, whilst in the

<figure><img src="image_8.png" /><figcaption>**Fig. 4.** –Results of vibrotactile stimulation at the perceived receptive field of the 10 FA1 units studied overlaid onto inflated cortical surfaces, gyri are represented as light grey surfaces sulci are dark grey. A) Results of the individual vibrotactile analysis, represented as a conjunction map, showing how many exhibited significant changes in beta band (13–30 Hz) power for a given location. B) The group-averaged activation image of the 10 vibrotactile images, with a statistical ($p < 0.05$) threshold applied, showing vibrotaction elicited a reduction in beta oscillatory power over the contralateral sensorimotor area. The location of the maximal reduction of beta power is represented with the magenta dot. C) Zoomed portion of panel A, with the Brodmann areas as defined by the Glasser atlas (Glasser et al., 2016) overlaid in grey and hand area derived from a functional atlas (Sengupta et al., 2018) in black. Here the region representing the largest overlap across units (orange; $N = 9$) is contained primarily in the Brodmann areas 2 and 3b, both S1 areas. D) A zoomed in portion of panel B, again with Glasser atlas (green) and hand area (black) overlaid. The peak reduction of beta power (represented as the magenta dot) is located in postcentral gyrus. E) Time-frequency spectrogram of the peak location from panel B, showing a distinct event-related behaviour in the beta band which extends beyond the initial 1 s of stimulation. Interactive versions of the cortical plots can be found at http://georgeoneill.github.io.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Comparison of the spatiotemporal properties between INMS and the vibrotactile experiments. A) A conjunction map of where both INMS and vibrotactile experiments showed significant changes in beta power for a given mechanoreceptive afferent, and how many experiments spatially overlapped for a given area. B) A zoomed map of the right hemisphere's central gyri showing the largest group-averaged changes in power for INMS (magenta) and vibrotactile (cyan) and where they overlap with the locations of the largest spatial agreement as depicted in panel A (yellow). Any location where all three conditions held are represented in black. Brodmann areas, as defined using the Glasser atlas, and an fMRI derived hand area have been overlaid for reference. C) The group-averaged beta power timecourses for INMS (blue) and vibrotactile (red) responses, showing similar response timecourses in the contralateral (right) S1.</figcaption></figure>

ipsilateral hemisphere only activation in a single subject showed co-
significance. Further, at the group level, permutation testing found no
significant differences between ĀINMS and ĀVibrotactile in the beta band
during stimulation at any location within the brain. Fig. 5B shows the
conjunction map within the right S1 and M1. The map shows three
features; locations where ĀINMS was within 90% of the most extreme
value (magenta), ĀVibrotactile within 90% of the most extreme value
(cyan), and areas from Fig. 5A where n = 8 (yellow). Overlaying these
highlights two regions in S1 where all three conditions were satisfied
marked in black), one within Brodmann Area 3b, and one which
straddles the borders of Brodmann areas 2 and 1, the latter of which
agreed with the peak locations of ĀINMS and ĀVibrotactile. However, it is
important to note that due to the nature of unfolding source recon-
structed data onto a cortical surface, that these two distinct regions may
be potentially the same source projected onto both sides of the gyrus.
Fig. 5C shows the group-averaged oscillatory power in only the beta band
and compares the INMS beta timecourse to that of the vibrotactile ex-
periments in the same location. The shaded areas represent the standard
error across subjects. Non-parametric randomisation testing of the

timecourses revealed no significant differences at any time point between
the INMS and vibrotactile experiments.

4. Discussion

Here, we aimed to show the feasibility of combining single-unit INMS
(using recently-developed stimulation hardware; Glover et al., 2017)
with MEG, to probe whether we could reliably capture the spatial,
spectral, and temporal properties of the cortical projections of individual
single-unit mechanoreceptive afferents in human glabrous hand skin.

Our results show that single-unit INMS of these quantal elements of
touch yielded a robust reduction of oscillatory power in the beta
(13-30 Hz) band, which could be localised spatially using MEG, with
80% of individual FA1 units showing significant oscillatory power
reduction in contralateral S1. Furthermore, we compared our results with
oscillatory responses to amplitude-matched vibrotactile stimulation of
the same receptive field location and showed that our ability to localise
the ERD of individual mechanoreceptive afferent stimulation using INMS
was as reliable as localising sources derived from vibrotactile stimulation
and that the group results from single FA1 mechanoreceptive afferent stimulation revealed similar spatiotemporal patterns to those from vibrotactile stimuli. In addition, we also observed beta ERD to SA1 and FA2 units (see Supplementary material), showing that this was not a specific feature to FA1 units. The fundamental finding that single-unit INMS ERD responses were robust and consistent with natural somatosensory stimuli permits us to more dynamically probe central nervous system responses in humans. This will allow us to address a number of questions about the processing of touch from the different classes of mechanoreceptive afferent to a range of input patterns, including how varying the stimulus frequency and patterning in the different receptor types may modify cortical responses.

In the current study we chose to assess beta power, which was motivated by a number of studies observing modulations in beta power to tactile stimulation using both intracranial (Brovelli et al., 2004; Crone et al., 1998; Witham and Baker, 2007) and extracranial modalities (examples include, but are not limited to Bauer et al., 2006; Cheyne et al., 2003; Gaetz and Cheyne, 2006; Pfurtscheller et al., 2002; Salenius et al., 1997; van Ede et al., 2010). Bauer and colleagues used a spatial-selective tactile task to investigate how attention affected oscillatory power during tactile stimulation. They found the reduction of beta power during stimulation was not significantly affected by attentional modulation, suggesting this was likely an automatic response to afferent stimulation, independent of processing demand (Bauer et al., 2006). The same study also noted that the amount of post-event rebound in the somatosensory cortex was reduced if participants were instructed to attend to the stimulus, compared to those who were asked to ignore. The sensation from microsimulation is subtle and required subjects to attend to it, thus we expected a high level of attention. It is thus interesting to note that stimulation in only a few of the units produced significant rebounds of beta power, and these all originated from the same participant, a participant well accustomed to single-unit INMS and so more familiar to the sensation, thus potentially requiring less attentional load. Another study has shown that anticipation may modulate event-related reductions in beta, in research employing a delay between a cue and the delivery of a stimulus, where there was contralateral ERD of the beta band, with the extent of the reduction monotonically increasing with delay time, implying an anticipatory element of somatosensory beta (van Ede et al., 2011). However, in the present study we used no explicit cues and therefore have confidence that our observations in beta power reduction are entirely due to stimulation.

Recently, there has been an emerging shift in how we interpret beta band oscillations in the somatosensory and motor cortices (van Ede et al., 2018). Rather than treating beta as a sustained oscillation, which is power modulated based on processing demand, it may be better to discuss beta in terms of transient (~50 ms) events, or bursts. These punctate events have been shown to be robust across multiple species (Sherman et al., 2016; Shin et al., 2017), and their occurrence in human experiments seems to imply these are inhibitory processes (Sacchet et al., 2015a; Shin et al., 2017). Furthermore, the probability of the occurrence of a burst follows the trial-averaged beta power in the sensorimotor cortex; a rarefication of bursting events occurs during the ERD period and maximal likelihood coincides with the post movement beta rebound (Little et al., 2018). Furthermore, there are some early indications that these beta events form long-range transient synchronisations across the brain (Sacchet et al., 2015b), which could make bursts a potential mechanism of transient functional networks in the brain (O'Neill et al., 2018). There is also evidence that this bursting behaviour may facilitate communication to other regions of the body, as bursts centred around 20 Hz have been shown to be coherent with EMG recordings in the hand during a grip force task (Bourguignon et al., 2017). To be able to fully understand the mechanisms of these bursts, we require experimental validations for predictions which mathematical models of this bursting mechanism (Jones et al., 2009) make. Given the precise nature of single-unit INMS, it is a suitable experimental candidate to explore such mechanisms.

In a previous study by Kelly et al. (1997), INMS was combined with EEG, where they identified sustained evoked potentials that followed modulations in stimulation frequency. They found that both single-unit INMS and vibrotactile stimulation could evoke clear frequency following responses (FFRs) and these were matched to the stimulation frequency (33 Hz; Kelly et al., 1997). In our current study, we did not observe any such FFRs for neither INMS nor vibrotactile stimulation, although we used a stimulus frequency of 60 Hz. We can identify two potential reasons why this was the case, which are both linked to the experimental design. The first can be attributed to the number of trials per experiment. In the Kelly et al. study they used 200 trials of 1 s stimulation, whereas here we used 80 trials. Our choice of 80 was to optimise the experiment for capturing the modulation in beta power, which in the sensorimotor system takes several seconds to return to baseline after an event (Barratt et al., 2017; Fry et al., 2016; Robson et al., 2016), and so to gain an adequate inter-trial interval, the total number of trials is fewer. Assuming zero trials rejected, the Kelly et al. study has a 58% signal to noise ratio (SNR) improvement over our paradigm and this may be crucial to see such effects. The second reason, which may be more important than number of trials, is the stimulation frequency. Our choice of 60 Hz stimulation was driven by this providing a clearly perceivable, constant sensation for subjects when stimulating both FA1 and SA1 afferents (Johnson, 2001). However, it is likely that this frequency of stimulation is too high to generate a (detectible) FFR within the somatosensory cortex with MEG. Previous literature on the non-invasive imaging mechanoreceptive FFRs (see Vlaar et al., 2015 for a brief review) implies that the bulk of studies have been performed at a stimulation frequency of less than 40 Hz, and those which use higher carrier frequencies (such as Tobi-matsu et al., 1999) used amplitude modulation below 40 Hz. Studies characterising the optimal frequencies show that the optimal frequency for an FFR is around 20–26 Hz (Jamali and Ross, 2013; Snyder, 1992; Tobimatsu et al., 1999). There has been a single report into the detection of subharmonics (Langdon et al., 2011), but their carrier frequencies were in this optimal spectral window.

In the present study, we did not perform a systematic analysis of the early somatosensory evoked fields since the experimental paradigm was designed to assess oscillatory responses. Ideally, single pulses would be used in an evoked fields analysis. In single pulse INMS, it is difficult for the subject to focus on such a brief sensation (i.e. to detect a single pulse coming from a single afferent), and such sensations have only been reported when stimulating FA1 units (Vallbo et al., 1984b). In this first study using MEG to study cortical responses to INMS, we wanted to ensure that we could gain a perceptual response to test, thus we delivered a longer train of pulses known to elicit perceivable sensations and detectable brain responses in EEG (Kelly et al., 1997). The main disadvantage of using such a train of pulses is that evoked fields from consecutive stimulation pulses will superimpose. Given that in our study we only had ~16.6 ms between pulses, it would be difficult to deconvolve and resolve specific waveforms (e.g. M20, M50). To specifically investigate evoked fields using trains of stimuli, ~500 ms intervals between pulses is required and will be of interest to address in future work. Nevertheless, we acknowledge whether single unit INMS can generate detectable evoked fields will be an interest to many, and have therefore performed a basic analysis of early evoked responses (see Appendix). Even with our suboptimal paradigm, we were able to localise the strongest source activations in the first 100 ms to S1 and S2. Therefore, this provides the opportunity of investigating whether differences between classes of mechanoreceptors are encoded in the earliest neural responses, or in later and slower oscillatory power modulations.

In our current study, we used bespoke flexible headcasts to improve the reliability of the spatial localisation of the oscillatory sources and comparison of INMS and vibrotactile responses between MEG scan sessions. The use of such headcasts has been shown to minimise head movement and the inter-session errors which arise from sensor-to-
anatomy co-registration (Liuzzi et al., 2017). Since it was unknown how spatiotemporal responses to single-unit INMS would manifest in the human cortex, we wanted to ensure that we maximized the robustness of localisation of the sources between INMS and vibrotactile sessions using optimal methods. However, given that one unit (Unit 2a) was captured without a headcast and was localised to somatosensory cortex (see Supplementary material), with a source amplitude similar to that of the vibrotactile stimulation, future studies may not necessarily require a headcast, as long as sufficient practices in reducing head motion and co-registration errors are adhered to. Further the use of the headcast required care when inserting the participant into the MEG sensor array once a unit was found, which has the potential of causing small electrode movements, which may lead to the loss of the single-unit nature of the stimulation. We originally intended to perform the search for mechanoreceptive afferents whilst the participant was within the MEG sensor array, however due to an unacceptable high level of capacitively-coupled mains-borne interference related to the MEG system earth in our installation, this interference obscured the neural responses, making effective identification of single-unit recordings near impossible whilst the head was inside the MEG sensor array. The reference and recording electrodes used for INMS have very different impedances which makes the Common Mode Rejection Ratio performance for the pre-amplifier quite poor, but usually in screened room environments this is not an issue. We ameliorated this problem by lowering the participant out of the MEG sensor array and cutting the power to the MEG gantry during the search, with

the former of these solutions necessitating careful entry and exit into the sensor array before and after data collection. It should be noted that cryogen-free on-scalp MEG solutions are starting to come into use (Boto et al., 2018, 2017), and since they provide of the order of 5-8 fold improvement in SNR over the central gyri and have sensors which move with the head, these will potentially provide flexibility for combined single-unit INMS without compromising on spatial specificity (Boto et al., 2016; liivanainen et al., 2017).

In conclusion, combining the spatiotemporal specificity of MEG with the stimulation specificity of single-unit INMS enables the interrogation of the central representations of different aspects of tactile afferent signalling within the human cortices. Our initial findings imply that beta power in the somatosensory cortex is significantly modulated from the stimulation of single mechanoreceptive afferents, and given the precise nature of the stimulation method, single-unit INMS may serve as a controlled means to explore the mechanisms behind beta bursting events.

## Acknowledgements

This work was supported by the Medical Research Council [grant number MR/M022722/1, MR/K005464/1], to STF and MEG UK partnership respectively; the Swedish Research Council [grant 2017-01717], to JW; and the Knut and Alice Wallenberg Foundation [project Neuro-SQUID], to JW.

# Appendix: A. preliminary analysis into the early onset responses of INMS.

As stated within the main manuscript, our experimental design means that our study is not particularly optimised for characterising the early onset somatosensory evoked fields; with only ~16.6 ms between stimulation pulses, it would be difficult to deconvolve the effects of the adjacent stimulations. This will be addressed in future study, but for the time being we present a preliminary investigation of the early onset responses from our 10 FA1 units.

## Methods

In this analysis we present the evoked response to the entire train of stimulation pulses, rather than individual stimulations. Our methods to analyse the data are in places different to the analysis in the main manuscript and so a description follows. After the removal of artefacts from the data, we filtered between 1 and 80 Hz. For source reconstruction, we used the same source locations on the cortical surface and same forward model as before (3-shell BEM), but instead opted to use minimum norm estimation rather than a beamformer for inverse modelling as they are better suited to localising evoked sources (Ou et al., 2009). In particular, we used the sLORETA variant (Pascual-Marqui, 2002), where the sensor level noise covariance matric was generated from data 2 s prior to the stimulation onset. For a given vertex location $r$, three sets of source reconstruction weights representing three orthogonal dipoles were generated, which we shall define as $w_{rx}$, $w_{ry}$ and $w_{rz}$ respectively. The sLORETA score is then calculated for given time point $t$ by

$$s(\mathbf{r}, t) = \sqrt{(w_{rx}b(t))^2 + (w_{ry}b(t))^2 + (w_{rz}b(t))^2}, \quad (A1)$$

where $b$ is the MEG sensor level data.

## Results

The results of this analysis can be found in Figure A1. Before focusing on the source level results, Figure A1A represents grand average sensor level plot, across all trials and subjects. We see there is a clear modulation of power oscillatory power within the 0-1 s epoch of microsimulation, both reflected in the higher amplitudes in the butterfly plot (coloured traces) and the increase in global field power (black trace) before eventually recovering back to baseline levels. Figure A1B is the source localised image of the earliest evoked responses in the contralateral hemisphere, averaged across both subjects and the first 100 ms of stimulation, with a threshold of 80% of the maximum sLORETA value applied. We observed two anatomical regions where activation is the strongest, the postcentral gyrus (S1) and the supramarginal gyrus (S2). Overlaid on top of the source plot in black is the fMRI derived hand region, which shows many of the peaks from the source plot contained within. The finding that the hand area intersected with this early source activation map and our areas of maximal beta power reduction (in Figs. 3 and 4) increasing our confidence that those results were due to the stimulation and not other errant effects. Figure A1C and A1D show the average sLORETA timecourses in S1 and S2 in both hemispheres. We found that there were two distinct events in the data, a first early component that peaked in intensity at ~200 ms, and a broader latent component towards the end of the stimulation. Note that the first component is also present in the main manuscript, represented as an increase in oscillatory power below 10 Hz (see Fig. 3E). We also saw a bilateral response, with diminished responses in the ipsilateral (left) S1 and S2, which also appear to be lagged in time. This bilateral S2 response has been observed in somatosensory stimulation before (Gao et al., 2015). The latent effect of the ipsilateral hemisphere can be
seen clearer in Figure A1D. From these results we are confident that with careful experimental design, evoked components related to single unit stimulation have the potential to be detected and characterised.

<figure><img src="image_10.png" /><figcaption>**Fig. A1.** Results of a preliminary analysis into the early onset responses from INMS. A) A sensor-level grand average of the 10 FA1 units. B) A grand average sensor level plot of the first 100 ms of stimulation, showing the strongest sources in contralateral S1 (postcentral gyrus) and S2 (supramarginal) regions. For reference, the black overlay is the fMRI derived hand area (Sengupta et al., 2018). C) The average sLORETA timecourses from the left and right S1 and S2 areas. D) The same timecourses zoomed in to the first 600 ms, which shows the delayed ipsilateral response from both S1 and S2. An interactive version of the cortical plot can be found at http://georgeoneill.github.io.</figcaption></figure>

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.neuroimage.2019.01.017.

## References

Ackerley, R., Kavounoudias, A., 2015. The role of tactile afference in shaping motor behaviour and implications for prosthetic innovation. Neuropsychologia 79, 192-205. https://doi.org/10.1016/j.neuropsychologia.2015.06.024.

Barratt, E.L., Tewarie, P.K., Clarke, M.A., Hall, E.L., Gowland, P.A., Morris, P.G., Francis, S.T., Evangelou, N., Brookes, M.J., 2017. Abnormal task driven neural oscillations in multiple sclerosis: A visuomotor MEG study. Hum. Brain Mapp. 38, 2441-2453. https://doi.org/10.1002/hbm.23531.

Bauer, M., Oostenveld, R., Peeters, M., Fries, P., 2006. Tactile Spatial Attention Enhances Gamma-Band Activity in Somatosensory Cortex and Reduces Low-Frequency Activity in Parieto-Occipital Areas. J. Neurosci. 26, 490-501. https://doi.org/10.1523/JNEUROSCI.5228-04.2006.

Boto, E., Bowtell, R., Krüger, P., Fromhold, T.M., Morris, P.G., Meyer, S.S., Barnes, G.R., Brookes, M.J., 2016. On the potential of a new generation of magnetometers for MEG: A beamformer simulation study. PLoS One 11. https://doi.org/10.1371/journal.pone.0157655.

Boto, E., Holmes, N., Leggett, J., Roberts, G., Shah, V., Meyer, S.S., Muñoz, L.D., Mullinger, K.J., Tierney, T.M., Bestmann, S., Barnes, G.R., Bowtell, R., Brookes, M.J., 2018. Moving magnetoencephalography towards real-world applications with a wearable system. Nature 555, 657-661. https://doi.org/10.1038/nature26147.

Boto, E., Meyer, S.S., Shah, V., Alem, O., Knappe, S., Kruger, P., Fromhold, T.M., Lim, M., Glover, P.M., Morris, P.G., Bowtell, R., Barnes, G.R., Brookes, M.J., 2017. A new generation of magnetoencephalography: Room temperature measurements using optically-pumped magnetometers. Neuroimage 149, 404-414. https://doi.org/10.1016/j.neuroimage.2017.01.034.

Bourguignon, M., Piitulainen, H., Smeds, E., Zhou, G., Jousmäki, V., Hari, R., 2017. MEG insight into the spectral dynamics underlying steady isometric muscle contraction. J. Neurosci. https://doi.org/10.1523/JNEUROSCI.0447-17.2017.

Brookes, M.J., Vrba, J., Robinson, S.E., Stevenson, C.M., Peters, A.M., Barnes, G.R., Hillebrand, A., Morris, P.G., 2008. Optimising experimental design for MEG beamformer imaging. Neuroimage 39, 1788-1802.

Brovelli, A., Ding, M., Ledberg, A., Chen, Y., Nakamura, R., Bressler, S.L., 2004. Beta oscillations in a large-scale sensorimotor cortical network: Directional influences revealed by Granger causality. Proc. Natl. Acad. Sci. Unit. States Am. 101, 9849-9854. https://doi.org/10.1073/pnas.0308538101.

Caetano, G., Jousmaki, V., Hari, R., 2007. Actor's and observer's primary motor cortices stabilize similarly after seen or heard motor actions. Proc. Natl. Acad. Sci. Unit. States Am. https://doi.org/10.1073/pnas.0702453104.

Cheyne, D., Gaetz, W., Garnero, L., Lachaux, J.P., Ducorps, A., Schwartz, D., Varela, F.J., 2003. Neuromagnetic imaging of cortical oscillations accompanying tactile stimulation. Cognit. Brain Res. 17, 599-611. https://doi.org/10.1016/S0926-6410(03)00173-3.

Cheyne, D.O., 2013. MEG studies of sensorimotor rhythms: A review. Exp. Neurol. https://doi.org/10.1016/j.expneurol.2012.08.030.

Crone, N.E., Miglioretti, D.L., Gordon, B., Sieracki, J.M., Wilson, M.T., Uematsu, S., Lesser, R.P., 1998. Functional mapping of human sensorimotor cortex with electrocorticographic spectral analysis. I. Alpha and beta event-related desynchronization. Brain 121, 2271.

Dale, A.M., Fischl, B., Sereno, M.I., 1999. Cortical surface-based analysis: I. Segmentation and surface reconstruction. Neuroimage 9, 179-194. https://doi.org/10.1006/nimg.1998.0395.
Fischl, B., Sereno, M.I., Dale, A.M., 1999. Cortical surface-based analysis: II. Inflation, flattening, and a surface-based coordinate system. Neuroimage. https://doi.org/10.1006/nimg.1998.0396.

Fry, A., Mullinger, K.J., O'Neill, G.C., Barratt, E.L., Morris, P.G., Bauer, M., Folland, J.P., Brookes, M.J., 2016. Modulation of post-movement beta rebound by contraction force and rate of force development. Hum. Brain Mapp. 37, 2493-2511. https://doi.org/10.1002/hbm.23189.

Gaetz, W., Cheyne, D., 2006. Localization of sensorimotor cortical rhythms induced by tactile stimulation using spatially filtered MEG. Neuroimage 30, 899-908. https://doi.org/10.1016/j.neuroimage.2005.10.009.

Gao, L., Sommerlade, L., Coffman, B., Zhang, T., Stephen, J.M., Li, D., Wang, J., Grebogi, C., Schelter, B., 2015. Granger causal time-dependent source connectivity in the somatosensory network. Sci. Rep. https://doi.org/10.1038/srep10399.

Glasser, M.F., Coalson, T.S., Robinson, E.C., Hacker, C.D., Harwell, J., Yacoub, E., Ugurbil, K., Andersson, J., Beckmann, C.F., Jenkinson, M., Smith, S.M., Van Essen, D.C., 2016. A multi-modal parcellation of human cerebral cortex. Nature 536, 171-178. https://doi.org/10.1038/nature18933.

Glover, P.M., Watkins, R.H., O'Neill, G.C., Ackerley, R., Sanchez-Panchuelo, R., McGlone, F., Brookes, M.J., Wessberg, J., Francis, S.T., 2017. An intra-neural microstimulation system for ultra-high field magnetic resonance imaging and magnetoencephalography. J. Neurosci. Methods 290, 69-78. https://doi.org/10.1016/j.jneumeth.2017.07.016.

Hämäläinen, M., Hari, R., Ilmoniemi, R.J., Knuutila, J., Lounasmaa, O.V., 1993. Magnetoencephalography theory, instrumentation, and applications to noninvasive studies of the working human brain. Rev. Mod. Phys. 65, 413-497. https://doi.org/10.1103/RevModPhys.65.413.

Hari, R., Salmelin, R., 1997. Human cortical oscillation: a magnetic visual view through the skull. Trends Neurosci. 20, 44-49. https://doi.org/10.1016/S0166-2236(96)10065-5.

Iivanainen, J., Stenroos, M., Parkkonen, L., 2017. Measuring MEG closer to the brain: Performance of on-scalp sensor arrays. Neuroimage 147, 542-553. https://doi.org/10.1016/j.neuroimage.2016.12.048.

Jamali, S., Ross, B., 2013. Somatotopic finger mapping using MEG: Toward an optimal stimulation paradigm. Clin. Neurophysiol. 124, 1659-1670. https://doi.org/10.1016/j.clinph.2013.01.027.

Johansson, R.S., Vallbo, A.B., 1979. Tactile sensibility in the human hand: relative and absolute densities of four types of mechanoreceptive units in glabrous skin. J. Physiol. 286, 283-300. https://doi.org/10.1113/jphysiol.1979.sp012619.

Johnson, K.O., 2001. The roles and functions of cutaneous mechanoreceptors. Curr. Opin. Neurobiol. https://doi.org/10.1016/S0959-4388(00)00234-8.

Jones, S.R., Pritchett, D.L., Sikora, M.A., Stufflebeam, S.M., Hamalainen, M., Moore, C.I., 2009. Quantitative Analysis and Biophysically Realistic Neural Modeling of the MEG Mu Rhythm: Rhythmogenesis and Modulation of Sensory-Evoked Responses. J. Neurophysiol. 102, 3554-3572. https://doi.org/10.1152/jn.00535.2009.

Kelly, E.F., Trulsson, M., Folger, S.E., 1997. Periodic microstimulation of single mechanoreceptive afferents produces frequency-following responses in human EEG. J. Neurophysiol. 77, 137-144.

Langdon, A.J., Boonstra, T.W., Breakspear, M., 2011. Multi-frequency phase locking in human somatosensory cortex. Prog. Biophys. Mol. Biol. 105, 58-66. https://doi.org/10.1016/j.pbiomolbio.2010.09.015.

Little, S., Bonaiuto, J., Barnes, G., Bestmann, S., 2018. Motor cortical beta transients delay movement initiation and track errors. bioRxiv. https://doi.org/10.1101/384370.

Liuzzi, L., Gascoyne, L.E., Tewarie, P.K., Barratt, E.L., Boto, G.R., E. Brookes, M.J., 2017. Optimising experimental design for MEG resting state functional connectivity measurement. Neuroimage 155, 565-576. https://doi.org/10.1016/j.neuroimage.2016.11.064.

Maris, E., Oostenveld, R., 2007. Nonparametric statistical testing of EEG- and MEG-data. J. Neurosci. Methods 164, 177-190. https://doi.org/10.1016/j.jneumeth.2007.03.024.

McGlone, F., Wessberg, J., Olausson, H., 2014. Discriminative and Affective Touch: Sensing and Feeling. Neuron. https://doi.org/10.1016/j.neuron.2014.05.001.

Meyer, S.S., Bonaiuto, J., Lim, M., Rossiter, H., Waters, S., Bradbury, D., Bestmann, S., Brookes, M., Callaghan, M.F., Weiskopf, N., Barnes, G.R., 2017. Flexible head-casts for high spatial precision (MEG). J. Neurosci. Methods 276, 38-45. https://doi.org/10.1016/j.jneumeth.2016.11.009.

Nichols, T.E., Holmes, A.P., 2002. Nonparametric permutation tests for functional neuroimaging: A primer with examples. Hum. Brain Mapp. 15, 1-25. https://doi.org/10.1002/hbm.1058.

O'Neill, G.C., Tewarie, P., Vidaurre, D., Liuzzi, L., Woolrich, M.W., Brookes, M.J., 2018. Dynamics of large-scale electrophysiological networks: A technical review. Neuroimage 180. https://doi.org/10.1016/j.neuroimage.2017.10.003.

Ochoa, J., Torebjörk, E., 1983. Sensations evoked by intraneural microstimulation of single mechanoreceptor units innervating the human hand. J. Physiol. 342, 633-654. https://doi.org/10.1113/jphysiol.1983.sp014873.

Oostenveld, R., Fries, P., Maris, E., Schoffelen, J.-M., 2011. FieldTrip: open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data. Comput. Intell. Neurosci. 1, 2011.

Ou, W., Hämäläinen, M.S., Golland, P., 2009. A distributed spatio-temporal EEG/MEG inverse solver. Neuroimage 44, 932-946. https://doi.org/10.1016/J.NEIROIMAGE.2008.05.063.

Pantazis, D., Nichols, T.E., Baillet, S., Leahy, R.M., 2005. A comparison of random field theory and permutation methods for the statistical analysis of MEG data. Neuroimage 25, 383-394. https://doi.org/10.1016/j.neuroimage.2004.09.040.

Pascual-Marqui, R.D., 2002. Standardized low-resolution brain electromagnetic tomography (sLORETA): technical details. Methods Find. Exp. Clin. Pharmacol. (lpiii) https://doi.org/841.

Pfurtscheller, G., Woertz, M., Müller, G., Wriessnegger, S., Pfurtscheller, K., 2002. Contrasting behavior of beta event-related synchronization and somatosensory evoked potential after median nerve stimulation during finger manipulation in man. Neurosci. Lett. 323, 113-116. https://doi.org/10.1016/S0304-3940(02)00119-2.

Robson, S.E., Brookes, M.J., Hall, E.L., Palaniyappan, L., Kumar, J., Skelton, M., Christodoulou, N.G., Qureshi, A., Jan, F., Katshu, M.Z., Liddle, E.B., Liddle, P.F., Morris, P.G., 2016. Abnormal visuomotor processing in schizophrenia. NeuroImage Clin 12, 869-878. https://doi.org/10.1016/j.nicl.2015.08.005.

Saal, H.P., Bemsaia, S.J., 2014. Touch is a team effort: Interplay of submodalities in cutaneous sensibility. Trends Neurosci. https://doi.org/10.1016/j.tins.2014.08.012.

Sacchet, M.D., LaPlante, R.A., Wan, Q., Pritchett, D.L., Lee, A.K.C., Hamalainen, M., Moore, C.I., Kerr, C.E., Jones, S.R., 2015a. Attention Drives Synchronization of Alpha and Beta Rhythms between Right Inferior Frontal and Primary Sensory Neocortex. J. Neurosci. 35, 2074-2082. https://doi.org/10.1523/JNEUROSCI.1292-14.2015.

Sacchet, M.D., LaPlante, R.A., Wan, Q., Pritchett, D.L., Lee, A.K.C., Hamalainen, M., Moore, C.I., Kerr, C.E., Jones, S.R., 2015b. Attention Drives Synchronization of Alpha and Beta Rhythms between Right Inferior Frontal and Primary Sensory Neocortex. J. Neurosci. 35, 2074-2082. https://doi.org/10.1523/JNEUROSCI.1292-14.2015.

Senalius, S., Schnitzler, A., Salmelin, R., Jousmäki, V., Hari, R., 1997. Modulation of human cortical rolandic rhythms during natural sensorimotor tasks. Neuroimage 5, 221-228. https://doi.org/10.1006/nimg.1997.0261.

Sanchez Panchuelo, R.M., Ackerley, R., Glover, P.M., Bowtell, R.W., Wessberg, J., Francis, S.T., McGlone, F., 2016. Mapping quantal touch using 7 Tesla functional magnetic resonance imaging and single-unit intraneural microstimulation. Elife 5, e12812.

Sekihara, K., Nagarajan, S.S., Poeppel, D., Marantz, A., 2004. Asymptotic SNR of scalar and vector minimum-variance beamformers for neuromagnetic source reconstruction. IEEE Trans. Biomed. Eng. 51, 1726-1734. https://doi.org/10.1109/TBME.2004.827926.

Sengupta, A., Schluppeck, D., Barratt, E.L., Besle, J., Francis, S.T., Sanchez-Panchuelo, R., 2018. A Probabilistic Atlas of Digit Somatotopy in the Human Primary Somatosensory Cortex. In: Proc. Intl. Soc. Mag. Reson. Med. Paris, France, p. 5412.

Sherman, M.A., Lee, S., Law, R., Haegens, S., Thorn, C.A., Hämäläinen, M.S., Moore, C.I., Jones, S.R., 2016. Neural mechanisms of transient neocortical beta rhythms: Converging evidence from humans, computational modeling, monkeys, and mice. Proc. Natl. Acad. Sci. Unit. States Am. 113, E4885-E4894. https://doi.org/10.1073/pnas.1604135113.

Shin, H., Law, R., Tsutsui, S., Moore, C.I., Jones, S.R., 2017. The rate of transient beta frequency events predicts behavior across tasks and species. Elife 6. https://doi.org/10.7554/eLife.29086.

Singh, K.D., Barnes, G.R., Hillebrand, A., 2003. Group imaging of task-related changes in cortical synchronisation using nonparametric permutation testing. Neuroimage 19, 1589-1601. https://doi.org/10.1016/S1053-8119(03)00249-0.

Smith, S.M., Nichols, T.E., 2009. Threshold-free cluster enhancement: Addressing problems of smoothing, threshold dependence and localisation in cluster inference. Neuroimage 44, 83-98. https://doi.org/10.1016/j.neuroimage.2008.03.061.

Snyder, A.Z., 1992. Steady-state vibration evoked potentials: description of technique and characterization of responses. Electroencephalogram. Clin. Neurophysiol. Evoked Potentials 84, 257-268. https://doi.org/10.1016/0168-5597(92)90007-X.

Stenroos, M., Nummenmaa, A., 2016. Incorporating and Compensating Cerebrospinal Fluid in Surface-Based Forward Models of Magneto- and Electroencephalography. PLoS One 11, e0159595. https://doi.org/10.1371/journal.pone.0159595.

Tobimatsu, S., Zhang, Y.M., Kato, M., 1999. Steady-state vibration somatosensory evoked potentials: Physiological characteristics and tuning function. Clin. Neurophysiol. 110, 1953-1958. https://doi.org/10.1016/S1388-2457(99)00146-7.

Torebjörk, H.E., Ochoa, J.L., 1980. Specific sensations evoked by activity in single identified sensory units in man. Acta Physiol. Scand. 110, 445-447. https://doi.org/10.1111/j.1748-1716.1980.tb06695.x.

Torebjörk, H.E., Vallbo, Å.B., Ochoa, J.L., 1987. Intraneural microstimulation in man: Its relation to specificity of tactile sensations. Brain 110, 1509-1529.

Troeiböker, L., López, J.D., Lutti, A., Bradbury, D., Bestmann, S., Barnes, G., 2014. High precision anatomy for MEG. Neuroimage 86, 583-591. https://doi.org/10.1016/j.neuroimage.2013.07.065.

Trulsson, M., Francis, S.T., Kelly, E.F., Westling, G., Bowtell, R., McGlone, F., 2001. Cortical responses to single mechanoreceptive afferent microstimulation revealed with fMRI. Neuroimage 13, 613-622.

Vallbo, Å.B., 1981. Sensations evoked from the glabrous skin of the human hand by electrical stimulation of unitary mechanosensitive afferents. Brain Res. 215, 359-363. https://doi.org/10.1016/0006-8993(81)90517-5.

Vallbo, Å.B., Hagbarth, K.-E., 1968. Activity from skin mechanoreceptors recorded percutaneously in awake human subjects. Exp. Neurol. 21, 270-289.

Vallbo, Å.B., Johansson, R.S., others, 1984a. Properties of cutaneous mechanoreceptors in the human hand related to touch sensation. Hum. Neurobiol. 3, 3-14.

Vallbo, Å.B., Olsson, K.A., Westberg, K.G., Clark, F.J., 1984b. Microstimulation of single tactile afferents from the human hand. Sensory attributes related to unit type and properties of receptive fields. Brain 107 (Pt 3), 727-749. https://doi.org/10.1093/brain/107.3.727.

van Ede, F., de Lange, F., Jensen, O., Maris, E., 2011. Orienting Attention to an Upcoming Tactile Event Involves a Spatially and Temporally Specific Modulation of Sensorimotor Alpha- and Beta-Band Oscillations. J. Neurosci. 31, 2016-2024. https://doi.org/10.1523/JNEUROSCI.5630-10.2011.
van Ede, F., Jensen, O., Maris, E., 2010. Tactile expectation modulates pre-stimulus A-band oscillations in human sensorimotor cortex. Neuroimage 51, 867–876. https://doi.org/10.1016/j.neuroimage.2010.02.053.

van Ede, F., Quinn, A.J., Woolrich, M.W., Nobre, A.C., 2018. Neural Oscillations: Sustained Rhythms or Transient Burst-Events? Trends Neurosci. (in press).

Vlaar, M.P., van der Helm, F.C.T., Schouten, A.C., 2015. Frequency Domain Characterization of the Somatosensory Steady State Response in

Electroencephalography. IFAC-PapersOnLine 48, 1391–1396. https://doi.org/10.1016/j.ifacol.2015.12.327.

Witham, C.L., Baker, S.N., 2007. Network oscillations and intrinsic spiking rhythmicity do not covary in monkey sensorimotor areas. J. Physiol. 580, 801–814. https://doi.org/10.1113/jphysiol.2006.124503.