Multisine frequency modulation of intra-epidermal electric pulse sequences: A novel tool to study nociceptive processing

A sustained sensory stimulus with a periodic variation of intensity creates an electrophysiological brain response at associated frequencies, referred to as the steady-state evoked potential (SSEP). The SSEPs elicited by the periodic stimulation of nociceptors in the skin may represent activity of a brain network that is primarily involved in nociceptive processing. Exploring the behavior of this network could lead to valuable insights regarding the pathway from nociceptive stimulus to pain perception. We present a method to directly modulate the pulse rate of nociceptive afferents in the skin with a multisine waveform through intra-epidermal electric stimulation. The technique was demonstrated in healthy volunteers. Each subject was stimulated using a pulse sequence modulated by a multisine waveform of 3, 7 and 13 Hz. The EEG was analyzed for the presence of the base frequencies and associated (sub)harmonics. Topographies showed significant central and contralateral SSEP responses at 3, 7 and 13 Hz in respectively 7, 4 and 3 out of the 9 participants included for analysis. As such, we found that intra-epidermal stimulation with a multisine frequency modulated pulse sequence can generate nociceptive SSEPs. The possibility to stimulate the nociceptive system using multisine frequency modulated pulses offers novel opportunities to study the temporal dynamics of nociceptive processing.


Introduction
Despite decades of research, it remains unclear how our brain creates the perception of pain based on a combination of peripheral nociceptive input and central nervous activity. To study the relation between peripheral nociceptive input and brain activity, researchers have extensively documented the brain potential evoked by a single painful stimulus. These pain evoked potentials provided one of the first tools to evaluate properties of nociceptive processing in both healthy and abnormal conditions. As such, pain evoked potentials have been used in a wide variety of contexts to study pain processing and modulation (e.g. (Liang et al., 2016;Manresa et al., 2018;Wager et al., 2006)). These pain evoked potentials were shown to be associated with stimulus saliency  rather than any specific effects of stimulus modality (e.g. nociceptive or somatosensory).
An alternative technique to characterize stimulus-evoked brain activity is the measurement of steady-state evoked potentials (SSEPs). Instead of evoking brain activity using a single transient stimulus, SSEPs rely on sustained activation of brain areas by a continuous sensory stimulus. The tonic pain evoked by the continuous nociceptive stimulus leading to a nociceptive SSEP bears more similarity to the continuous and dynamic pain experienced by patients in daily life than the transient stimuli used to evoke brain potentials, and reduces the effect of saliency. By modulating the intensity of this stimulus with one or multiple frequencies, evoked brain activity can be observed in the electroencephalogram (EEG) at these frequencies and their harmonics. Visual and auditory SSEPs have been widely used in cognitive and clinical neuroscience (Vialatte et al., 2010;Norcia et al., 2015). Sensory transmission and processing can be further characterized by evoking SSEPs using multisine waveforms, perturbing the sensory system at multiple frequencies simultaneously. Some of the useful properties that could be quantified using this technique include the delay, nonlinearity and signal-to-noise ratio of a sensory system (Yang et al., 2016). Recent studies have used the technique to observe fundamental properties of the sensorimotor system (Yang et al., 2017), and to demonstrate significantly altered transmission of sensory signals after stroke (Vlaar et al., 2017). Another important topic of current SSEP research is the (non)linearity of sensory processing. Modelling studies suggest that brain responses are highly nonlinear (Roberts and Robinson, 2012;Spiegler et al., 2011) and investigating these nonlinear relations between sensory input and brain activity is essential for our understanding of sensory systems. Such nonlinear relations have been shown clinically relevant for motor disorders (Sanger et al., 2002), migraine (Nyrke and Lang, 1982) and epilepsy (Kalitzin et al., 2002). The (non)linearity of nociceptive processing and its potential clinical applications remains relatively unexplored, which is an area where multisine SSEP techniques might prove valuable as a relatively cheap and non-invasive technique to explore this topic.
Recently, Mouraux et al. (2011) showed that it was possible to measure SSEPs related to specific nociceptive stimulation by applying a periodic sequence of laser pulses. This lab showed that the technique, also referred to as 'frequency tagging', could be used to dissociate cortical responses to nociceptive and tactile stimuli (Colon et al., 2014). More recently, the same lab has also demonstrated SSEPs generated by slow periodic variation of the heat in a contact thermode (Colon et al., 2017;Mulders et al., 2020), eliciting SSEPs with a similar topography, but limited to ultra-low stimulation frequencies. Nevertheless, a limitation of both thermode stimulation and laser stimulation is that nociceptive nerve fibers are activated through a nonlinear thermal transduction process (Xu et al., 2010), which could limit the observability of central nociceptive processing.
The thermal transduction process can be bypassed by intraepidermal electric stimulation, which has also been shown to preferentially activate nociceptive afferents, provided that stimulus intensity remains at or below twice the detection threshold . However, the benefit of bypassing transduction processes comes with an additional risk of eliciting stimulus artefacts in the EEG. A first study used square-wave modulated pulse sequences to elicit SSEPs through intra-epidermal stimulation (Colon et al., 2012a) and successfully demonstrated the potential of intra-epidermal electric stimulation by measuring SSEPs with a large signal-to-noise ratio on a range of frequencies from 3 to 43 Hz. Another study successfully evoked SSEPs at 31 and 37 Hz by stimulation of superficial skin afferents using a concentric planar electrode (Blöchl et al., 2015). However, both studies showed distinctly different SSEP topographies and lacked (description of) rigorous checks for stimulus artefacts and a description of the variability of elicited SSEPs among participants, making it difficult to anticipate potential applications. Furthermore, the square wave pulse train modulation used in these studies cannot be adapted for perturbing multiple frequencies simultaneously and already includes spectral peaks at harmonics in the input signal, which compromises usage of most system identification methods.
In this study, we developed a new procedure for eliciting multifrequency nociceptive SSEPs through intra-epidermal stimulation. Accurate frequency modulated stimulation of nociceptive nerve fibers and measurement of the generated SSEPs creates many challenges with respect to hardware, stimulation procedures, recording procedures and analysis. As such, custom hardware was developed to modulate electric pulse sequences with sinusoid and multisine waveforms. We evaluated the performance and the limits of this method, and outlined a procedure to accurately quantify and map stimulus artefacts. Subsequently, the effectiveness of multisine frequency modulation of a pulse sequence for generating and studying nociceptive SSEPs was evaluated in an experiment on ten participants. We tested our hypothesis that the technique would generate significant peaks in the power spectrum at the base frequencies and potentially at some of the harmonics. Furthermore, we used the evoked spectral components to estimate system delay and explore system nonlinearity.

Participants
A group of ten healthy men (aged 23-27 years, nine right-handed) participated in this study. All participants provided written informed consent before participation. All experiments were approved by the local ethics committee and in accordance with the declaration of Helsinki.

Procedure
Experiments were performed in a dim and silent room shielded from external electromagnetic interference. Participants were seated upright in a chair facing a single neutral image. Stimulation was applied on the dorsum of the right hand on five different locations in separate stimulation blocks of 20 sequences each. Each pulse sequence had a duration of 8.5 s (Fig. 1). Before every stimulation block, the detection threshold to a single 0.5 ms pulse was determined using a staircase paradigm and the pulse amplitude was set to twice this detection threshold. After the first stimulus sequence of each block, participants were asked to rate the pain on a visual analog scale (VAS) ranging from 0 (no sensation at all) to 5 (painful) to 10 (worst pain imaginable). The location order was chosen such that none of the locations were located directly next to each other to limit local (de)sensitization of the skin.

Nociceptive stimulation
Participants were electrically stimulated using intra-epidermal electric stimulation at twice the detection threshold with a current controlled stimulator (AmbuStim, University of Twente, Enschede, the Netherlands). Cathodic square-wave pulses (pulse width: 0.5 ms) were applied using an electrode comprised of 5 microneedles in a layer of flexible silicone (Steenbergen et al., 2012) to selectively stimulate nociceptive afferent nerve fibers in the epidermis. Electrode dimensions are displayed in Fig. 2. Each microneedle protrudes 0.5 mm from the electrode surface. Experiments  and simulations (Motogi et al., 2016;Poulsen et al., 2020) show that intra-epidermal stimulation preferentially activates nociceptive afferents, provided that the stimulus intensity remains at or below twice the detection threshold.
A validation study of the electrode used in this study showed that Fig. 1. Experiment outline. Participants were stimulated in 5 blocks of 20 sequences with a randomized interval of 10-15 s. For every block, the intraepidermal stimulation (IES) electrode was moved to a different location on the dorsum of the right hand. Participants were allowed a small break in between every stimulation block. Before the start of every block, the detection threshold was measured which was used to set the amplitude of the pulse sequence to twice the detection threshold.
stimulation resulted in a sharp pricking sensation (Steenbergen et al., 2012). Average response time and evoked potential latency recorded using this electrode were similar to the average response time and evoked potential latency in previous studies using intra-epidermal stimulation (van den Berg and Buitenweg, 2021).

Multisine frequency modulation
Stimulation was controlled by a microcontroller connected to the trigger input of the stimulator. Each trigger pulse generated by the microcontroller resulted in a single stimulation pulse. Frequency modulated trigger sequences were generated by modulating the inter-pulse interval, which was computed based on (multi)sine frequency modulation, see Eq. (1).
Modulation frequencies (f 1 , f 2 , f 3 ) were chosen such that measured SSEPs are representative of the behavior of the studied sensory system. In a previous study using square wave intra-epidermal stimuli (Colon et al., 2012a) brain responses were measured in a range from 3 to 43 Hz, where lower frequencies resulted in a more consistent response. To improve signal-to-noise ratio, frequencies with a large interference of EOG artefacts and alpha waves should be avoided. For multisine modulation, frequencies should be chosen such that the number of overlapping (sub)harmonics is minimized in order to apply nonlinear system identification techniques to the measured SSEP (Yang et al., 2016). In this study, 3, 7 and 13 Hz were used as modulation frequencies.
To avoid any transient brain activity to individual pulses within a sequence, the maximum inter-pulse interval should be well below the minimum detectable inter-stimulus interval. For nociceptive laser stimuli, stimuli can be individually perceived with an inter-stimulus interval as small as 200 ms (Lee et al., 2009). As this interval could be even shorter for intra-epidermal stimuli, we chose to limit the minimum pulse frequency to 20 Hz (i.e. 50 ms inter-pulse interval). Furthermore, we chose a maximum pulse frequency of 200 Hz (i.e. 5 ms inter-pulse interval) to limit the effects of peripheral nerve repolarization on measured SSEPs. The modulation amplitude (A1 A2 A3) is limited by practical constraints. The minimum pulse frequency (min(f pulse (t))), which is dependent on the offset (C offset ) and modulation amplitudes (A 1 ,A 2 ,A 3 ), cannot be smaller than can be measured given the duration of the stimulus sequence (T stim ), i.e. min(f pulse (t)) > 1 Tstim . Observed power increases with modulation frequency. The observed power is lower, and less accurate, as the modulation frequencies (f 1 ,f 2 ,f 3 ) get closer to the minimum pulse frequency (Fig. 3). For a multisine waveform the observed power at each frequency is dependent on the combination of frequencies (f 1 , f 2 , f 3 ), modulation amplitudes (A 1 , A 2 , A 3 ) and phases (ϕ 1 , ϕ 2 , ϕ 3 ). Simulations indicated that the selected frequencies 3, 7 and 13 Hz are a good tradeoff between observed power and the physiological constrains (20 Hz < f pulse < 200 Hz) with an offset (C offset ) of 110 Hz, a modulation amplitude of 30 Hz for each frequency and phases of + 1 3 π and − 1 3 π for 7 and 13 Hz respectively (Fig. 3). Also note that as the total available bandwidth (90 Hz) is divided over three modulation frequencies, this leads to a reduction of power to one-ninth of the power when the full bandwidth is used for a single frequency (Fig. 3).  Simulation results of the observed power of a modulated pulse sequence when the modulation frequency is increased from 0 to the minimum pulse frequency (left), when the modulation amplitude is increased from 0 to the frequency offset when modulating at a single frequency (middle left) or at a combination of three frequencies (middle right). The observed power increases as the modulation amplitude increases. However, the observed power decreases as the stimulation frequency gets closer to the minimum pulse frequency. The last figure (right) shows the power spectrum of pulse sequences modulated at a frequency of 3 Hz with a modulation amplitude of 90 Hz and a modulation amplitude of 30 Hz. As the total modulation amplitude is divided over three frequencies in a multisine waveform, the total modulation amplitude per frequency is reduced to one-third and the power is reduced to one-ninth.

EEG recording
Scalp EEG was recorded with a sampling rate of 1024 Hz using a REFA amplifier (TMSi B.V., Oldenzaal, the Netherlands) and 128 Ag/ AgCl electrodes located according to the international 10/5 system (Oostenveld and Praamstra, 2001) with a common average reference. Electrodes were gelled with an impedance below 10 kΩ. A tubular net bandage was placed over the EEG cap for improved fixation of the electrodes and reduction of potential movement artefacts. Impedances were checked before each stimulation block and decreased if necessary by adding gel. Channels of which the impedance did not decrease to a value below 10 kΩ were disconnected.

Data preprocessing
The recorded EEG was pre-processed using EEGlab (Delorme and Makeig, 2004), a Matlab toolbox for EEG signal processing and used for identification of stimulation artefacts and SSEPs. For identification of stimulation artefacts, the EEG recording was high-pass filtered with a cutoff frequency of 60 Hz, and extracted in epochs from -10 to 30 ms with respect to each pulse applied during the experiment. For identification of SSEPs, the EEG recording was high-pass filtered with a cutoff frequency of 0.5 Hz and low-pass filtered with a cutoff frequency of 40 Hz. Channels were re-referenced to common average after removing channels with flat or excessive EMG activity. Epochs were extracted from 0.5-8.5 s with respect to stimulation onset. Epochs with excessive EMG activity or eye movement artefacts were removed by visual inspection. Subsequently, any residual contamination of the EEG by eye blinking, eye movements or EMG activity was removed using adaptive mixture independent component analysis (Palmer et al., 2008).

Identification and analysis of stimulation artefacts
To inspect the stimulation artifact, the average over all epochs in all sequences and blocks was computed for each participant, with approximately 80 000 epochs per participant as an epoch was extracted around each pulse. Consequently, we were able to detect and map the stimulation artefact on a nanovolt scale. Based on this analysis, one participant was excluded due to an excessive stimulation artefact, defined as a stimulation artefact larger than 100 nV at the Cz channel.

Identification and analysis of steady-state evoked potentials
For spectral analysis, epochs were divided in segments of two seconds, and each segment was Fourier transformed. To identify potential SSEPs on an individual level the T 2 circ value (Victor and Mast, 1991) was computed for every channel on every stimulated frequency using Eq. (2).
The T 2 circ is the ratio between the power ( ) and the variance (σ 2 (f)) of the segments in the frequency domain (X m (f)) multiplied by the number of segments (M) minus one, and is therefore a scaled version of the signal-to-noise ratio. This value is distributed according to the F-statistic and can be statistically tested accordingly (Victor and Mast, 1991).
On a group level T 2 circ was averaged over all participants to determine the average scalp distribution of observed SSEPs. In each subject and on a group level, the 7 central midline electrodes (C5, C3, C1, Cz, C2, C4, C6) were tested for significance by testing the T 2 circ against the F-statistic with a significance level of 0.05 as these are the locations where most nociceptive SSEP activity was observed in previous studies (Blöchl et al., 2015;Colon et al., 2012b;Mouraux et al., 2011). Subsequently, we computed the power [|X(f) ⃒ ⃒ ⃒ 2 ], phase [Arg(X(f))] and noise level [ σ 2 (f ) M ] of time-locked activity at all (sub)harmonic frequencies for all contralateral central midline electrodes with a significant T 2 circ at one of the base frequencies.
To estimate the overall time delay, the time delay of each significant frequency at the corresponding electrode was estimated based on the phase. The time delay (τ) was computed based the phase delay (ϕ) using Eq. (3). As phases wrap over a period of 2π and SSEP responses could be both positive and negative (Norcia et al., 2015), time delay estimates are repeated at k 2f , in which f is the stimulated frequency and k is an integer number

Pain rating
Nociceptive detection thresholds remained between 0.1 and 0.5 mA ( Fig. 4) with an average of 0.43 ± 0.17 mA. Pain ratings were reported on a VAS Scale ranging from 0 (no sensation at all) to 5 (painful) to 10 (worst pain imaginable). The VAS line mentioned 0, 5 and 10 and measured 12 cm in length. One of the nine included participants did not report pain ratings. Participants reported a continuous mild sensation with an average VAS score of 3.0 ± 1.9 (Fig. 4). A single subject did report a strong sensation of pain (VAS > 7) in the 1st, 3rd and 5th block. Participants also reported that the perceived intensity of each sequence decreased across sequences within each block. Redness of the skin around the needle locations was observed after each stimulation block.

Stimulation artefacts
A stimulation artefact was observed in all participants (Fig. 5). In 8 out of 9 participants, the stimulation artefact was concentrated around the ground electrode on the right mastoid, indicating that the observed stimulation artefact was caused by displacement currents toward this electrode. In all included participants the stimulation artefact at Cz was limited to a maximum of 50 nV and occurred between 0 and 10 ms after each pulse. The topographical distribution of the artefact is ipsilateral to the side of stimulation and therefore contralateral to any expected brain activity.

Multisine SSEP topographies
The group and individual level topographies of the T 2 circ statistic at 3, 7 and 13 Hz are shown in Fig. 6. At 3 Hz, a total of 7 participants showed a significant SSEP response on at least one of the central midline electrodes. At 7 Hz a significant SSEP response was observed in 4 participants and at 13 Hz only in 3 participants. At 13 Hz, significant electrodes were located ipsilateral in 2 participants. On a group level, significant 3 Hz and 7 Hz spectral components were mostly observed central/midparietal and slightly shifted towards the contralateral side, centered around the Cz, C1, C3 and C5 channels. At 13 Hz, the group level topography shows activity in similar areas as 3 and 7 Hz in addition to other regions ipsilateral with respect to stimulation.

Multisine SSEP spectra
Of the central midline electrodes, a significant spectral component of one or more of the base frequencies was observed at C5, C3, C1, Cz, C2 and C6. At the contralateral electrodes (C5, C3, C1 and Cz), the spectra were analyzed to investigate the base frequencies and (sub)harmonics, see Fig. 7. At 3 Hz, 3 participants had a significant spectral component at Cz, 4 participants had a significant spectral component at C1 and C3, and 2 participants had a significant spectral component at C5. Furthermore, the group level averages at Cz, C1 and C3 had a significant spectral component (p < 0.05). At 7 Hz, 3 participants had a significant spectral component at C5 and with a significant group level average (p < 0.05) at the same electrode. At 13 Hz, 1 participant had a significant spectral component at Cz, C1 and C3 and no significant spectral component was observed at C5. At all other frequencies, a few significant individual spectral components were incidentally found on 2nd and 3rd order harmonics at both electrodes.

Time delay
Time delay was estimated for 3 and 7 Hz in individual participants and on a group level at the electrodes with significant SSEP on a group level (C2, Cz, C1 and C3 for 3 Hz, C5 for 7 Hz). As for 13 Hz, no electrodes were significant on group level, this frequency was not used for time delay estimation. For the group level estimate, only the participants with a significant spectral component at those frequencies were averaged.
The resulting estimates are shown in Fig. 8. Multiple subject and group level estimates are displayed in each column due to the 2π wrapping effect. Group level estimates align around an average time delay of 168 ms.

Discussion
This study outlined a method to directly stimulate nociceptive afferents in the skin to evoke SSEPs in order to study nociceptive processing. In contrast to previous studies, using square waveforms at a single frequency to modulate electrical stimulation (Colon et al., 2012a;Blöchl et al., 2015), we used a multisine waveform (Fig. 9). The use of specific combinations of base frequencies in this multisine waveform allows for system identification techniques to explore system properties such as delay, signal-to-noise ratio and (non)linearity.

Technical challenges
A number of challenges were addressed to reliably measure SSEP brain activity in response to intra-epidermal electric stimulation. Low pulse amplitudes were a requirement for preferential activation of superficial nociceptive afferents in the skin leading to a low signal-to-noise ratio of any observable brain activity. In addition, the total available bandwidth for frequency modulation was divided in three to achieve multisine modulation, leading to a ninefold decrease of the power of single stimulation frequencies in the input signal. As such, the implementation of the paradigm in terms of hardware and procedures was optimized to reduce noise and augment any potential SSEP activity. This involved 1) optimizing the temporal accuracy of stimulation hard-and Fig. 4. Nociceptive detection thresholds (NDTs) in mA and pain ratings on a visual analog scale (VAS) ranging from 0 (no sensation at all) to 5 (painful) to 10 (worst pain imaginable) at each stimulation block (see also Fig. 1). Detection thresholds remained between 0.1 and 0.5 mA. Participants were stimulated at twice the NDT. In general, participants reported a mild sensation with an average score of 3 in response to this stimulation.

Fig. 5.
Maps and amplitude of the stimulation artefact. Most stimulation artefacts were centered around the ground electrode on the right mastoid. As such, the observed stimulation artefact was likely due to displacement currents. The stimulation artefact at Cz remained limited to below 50 nV in all included participants and occurred between 0 and 10 ms after each pulse. software, 2) strategically choosing stimulation parameters, 3) reducing stimulation artefacts by strategic electrode placement and lowering electrode impedances and 4) reducing EOG, EMG and movement artefacts as much as possible. We observed that incorrect implementation of any of these steps could lead to a substantial loss of power of the observed SSEP and an increase of background noise during several pilot recordings.
Special care was taken to identify and reduce potential stimulation artefacts. Stimulation artefact would have a similar frequency content as the input signal and overlap with potential brain activity related to stimulation. Three potential sources of electric stimulation artefacts are identified in literature: volume conduction current, displacement current and electromagnetic coupling (McLean et al., 1996). In this study, the volume conduction current artefact was reduced by limiting the electric field of the stimulation (i.e. using multiple intra-epidermal needles, scrubbing the skin, using a low stimulation intensity). Electromagnetic coupling was reduced by lowering the impedance of EEG scalp electrodes as much as possible, and keeping stimulation and EEG electrode leads as far apart as possible. The displacement current artefact (Fig. 10) was reduced by scrubbing the location of stimulator ground before application and using a large surface TENS electrode as a ground proximal with respect to the stimulation electrode. As displacement currents could nevertheless occur, the EEG ground was placed on the mastoid bone ipsilateral to stimulation to prevent any interference with brain activity, which was expected on central or contralateral locations.
Stimulation artefacts are usually much faster than the dynamics of neural systems. While an electric stimulation artefact usually occurs in the range of a few milliseconds (McLean et al., 1996) a neural response takes several tens to hundreds of milliseconds. We made use of this property to design an accurate procedure for identification of potential stimulus artefacts by averaging the EEG with respect to each pulse. The resulting waveforms showed that a small stimulation artefact remained in most participants with an amplitude of up to 50 nV at Cz and concentrated around the EEG ground electrode on the right mastoid. Fortunately, it is unlikely that the stimulus artefact interferes with any observed central and contralateral brain activity due to the small amplitude at Cz and its ipsilateral distribution.

Multisine SSEP responses
We observed brain activity in response to multisine intra-epidermal stimulation with significant spectral components at the base frequencies 3, 7 and 13 Hz. On average, stimulation was perceived as a mild sensation, rated 3.0 on a VAS. Stimulation caused a slight redness of the skin around needle locations after each stimulation block, which was potentially caused by the release of substance P by stimulated nociceptors. No other signs of potential tissue damage were observed. Fig. 6. The group and individual level topographies of the T 2 circ statistic at 3, 7 and 13 Hz. The SSEP response was tested for significance at the central midline electrodes (C5, C3, C1, Cz, C2, C4, C6). Significant electrodes (p < 0.05) are indicated by green dots. At 3 Hz, a significant SSEP response was observed in 7 participants on at least one of the central midline electrodes. The group level average has a central/mid-parietal distribution and is significant at C2, Cz, C1 and C3. At 7 Hz, a significant SSEP response was observed in 4 participants on at least one of the central midline electrodes. The group level average also has a central/mid-parietal distribution and is significant at C5. At 13 Hz, a significant SSEP response was observed in 3 participants on at least one of the central midline electrodes. The group level average shows activity at the same area as 7 Hz stimulation and several other areas, with significance at C6.
Crucially, the stimulus intensity was high enough to be clearly perceived by most participants, but low enough to allow the participant to relax during the stimulation and EEG recording. At 3 Hz, we observed significant central midline electrodes for a majority of the participants and a clear group level topography. At 7 Hz and 13 Hz, we observed significant central midline electrodes in four and in three participants respectively. There is a considerable variation in topographies among individual participants, which could be due to a larger amount of background activity in some of the participants. However, this might also reflect anatomical variation and the large diversity in which participants tend to respond to nociceptive stimulation.
To the best of our knowledge, this is the second study to demonstrate SSEP responses to intra-epidermal stimulation. The first demonstration of this technique was done by Colon et al. (2012b) where SSEP responses were shown at 3, 7, 13, 23 and 43 Hz. Similar to the study of Colon et al., we observed SSEP responses to 3, 7 and 13 Hz intra-epidermal stimulation. However, the observed signal topographies are markedly different. In the study of Colon et al. a frontal topography was observed which centers around AFz as the stimulation frequency is increased. Another recent study by Blöchl et al. (2015) used a concentric planar electrode in an attempt to generate nociceptive specific SSEPs with 31 and 37 Hz stimulation. They reported contralateral activation similar to the topographies observed in this study for 7 and 13 Hz stimulation.
In both these studies it remains unclear whether stimulus artefact could have contributed to observed topographies. In the first study, the   9. Examples of frequency modulated intraepidermal stimulation using pulse sequences. The aim is to preferentially stimulate nociceptive afferents in the skin and therefore measure SSEPs related to nociception. Earlier studies aimed at generating SSEPs by stimulation of nociceptors using square wave modulation (left). Theoretically, this type of stimulation does not only elicit the stimulated frequency, but also its harmonics. In this work, we proposed to stimulate nociceptors using sinusoid or multisine frequency modulated pulse sequences. Fig. 10. Examples of two stimulation artefacts observed in pilot sessions. In A) and B) the artefact was concentrated around the EEG ground electrode, which was located on the forehead and right mastoid respectively. As such, these artefacts were likely caused by the displacement current. fact that the topography was centering around the ground of the EEG cap as stimulation frequency increased might suggest that not all stimulation artefact was removed by the independent component analysis algorithm used for artefact removal. Furthermore, the average stimulus intensity in that study was 0.23 (±0.08) mA applied on a single intra-epidermal needle, while the average stimulus intensity in this study was 0.43 (±0.17) mA applied over five intra-epidermal needles, effectively resulting in a much lower current density and smaller electric field in this study. In the second study, no method of dealing with stimulation artefacts was reported.
A recent study to observe purely nociceptive SSEP topographies was done by Mouraux et al. (2011). They used laser stimulation, which selectively activates nociceptive afferents in the skin and generates no electric stimulation artefacts. In response to 7 Hz periodic stimuli, they reported a topography centered around Cz similar to the topography observed in response to 3 Hz stimulation in this study. Although stimuli where applied through intra-epidermal needles, the observation of contralateral activation in response to 7 Hz and 13 Hz stimuli in the current study could imply concurrent activation of tactile large diameter Aβ-fibers in the skin. As it remains uncertain whether this difference is caused by concurrent activation or simply by the fact that the brain responds differently to a different stimulus modality, the intra-epidermal stimulation in this study should be considered 'preferentially' but not 'specifically' nociceptive.

System behavior
To the best of our knowledge, this is the first study demonstrating SSEP responses to multisine frequency modulated electric stimulation of superficial afferents in the skin. Compared to previous studies using square wave modulation at a single frequency, the multisine modulation used in this study allows for exploration of system (non)linearity by observation of the combination of (sub)harmonics associated with the three stimulation frequencies in the power spectrum. The observed brain activity also allows for estimation of other properties such as the delay and order of the stimulated system. An estimate of time delay was obtained by unwrapping the phase delay of the SSEP at multiple frequencies (Norcia et al., 2015). The average time delay of participants with significant spectral components at the significant base frequencies aligned around a latency of 168 ms. This is markedly later than the latency observed in an earlier study using laser stimulation by Mouraux et al. (2011), where delays between 30 and 50 ms were observed in response to 7 Hz stimulation. However, it should be noted that Mouraux et al. used the average potential over each stimulation period to compute the time delay. This implies that the true value of the time delay observed in that study could also be ± k larger due to phase wrapping (where k = 0, 1, 2, …), indicating that their time delay could actually be similar to the one measured in this study. The observed delay corresponds with the latency of the N2 component observed in evoked potentials in response to intra-epidermal electric stimuli during earlier studies (van den Berg et al., 2020;van den Berg and Buitenweg, 2021), and supports that the observed signal is indeed related to brain activity. The distribution of signal power in base frequencies and (sub)harmonics provides information on the linearity of the system. Studying system (non)linearity could have important implications for the way we model nociceptive processing. While base frequencies were significant, we did not find any significant (sub)harmonics on a group level. As such, we did not find any indications that nociceptive processing is nonlinear. Nevertheless, further studies replicating these findings in a larger sample size are required to determine if nociceptive processing could be modeled as a linear system.

Conclusion
Intra-epidermal stimulation of superficial nociceptive afferents in the skin using a multi-sine modulated pulse sequence of 3, 7 and 13 Hz, results in SSEPs at central electrodes. Significant SSEPs at the base frequencies 3, 7 and 13 Hz were found in a majority of participants. Such multisine SSEPs can be used to study the temporal dynamics of nociceptive processing in terms of delay and nonlinearity. Phase analysis indicated an average time delay of 168 ms. No indications for nonlinearity of nociceptive processing were observed in the current exploratory dataset.