Three age groups of subjects participated in the study. Thirty-one healthy typically developing preschoolers aged 4 to 5 years (right-handed, 16 males) and 30 healthy typically developing schoolchildren aged 8 to 10 years (right-handed, 15 males) participated in the study. Also, 30 healthy adult volunteers between 20 and 40 years of age (right-handed, 14 males, mean age 29.3) participated in the experiment. All participants took the Wechsler intelligence test (see Table 1).

All participants (and parents of children) were interviewed before study by clinical psychologist and physician to exclude the development, psychological or neurological disorders, to check handedness, and to describe the design of the experiment and the type of stimuli. Exclusion criteria were as follows: history of anxiety, depressive, personality, or psychotic disorders, post-traumatic stress disorder (PTSD), sleep disorders, attention-deficit hyperactivity disorder (ADHD), and other trauma-related mental disorders; autism/social communication disorders; language dysfunction; tactile hypersensitivity; and ethnic minorities, abnormalities during EEG registration, left-handed or ambidexter.

Prior to the study, both children and their parents received an introduction to the experimental procedures and had the opportunity to familiarize themselves with tactile stimulation several days before the experiment. If any of the participants showed intolerance to cold or other tactile stimuli, the study was not conducted for their safety and comfort. The study was performed following the ethical standards reported in the 1964 Declaration of Helsinki and its later amendments or comparable ethical standards. The Ethics Committee approved the study protocol of the Pushkin State Russian Language Institute (Moscow, Russia) (protocol code 17-3-24-118, date of approval 15.07.22). All individuals provided written informed consent after receiving a complete description of this study.

Subjects were presented with 3 types of stimuli (Fig. 1A):

(1) slow stroking with a soft brush (2–5 cm/s) with medium pressure;

(2) slow stroking with a bristle brush (2–5 cm/s) with medium-strong pressure;

(3) slow stroking with a specialized ice brush (2–5 cm/s).

The stimuli were selected basing on our previous studies, demonstrated which type of stimuli induced pleasant or unpleasant response in children in adults [21] and the pilot study described below.

The pilot study consisted only on stimuli assessment and recruited 38 healthy adults from 20 to 40 years old (20 males, mean age 30.6). Participants assessed the stroking of nine stimuli types (three pleasant, three unpleasant and three neutral), which included two surfaces of soft brush, two surfaces of bristle brush, neurological wheel, hot bags, ice, vibration, textile surface. The stimuli were presented with slow velocity and varied pressure (from medium to strong). The stimuli were assessed by the following scale: pleasantness, pain, arousal, and tickling. As a result, three stimuli were selected: soft brush with artificial squirrel texture with medium pressure (most pleasant stimulus with mean pleasantness 7.7), bristle artificial plastic brush with medium-strong pressure (as most neutral stimuli with mean pleasantness 4.9) and the ice with medium pressure (as most unpleasant stimulus with mean pleasantness 2.1)

After the main phase of the study and a short break, each stimulus was presented once, and subjects (including both adults and children) were asked to rate the following parameters on a scale of 1 to 9: pleasantness, pain, arousal, and tickling.

EEG was acquired using a 19-channel EEG amplifier Encephalan with the recording of polygraphy channels (Poly4, Medicom MTD, Taganrog, Russia). The sampling rate was 250 Hz. The amplifier bandpass filter was nominally set to 0.05–70 Hz. AgCl electrodes (Fp1, Fp2, F7, F3, Fz, F4, F8, T3, C3, Cz, C4, T4, T5, P3, Pz, P4, T6, O1, and O2) were placed according to the International 10–20 system. The electrodes placed on the left and right mastoids served as joint references under unipolar montage. The vertical electrooculogram (EOG) was recorded with AgCl cup electrodes (Poly4, Medicom MTD, Taganrog, Russia) placed 1 cm above and below the left eye, and the horizontal EOG was acquired by electrodes placed 1 cm lateral from the outer canthi of both eyes. The electrode impedances were kept below 10 k$\mathrm{\Omega }$.

Each subject’s EEG recordings consisted of 3 minutes of resting state EEG slices with eyes closed (approximately 30 seconds) and eyes open (2–2.5 minutes), followed by approximately 30 minutes of EEG with tactile stimulation. To eliminate unwanted artifacts such as oculomotor, electrocardiogram (ECG), and muscle interference, independent component analysis was applied to the EEG recordings using EEGLab 14_1_2b software (Matlab toolbox, Swartz Center for Computational Neuroscience, San Diego, CA, USA) [36]. Visual analysis was then performed and any remaining artifacts were manually removed, ensuring that no more than 5 seconds of the background recording and no more than 1% of the EEG recordings contained any labels.

The resulting data were fragmented according to stimulus type (labeled as Presentation markers). Further analysis of the EEG data included these fragments. To calculate linear and nonlinear EEG parameters, 10-second intervals for each type of tactile stimulation were used (10 $\times$ 20-second intervals of each type). The EEG was also analyzed in the resting state with eyes open (all measurements were made on average in the resting state before stimulation)—the interval of the analyzed fragment was 2–2.5 minutes.

All EEG characteristics (see below) were calculated for each electrode and then averaged over the frontal (F3, F4, Fz), central (C3, C4, Cz) and parietal areas (P3, P4, Pz, O1, O2).

Power Spectral Density (PSD) of EEG was analyzed using Matlab R2018b program (Matsoft, Natvik, MA, USA). Fast Fourier transform (FFT) was used to analyze the power spectra of the studied EEG fragments. The EEG spectrum was estimated for each studied condition. The obtained power spectra were averaged in 1 Hz steps in the range from 2 to 20 Hz (2–3 Hz, 3–4 Hz, … 19–20 Hz). Logarithmic power values were used for further statistical analysis.

Fractal dimension D2 was calculated for the entire frequency range (2–20 Hz). The mathematical computation of fractal dimension D2 was performed using the Higuchi method [37], based on the transition from the original sequence to sequences obtained by decimation, i.e., taking every k-th observation. For each selected value of k, k such independent series were generated, differing in the shift of the starting point. For each of them, the curve length was computed as the sum of absolute differences between data points, normalized to the length of the original curve. The arithmetic mean of these calculated lengths was then computed across all k series taken with a step size of k, resulting in an averaged estimation of the length L(k). The FD was calculated using tool prepared in C++ environment.

For series with a fractal dimension D, the function for series with a fractal dimension D, the function $L(k)\mathrm{~}{k}^{-D}$is used. To estimate the value of D, linear regression of the quantity logL(k) on k is employed.

The Envelope Mean Frequency (EMF) was calculated using tool prepared in C++ environment. We calculated the EEG signal envelope using the Hilbert transform and analyzed the envelope frequency in the 2–20 Hz frequency range. To calculate the Hilbert envelope frequency, the Hilbert transform y(t) of the function x(t) was first calculated using the following formula.

$$y(t)=\frac{1}{\pi }\left({\int }_{t-1/\epsilon }^{t-\epsilon }\frac{x(\tau )}{t-\tau }𝑑\tau +{\int }_{t+\epsilon }^{t+1/\epsilon }\frac{x(\tau )}{t-\tau }𝑑t\right)$$

The envelope was calculated as

$$A(t)=\sqrt{x^2(t)+y^2(t)}$$

and then was subjected to Fourier transformation to calculate its spectrum. In addition, the average frequency of the Hilbert envelope will be calculated as follows:

$$\underset{¯}{f}=\frac{{\sum }_{i=1}^{k}a(i)f^2(i)}{{\sum }_{i=1}^{k}a(i)f(i)}$$

This frequency ($\underset{¯}{f}$) represents the average frequency using instantaneous information (frequency f and amplitude a within the range from 1 to k).

We calculated the HC parameter in the 2–20 Hz range:

HC = $\frac{\text{mobility}\left(y^{\prime }(t)\right)}{\text{mobility}(y(t))}$, $\mathrm{\copyright }$™$\mathrm{\Sigma }$ mobility = $\sqrt{\frac{\mathrm{var}\left(y^{\prime }(t)\right)}{\mathrm{var}(y(t))}}$, where y(t) is EEG signal, while y’(t) represents the variability of the EEG signal. The HC was calculated using tool prepared in C++ environment.

Statistical analysis was performed using Statistica13 software (StatSoft Inc, Hamburg, Germany). To compare the effect of tactile stimulation (three types of stimuli and resting state) on EEG dynamics between groups, we used repeated measures analysis of variance (ANOVA) with Bonferroni correction (p 0.05). We examined the normality of subjective stimulus ratings. We analyzed the EEG data for three brain areas: the frontal (F3, F4, Fz), central (C3, C4, Cz) and parietal areas (P3, P4, Pz, O1, O2) as independent sources and prepared the repeated measures ANOVA to find differences at the topographic level. The Bonferroni correction for 3 brain areas was also applied. If ratings passed the Kolmogorov-Smirnov normality test, we further conducted a repeated-measures ANOVA to determine group differences in subjective stimulus ratings. However, if subjective ratings did not meet the normality test in either group, we used the Mann-Whitney method. The Spearman rank correlation test was used to analyze the relationship between subjective evaluations and EEG parameters. We also applied correction for multiple comparisons.

We found no significant group differences between the subjective evaluation of stimuli on the scales: the pleasantness of each stimulus (the analysis was made separately for subjective assessment of pleasant, unpleasant, and neutral stimuli) differed neither between children, nor between children and adults (p $>$ 0.33). The absence of significant differences between groups of participants was also revealed for scales pain, arousal, and tickling (p $>$ 0.45). At the same time, we noted a greater variability in the evaluation of pleasant and neutral stimuli in children compared to adults on all scales except pleasantness. The mean scores ($\pm$ standard deviation (SD)) are presented in Table 2, the SD values differed between groups marked with bold font.

The higher SD in the perception of pleasant tactile stimulation (CT-optimal stroking) was accompanied by a deviation in subjective evaluations among adults. Fig. 2 depicts the distribution of subjective assessments in children in adults for the scale pleasantness. While in children the distribution of subjective pleasantness from stroking with a soft brush was normal (Fig. 2), in adults’ subjective evaluations of pleasant stimuli could be divided into pleasant and neutral (Fig. 2c). Moreover, the variability in subjective pleasantness of pleasant stimulus in healthy adults was associated to the EEG dynamics (Fig. 2d,e), namely the decrease in PAF during pleasant stimulation compared to the resting state (see section “Correlation Analysis”).

An increase in the theta rhythm (4–8 Hz) compared to the resting state was found in the central regions of both schoolchildren and adults, but only when an unpleasant tactile stimulus was presented. In contrast, preschoolers showed a decrease in the theta rhythm specifically for unpleasant stimuli (F(4, 168) = 15.645, p = 0.0000).

A significant decrease in alpha rhythm power in the range of 8–10 Hz was observed in schoolchildren when presented with all types of stimuli in frontal and central areas. However, in preschoolers, a significant decrease was observed only for pleasant and unpleasant tactile stimulation (F(4, 168) = 14.109, p = 0.0000). In adults, a significant decrease in alpha rhythm frequency in the range of 8–10 Hz was observed during both neutral and unpleasant tactile stimulation compared to the resting state in frontal, central and parietal areas (Fig. 3).

In schoolchildren and preschoolers, a significant increase in EEG power at 10–12 Hz was observed only during unpleasant stimuli compared to the resting state in central and parietal areas. A decrease in alpha2-rhythm power with same topography was observed in schoolchildren only when pleasant tactile stimuli were presented. In adults, EEG power at 10–12 Hz was significantly higher only for unpleasant stimuli, and a decrease in alpha2 rhythm was observed for pleasant and neutral stimuli (F(4, 168) = 12.346, p = 0.0000) incentral and parietal areas.

An increase in beta rhythm in the range of 16–20 Hz was observed in preschoolers only when presented with pleasant and unpleasant tactile stimuli compared to the resting state in central and frontal areas. In schoolchildren and adults, this increase in beta rhythm was observed for all stimulus types (Wilcoxon Matched Pairs Test, p 0.0044).

In adult volunteers, a significant increase in PAF in central and parietal regions was observed during the presentation of neutral and unpleasant tactile stimulation compared to the resting state. The data were depicted in Fig. 4a. However, no significant differences in PAF were found when pleasant tactile stimulation was presented (mixed stimuli group (3) stimuli (4) effect F(6, 264) = 18.234, p 0.0001, partial eta-square 0.29; post hoc Bonferroni p 0.0088). A significant increase in PAF was observed in school children for all stimulus types (p 0.0029) in frontal and central areas, whereas in preschoolers a significant increase in PAF in the same brain areas was found only when unpleasant tactile stimulation was presented (post hoc Bonferroni p = 0.0021).

In adults, as in schoolchildren, there was a decrease in EEG power in the central and parietal regions during all types of tactile stimulation (F(6, 264) = 14.466, p 0.0001, partial eta-squared 0.25; post hoc Bonferroni p 0.0037). However, a significant decrease in EEG power in preschoolers (only in the parietal area) was observed only during unpleasant tactile stimulation (post hoc Bonferroni p = 0.0015), while no significant differences were found when neutral and pleasant stimuli were presented (Fig. 4b).

Healthy adult participants demonstrated the increase of the FD during pleasant and unpleasant tactile stimulation compared to the resting state in the central and parietal regions (Fig. 4c). Similar (but less pronounced for the pleasant stimuli) changes of the FD were observed in schoolchildren (F(6, 264) = 7.941, p 0.0001, partial eta-squared 0.19; post hoc Bonferroni p 0.0032). However, in preschoolers, only unpleasant tactile stimulation caused an increase in FD compared to the resting state (post hoc Bonferroni p = 0.0008). The topography of the significant EEG changes in children and adults was the same.

Adults showed an increase in the signal complexity index - EEG entropy (HC) for pleasant tactile stimulation and a decrease for unpleasant stimuli in central and parietal regions (F(6, 264) = 6.620, p 0.0001, partial eta-square 0.13; post hoc Bonferroni p 0.0009). Schoolchildren showed a decrease in this index for both pleasant and unpleasant stimuli (post hoc Bonferroni p 0.0026) in the central and parietal regions while in preschoolers significant changes were observed only for unpleasant stimuli (post hoc Bonferroni p = 0.0017) in the parietal area. The results were depicted in Fig. 4d.

Correlation analysis revealed a significant association between the decrease in PAF during pleasant stimulation compared to the resting state in the group of healthy adults (Spearman rank correlation r = –0.65, p = 0.0001) in central and frontal regions (Fig. 2d). No significant correlation between PAF dynamics and subjective pleasantness of stimuli was found in any of the pediatric groups.