Systemic neurophysiological entrainment to behaviorally relevant rhythmic stimuli

Physiological Reports. 2024;12:e70079. | 1 of 16 https://doi.org/10.14814/phy2.70079 wileyonlinelibrary.com/journal/phy2 1 | INTRODUCTION Oscillations are fundamental organizing phenomena in biological systems, providing functional advantages like predictive capabilities, energyefficiency, and enhanced communication through noise resilience and desensitization prevention (Glass, 2001; Rapp, 1987; Xiong & Garfinkel,2023). In the human body, numerous physiological processes, including heartbeat, respiration, and neural activity, exhibit oscillatory behavior and often interact with Received: 23 July 2024 | Revised: 9 September 2024 | Accepted: 24 September 2024 DOI: 10.14814/phy2.70079 ORIGINAL ARTICLE Systemic neurophysiological entrainment to behaviorally relevant rhythmic stimuli ManuelMuñoz- Caracuel1,2 | VanesaMuñoz1 | Francisco J.Ruiz- Martínez1 | Antonio J.VázquezMorejón2,3 | Carlos M.Gómez1 This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2024 The Author(s). Physiological Reports published by Wiley Periodicals LLC on behalf of The Physiological Society and the American Physiological Society. 1Department of Experimental Psychology, University of Seville, Sevilla, Spain 2Mental Health Unit, Hospital Universitario Virgen del Rocio, Seville, Spain 3Department of Personality, Evaluation and Psychological Treatments, University of Seville, Sevilla, Spain Correspondence Manuel Muñoz- Caracuel, Human Psychobiology Laboratory, Experimental Psychology Department, University of Seville, Calle Camilo José Cela, S/N, Sevilla 41018, Spain. Email: [email protected] Funding information MEC | Agencia Estatal de Investigación (AEI), Grant/Award Number: PID2022- 139151OB- I00; Consejería de Economía, Innovación, Ciencia y Empleo, Junta de Andalucía (Ministry of Economy, Innovation, Science and Employment, Government of Andalucia), Grant/Award Number: P20_00537 Abstract Physiological oscillations, such as those involved in brain activity, heartbeat, and respiration, display inherent rhythmicity across various timescales. However, adaptive behavior arises from the interaction between these intrinsic rhythms and external environmental cues. In this study, we used multimodal neurophysiological recordings, simultaneously capturing signals from the central and autonomic nervous systems (CNS and ANS), to explore the dynamics of brain and body rhythms in response to rhythmic auditory stimulation across three conditions: baseline (no auditory stimulation), passive auditory processing, and active auditory processing (discrimination task). Our findings demonstrate that active engagement with auditory stimulation synchronizes both CNS and ANS rhythms with the external rhythm, unlike passive and baseline conditions, as evidenced by power spectral density (PSD) and coherence analyses. Importantly, phase angle analysis revealed a consistent alignment across participants between their physiological oscillatory phases at stimulus or response onsets. This alignment was associated with reaction times, suggesting that certain phases of physiological oscillations are spontaneously prioritized across individuals due to their adaptive role in sensorimotor behavior. These results highlight the intricate interplay between CNS and ANS rhythms in optimizing sensorimotor responses to environmental demands, suggesting a potential mechanism of embodied predictive processing. KEYWORDS embodiment, entrainment, oscillations, predictive processing, rhythm 2 of 16 | MUÑOZCARACUEL etal. external environmental rhythms, as seen in music, dance, language, sports, and other human activities that require adaptive responses to rhythmic patterns (Charalambous & Djebbara,2023; Greenfield etal.,2021; Kotz etal.,2018). Physiological oscillations vary across timescales, but adaptive behavioral control may rely on their temporal alignment with external inputs, facilitating efficient responses to environmental demands (Penzel et al., 2017; Pezzulo et al., 2015). In this context, when physiological oscillations synchronize in phase and/or frequency with an external oscillatory source, they are said to be entrained to that rhythm (Lakatos etal.,2019; Pikovsky etal.,2001). Entrainment has been proposed to underlie cognitive processes like selective attention by aligning neural phases of high excitability with periodic features of relevant external events, thereby enhancing neural gain for attended sensory inputs (Henry etal.,2014; Lakatos etal.,2008, 2013, 2019; Obleser & Kayser,2019; Schroeder & Lakatos,2009; Zion Golumbic etal.,2013). Neural entrainment has also been shown to improve behavioral performance, reducing reaction times, and increasing accuracy in auditory sensorimotor tasks (Nozaradan etal.,2016; Stefanics etal.,2010). However, less in known about autonomic entrainment. While heart rate has shown inconsistent evidence of entrainment with rhythmic external stimuli (Mütze etal.,2020), respiration has been shown to entrain to auditory external stimuli that involve motor responses, though not during passive listening or with complex rhythms (Haas etal.,1986; Wilke etal.,1975). Respiration is recognized as an active sensory selection process modulating neural synchronization across widespread cortical regions to enhance information processing and cognitive performance (Brændholt etal.,2023; Herrero etal.,2018; Kluger etal., 2021; Maric et al.,2020; Perl et al., 2019; Zelano et al., 2016). Moreover, the breathing phase has been shown to systematically influence cognitive performance (Perl etal.,2019; Zelano etal.,2016), highlighting the bidirectional communication between body–brain signals aimed at adapting to taskspecific demands (Criscuolo etal.,2022). Additionally, from the predictive processing perspective (Clark,2013; Friston,2010; Hohwy,2013), the interplay between bodily and neural rhythms in response to external sensory inputs serves a predictive timing function. It aligns information sampling with bodily rhythms to minimize discrepancies between expected (preferred) and sensed body states through recurrent message passing among different hierarchical levels. This preserves organism integrity and supports adaptive goaldirected behavior (Allen & Friston,2018; Engel etal.,2001; Owens etal.,2018; Parr etal.,2022; Pezzulo etal.,2015; Seth & Friston,2016). Despite evidence of neural entrainment with periodic external events and the role of respiration in brain function, the intricate interplay between brain–body rhythms and external inputs remains poorly understood. Active engagement with sensory stimuli is believed to optimize bodily states for higherlevel information processing (Criscuolo etal.,2022; Parviainen etal.,2022), but integrated experimental findings using multimodal neurophysiological measures are scarce. Moreover, most studies exploring bodily entrainment to external stimuli have predominantly utilized simple repetitive stimuli to assess temporal predictive abilities. Incorporating experimental designs that account for the inherent complexities of discriminatory sensorimotor tasks could significantly contribute to understanding the multifaceted dynamics governing brain–body interactions with rhythmic cues. This study aims to shed new light on the complex physiological dynamics of brain and body rhythms in response to rhythmic auditory stimulation, and to elucidate their behavioral implications by differentiating between passive and active sensory encoding processes. Simultaneous recordings of the central and autonomic nervous systems (CNS and ANS) were used during an experimental paradigm featuring temporally invariant rhythmic stimulation. Distinct auditory sequences were presented to dynamically modulate attentional demands on a trialby- trial basis during a cueing task, fostering adaptive speedaccuracy responses that better reflect dynamical sensorimotor contexts. We hypothesize that active engagement with rhythmic auditory stimuli will lead to greater entrainment of CNS and ANS physiological rhythms to the external rhythm, compared to passive listening and baseline conditions. Furthermore, during the active condition, we expect higher accuracy and reduced reaction times for the more prevalent phases of bodily oscillations at the target onset. 2 | MATERIALS AND METHODS We employed a systemic approach to examine various brain and body neurophysiological measures in response to rhythmic auditory stimulation under different attentional conditions. Our analysis included power spectral analysis and phase coherence to compare the alignment of physiological oscillations with the external rhythm across conditions. In addition, we used circular statistics and linear mixed models to explore detailed patterns and their behavioral implications. 2.1 | Participants A sample of 32 healthy participants (12 males and 20 females, 29 righthanded and 3 lefthanded) aged between 19 and 36 years old (mean = 28.31 ± 3.87 SD) were recorded. Two participants were excluded from the analysis due to technical issues during the EEG recording. The 2051817x, 2024, 19, Downloaded from https://physoc.onlinelibrary.wiley.com/doi/10.14814/phy2.70079 by Readcube (Labtiva Inc.), Wiley Online Library on [21/10/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 3 of 16 MUÑOZCARACUEL etal. study was approved by the Bioethical Committee of the Junta de Andalucía. Participants did not report any neurological disease or auditory impairment. The experiments were conducted with the informed and written consent of each participant, following the Helsinki Protocol. 2.2 | Signal acquisition A multimodal signal acquisition was performed, covering different measures of the central and peripheral nervous system at the same time (Figure1a). Electroencephalography (EEG) and intracerebral Functional Near Infrared Spectroscopy (fNIRS- i) recorded the brain electrical activity and oxygenation changes, respectively. Electrodermal activity (EDA), Electrocardiogram (ECG), Plethysmography of the Peripheral Pulse (PPG), Respiratory effort (RSP), and extracerebral Functional Near Infrared Spectroscopy (fNIRS- e) were recorded to assess physiological peripheral signals controlled by Autonomic Nervous System (ANS) activity. FIGURE 1 (a) Multimodal signal acquisition representation (created in BioRe nder. com). (b) Probe layout of fNIRS optodes and EEG electrodes. Red circles indicate fNIRS sources, blue circles fNIRS detectors, and green circles EEG channels. Yellow linear segments between sources and detectors indicate fNIRS standard separation channels. Overlapping small blue circles indicate the location of fNIRS short separation channels. (c) Graphical representation of the experimental design. Three conditions: Baseline, Passive (auditory stimuli + muted film) and Active (auditory stimuli + pressing either the up or down arrow of the keyboard). (d) Auditory stimulation exemplification for the passive and active condition, respectively. S1 and S4 indicate the first (cue) and fourth stimulus (target), respectively, of each auditory sequence. Linear auditory sequences were occasionally broken (20% incongruent trials) in the last stimulus (S4). In the active condition, the expected responses after the arrival of S4 are indicated by black arrows, resulting from the comparison of the auditory frequency of S4 with that of the preceding stimulus (S3) within each trial. 2051817x, 2024, 19, Downloaded from https://physoc.onlinelibrary.wiley.com/doi/10.14814/phy2.70079 by Readcube (Labtiva Inc.), Wiley Online Library on [21/10/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 4 of 16 | MUÑOZCARACUEL etal. 2.2.1 | EEG The EEG was recorded using active electrodes (ActiCAP) from 13 scalp sites mainly distributed along the midline (Figure1b) with a Brain Vision V- Amp DC amplifier (Brain Products, Munich, Germany). The electrode impedance was kept below 20 kΩ. The left mastoid was used as the reference. DC amplification gain was 20,000, and the sampling rate was 1000 Hz. Data acquisition was done using BrainVision Recorder 1.20 (Brain Products). 2.2.2 | fNIRS The fNIRS signal was recorded using a NIRScoutXP device (NIRx Medical Technologies, Glen Head, NY, USA) with 16 LED sources and 30 detectors (14 standard separation detectors +16 short separation detectors) placed along the frontotemporal areas of both sides of the scalp (Figure1b), obtaining a total of 58 channels (42 standard separation channels, 16 short separation channels). Intracerebral (cortical) oxygenation changes (fNIRS- i) were captured by the standard separation channels, whereas extracerebral oxygenation changes (fNIRS- e) on the scalp surface were captured by the short separation channels. The sourcedetector distance was 3 and 0.8 cm for the standard and the short separation channels, respectively. The sampling rate was 3.91 Hz. Data acquisition was accomplished with NIRStar v.14.2 software. 2.2.3 | Peripheral signals The autonomic peripheral signals were recorded using an MP160 (BIOPAC Systems, Goleta, CA, USA) with four amplifier modules (PPG- 100C, EDA- 100C, ECG- 100C, and RSP- 100C), to obtain the PPG, EDA, ECG, and RSP signals, respectively. The peripheral pulse was recorded through a plethysmograph (TSD200 transducer with a wavelength of 860 nm ± 60 nminfrared light) placed on the index finger of the left hand. EDA was obtained using bipolar Ag- AgCl finger electrodes placed on the ring and middle fingers of the left hand. ECG was acquired using three Ag- AgCl lead electrodes (positive, negative, and ground, TSD203 transducer) placed on the left wrist, right wrist, and right ankle, respectively. RSP was recorded through a transducer band placed on the chest (TSD201 transducer), which allows measurement of changes in thoracic circumference that occur with breaths. The amplification gain was set as 100 in PPG, 5 in EDA, 1000 in ECG, and 10 in the RSP amplifiers. The sampling frequency was 1000 Hz. Data acquisition was performed using AcqKnowledge v.5.0.1 software (BIOPAC Systems). 2.3 | Experimental paradigm Participants were comfortably seated in a quiet room in front of a table equipped with a screen and a computer keyboard. The experiment was divided in three consecutive stages or conditions: baseline, passive and active condition (Figure1c). During the baseline condition, participants were asked to stay as quiet as possible for 3 min, looking at a fixed cross in the middle of the screen without any other requirement. In the passive and the active conditions, the same auditory stimulation was used (Figure1d). The stimuli consisted of trials of four consecutive pure auditory tones, presented in an ascending or descending frequency pattern (50%–50%, randomly counterbalanced within participants) that was occasionally broken (incongruent trial) in the last tone of the trial (S4). The auditory frequencies used for the tones were 600, 900, 1200, and 1500 Hz. Each tone lasted 200 ms (with a 20 ms risefall time) and the interstimulus interval (ISI) was 100 ms, resulting in 1100 ms of auditory stimulation, which was followed by 1900 ms of silence, comprising a total of 3000 ms (0.33 Hz) within each trial. Each condition consisted of 240 trials, with an 80%–20% ratio between congruent and incongruent trials, respectively. The experiment was programmed using PsychoPy v.2021.1.2 free software package (Peirce,2009). The auditory stimuli were edited using Audacity free software package, v.3.0.0. Both the passive and active conditions lasted 12 min and shared the same auditory paradigm, but the task instructions in each condition were different. In the passive condition, the participants were asked to watch a muted film and ignore the auditory tones, whereas in the active condition the participants were asked to pay attention to the tones and to respond in each of the trials by pressing the up or down arrow button on the keyboard, depending on whether the fourth auditory tone of the sequence (S4, target stimulus) had a lower or higher auditory frequency than the previous one (S3), trying to combine accuracy and speed in their responses (within a 1900 ms time window before the next trial). The active task was divided into four blocks of 3 min (60 trials), separated by a short break, to prevent the effects of fatigue. There was a training period with an explanatory video, followed by a test trial with visual feedback. The order of conditions, first passive and then active condition, was intentionally kept constant across participants to avoid possible interference or conditioning factors from the active to the passive condition in the auditory processing system. Auditory stimuli were presented through two speakers (Dell, model A215), positioned on either side of the participant's head, at a comfortable volume of 60 dB (measured with a Velleman- DVM1326 level soundmeter). 2051817x, 2024, 19, Downloaded from https://physoc.onlinelibrary.wiley.com/doi/10.14814/phy2.70079 by Readcube (Labtiva Inc.), Wiley Online Library on [21/10/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 5 of 16 MUÑOZCARACUEL etal. 2.4 | Signal preparation and preprocessing 2.4.1 | EEG EEG recordings were imported into EEGLAB v.14.1.1 and MATLAB 2022a software. To correct the EEG for eye blinking, ocular movements and muscle artifacts, an Independent Component Analysis (ICA) was computed. These components were removed, and the EEG signal was posteriorly reconstructed (mean = 3.5 ± 0.93 SD, range 2–7). The signal amplitudes in the electrodes Fp2, Fz9, and Fz10 were used to detect eye blink artifacts and ocular movements, while the 6 channels in the midline (AFz, Fz, FCz, Cz, CPz, and Pz) were the ones accounted for subsequent analysis. 2.4.2 | fNIRS fNIRS raw data were imported into HOMER2 v.2.8 toolbox (Huppert etal.,2009) and MATLAB 2022a software. For the reduction of the signal motion artifacts, the function hmrMotionCorrectionWavelet was applied, with an interquartile range of 1.5. This function has proven to be quite robust in decreasing motion artifacts through wavelet decomposition (Brigadoi et al., 2014; Cooper etal., 2012; Molavi & Dumont, 2012). Additionally, the enPCAFilter function was applied, with one principal component being removed from the data (mean explained variance of the removed first component: 37% ±14%). The function hmrSSR was also applied, which regresses the activity in the short separation channels from the activity of the standard channels to eliminate the contribution of the extracerebral signal from the intracerebral (cortical) signal. The oxyhemoglobin (HbO) concentration changes were obtained employing the modified Beer–Lambert law with a differential partial pathlength factor of 6.0 for 760 nm and 5.0 for the 850 nm wavelength (F. Scholkmann etal.,2010). 3D brain images of the fNIRS channels layout were extracted with the AtlasViewer toolbox (Huppert etal.,2009). 2.4.3 | Peripheral signals Peripheral signals were imported into AcqKnowledge v.5.0.1 (BIOPAC Systems), Ledalab (Benedek & Kaernbach,2010) and MATLAB 2022a software. Signals were visually inspected to check that they were mostly free of movement artifacts. The skin conductance response (SCR) from the EDA signal was extracted with the Ledalab software (Benedek & Kaernbach, 2010) using the Continuous Decomposition Analysis (CDA) function. HR was extracted with the AcqKnowledge software from the R peaks detection of the ECG signal. The ECG and PPG signals were used to calculate the Pulse Transit Time (PTT), a proxy measure of the arterial pressure that can be calculated from the time interval between the R- wave and the peak of the pulse wave at the finger. RSP and PPG raw signals were directly extracted from AcqKnowledge. Respiratory rate (RSP rate) was extracted from the RSP peak detection in AcqKnowledge. All data, except for fNIRS, were downsampled to 100 Hz and detrended by subtracting a bestfit third order polynomial. The HR and PTT tachograms were additionally smoothed using the smooth. m function of MATLAB with the ‘loess’ method, applying a smoothing parameter of 0.01. Data from the passive condition were artificially segmented into four blocks, aligning them with the structure of the active condition for comparative analysis. Consequently, a total of nine blocks (comprising one baseline, four passive and four active conditions), each consisting of 3 min of continuous data, were extracted for each signal and participant for subsequent data analysis. 2.5 | Power spectral density Power spectral analysis for all different signals, participants, blocks and channels (in the case of EEG and fNIRS) was performed using a custom MATLAB script. First, the data matrices were zeropadded by a factor of 10 to increase the frequency resolution of the spectral analysis. Following zeropadding, the FFT was applied to the data. From the FFT result, the onesided absolute values were obtained, since the FFT yields a complex result with both positive and negative frequencies. These absolute values were then squared and normalized by the sampling frequency and number of data points to obtain the power spectral density (PSD). Finally, the resulting power was doubled, except for the DC and Nyquist frequencies. To compare conditions with absolute power differences and focus on relative power changes, the PSD results were normalized dividing the value at each frequency point by the total power. The PSD results of the four data blocks in both passive and active conditions were averaged, resulting in a single PSD signal per condition and participant, covering frequencies from 0 to the Nyquist frequency (50 Hz). Finally, to capture the spectral power at the frequency of interest (auditory stimulation at 0.33 Hz) and compare it between conditions, the integral of the PSD closely around this frequency (0.33 ± 0.005 Hz) 2051817x, 2024, 19, Downloaded from https://physoc.onlinelibrary.wiley.com/doi/10.14814/phy2.70079 by Readcube (Labtiva Inc.), Wiley Online Library on [21/10/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 6 of 16 | MUÑOZCARACUEL etal. was calculated for each condition and participant and exported for statistical analysis. 2.6 | Phase coherence Phase Coherence analysis between all recorded signals and the auditory stimulation for each participant, block, and channel (in the case of EEG and fNIRS) was performed using a custom MATLAB script, enabling the measurement of the consistency of phase differences over time. First, all data matrices were resampled to a common frequency rate of 20 Hz. Then, a zerophase 8th order bandpass IIR filter between 0.2–0.5 Hz was used. The reason for using this frequency range was to capture the frequency of interest (auditory stimulation at 0.33 Hz) at the same time that a potential confounder such as the respiratory frequency band (commonly reported around 0.15–0.5 Hz) containing it. The commonly used lower limit of 0.15 Hz for the respiratory frequency band was not used because it exceeded half the frequency of interest (0.33/2 = 0.165 Hz) and could introduce some harmonics that could alter the phase analysis. To enable the computation of phase coherence between physiological signals and stimuli, a cosine wave at 0.33 Hz was generated from the onset of auditory stimulation in each block. This cosine wave was then included as an additional signal for coherence analysis, representing the auditory stimulation. The continuous analytic signal of all data signals was obtained through the Hilbert transformation, and the phase angle was computed from these signals. Then, the phase difference between each pair of data signals were obtained, converted into complex space through the Euler's formula (eⁱx = cos x + i sin x), averaged, and transformed to absolute values to finally obtain the phase coherence between all pair of signals, ranging its value from 0 (no coherence) to 1 (maximum coherence). Finally, the coherence results from the four data blocks in both passive and active conditions were independently averaged to obtain a single value for each condition, facilitating a direct comparison among the baseline, passive, and active conditions. Likewise, coherence results from all EEG, fNIRS- i and fNIRS- e channels were also averaged for each data signal separately. 2.7 | Phase angle As an additional analysis to explore the specific contribution of distinct phases of the recorded signals to the phase coherence results in the active condition, the phase angles at the arrival of every first stimulus (S1, cue stimulus), fourth stimulus (S4, target stimulus), and response (R), as the drivers of interest, were extracted from the Hilbert transform of the preprocessed signals for each participant, converted to complex space through the Euler's formula and averaged across trials. The resultant mean value of the phase angles for each signal, participant and type of driver (S1, S4, and R) was then exported for subsequent statistical analysis. In the case of EEG, the FCz electrode was the one accounted for phase angle analysis as a central, equidistant, and relatively artifactfree reference point. In the case of fNIRS- i, the channel between source number 14 and detector number 13 (Figure1b), on the right superior temporal gyrus, was the one considered due to the auditory nature of the experiment, whereas in the case of fNIRS- e, the short separation channel adjacent to source number 14 (Figure1b), on the right temporal region, was the one exported for phase angle statistical analysis. Complementarily, to analyze potential linearcircular relationships between reaction times and phase angles, these measures were extracted for each trial, type of trial (congruent or incongruent S4), participant and recorded signal. Subsequently, sine and cosine transformations were applied to the circular variable (phase angles), which allows to capture the periodic nature of the circular variable for statistical analysis. The relationship of hits and errors with phase angles was not finally analyzed due to the marginal existence of errors in the sample, as described in the results section. 2.8 | Statistical analysis Statistical analyses were carried out using SPSS 29.0.1 and R 4.3.2 (R Core Team, 2023) software packages. We first tested the normality of data distribution using a Shapiro– Wilk test and, when the normality assumption was violated, we applied appropriate nonparametric tests. To assess statistical differences between the processed signals, repeatedmeasures analyses of variance (ANOVA) were performed with condition as a factor with three levels: baseline, passive and active. Greenhouse–Geisser correction was used when sphericity was violated. For EEG, fNIRS- i and fNIRS- e, the channels were included as a factor, in addition to condition. Posthoc contrasts (twotailed paired ttest) were carried out when ANOVA differences reached significant results (p < 0.05). To assess the circular uniformity between phase angles of participants for each type of driver (S1, S4, and R) and processed signal, the Rayleigh's test (Jammalamadaka & SenGupta,2001) was performed using the circular package implemented in R. A linear mixedeffects (LME) model was utilized, using the lmer function in the lmerTest and lmer4 packages implemented in R, to examine the association between reaction times and sine/cosinetransformed phase angles, accounting for potential participantspecific variability. The model included reaction time as the dependent variable, with sine and cosinetransformed phase angles and type of trial (congruent/incongruent) as fixed 2051817x, 2024, 19, Downloaded from https://physoc.onlinelibrary.wiley.com/doi/10.14814/phy2.70079 by Readcube (Labtiva Inc.), Wiley Online Library on [21/10/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 7 of 16 MUÑOZCARACUEL etal. effects. Additionally, a random intercept was incorporated for each participant to account for individual differences in reaction time. In all cases, the p values were adjusted for multiple comparisons when necessary by applying the False Discovery Rate (FDR) correction using the Benjamini- Hochberg procedure. 3 | RESULTS 3.1 | Behavioral task In the active condition task, the average hit rate was 236 out of 240 (98%) ± 5.54 SD, with half of the sample committing either just one error or none. Therefore, further analysis regarding accuracy performance relationship with phase angle could not be conducted. The mean reaction time was 692 ms (667 ms for congruent trials; 764 ms for incongruent trials) ± 202 ms SD. Differences in reaction time between type of trials reached significant results [t (1, 31) = 7.27; p < 0.001]. 3.2 | Power spectral density The PSD for the different data signals around the frequency of interest (0.33 Hz) is shown in Figure2a. The repeatedmeasures ANOVA found significant statistical differences between conditions for all the recorded signals (Table1). The additional repeatedmeasures ANOVA, incorporating channels as a factor, revealed a significant main effect of channels on the PSD around 0.33 Hz for fNIRS- e [F (6.029, 186.902) = 5.813; p < 0.001; η2 = 0.158], fNIRS- i [F (6.644, 205.964) = 7.462; p < 0.001; η2 = 0.194], and EEG [F (2.977, 92.291) = 5.415; p = 0.002; η2 = 0.149]. These results indicate that PSD varies across different fNIRS- i, fNIRS- e, and EEG channels. A significant interaction between condition and channels were also found in the case of EEG [F (3.124, 96.843) = 3.614; p = 0.015; η2 = 0.104]. However, this interaction and the effect of channels in EEG and fNIRS- e will not be further examined in this study, as our focus is on fNIRS- i for localization topography due to its superior spatial resolution. A cortical topographical visualization of the PSD of fNIRS- i results around 0.33 Hz is depicted in Figure3a, which illustrates paired ttests (FDR corrected for 126 comparisons: 3 conditions × 42 channels) between conditions across all fNIRS- i channels. 3.3 | Phase coherence The coherence of the different data signals with the auditory stimulation signal (0.33 Hz) is shown in Figure2b. The repeatedmeasures ANOVA showed a main effect of condition on phase coherence with the auditory stimulation signal (0.33 Hz) for all the recorded signals (Table2). Additional repeatedmeasures ANOVA with channels as an added factor found a main effect of channels on phase coherence with the auditory stimulation signal (0.33 Hz) for fNIRS- i [F (13.698, 397.236) = 1.895; p < 0.001; η2 = 0.061] and EEG [F (2.800, 81.214) = 6.923; p < 0.001; η2 = 0.193]. Significant interactions between conditions and channels were also found in fNIRS- i [F (16.654, 482.971) = 2.070; p = 0.008; η2 = 0.067] and EEG [F (4.125, 120.204) = 4.429; p = 0.002; η2 = 0.132]. These results indicate that coherence with the auditory stimulation signal (0.33 Hz) varies across fNIRS- i and EEG channels, but only under specific conditions. Figure3b displays paired ttests (FDR corrected for 126 comparisons: 3 conditions × 42 channels) comparing conditions for the coherence between the HbO signal of all fNIRS- i channels and the auditory stimulation signal at 0.33 Hz, providing a topographical visualization of the coherence results on the cortical surface. 3.4 | Phase angle The circular representations of the phase angles in the active condition for each processed signal and type of driver (S1, S4, and R) are shown in Figure4. Rayleigh's test revealed significant unimodal clustering among participants (nonuniform circular distribution) around a specific mean angle direction (θ) for the signals specified in Table3 (FDR correction applied to 24 comparisons: 3 types of drivers × 8 signals). These results indicate that for certain signals, specific phase angles were consistently prioritized across participants at the onsets of stimuli and responses. A mixedeffects model was applied to all recorded signals. FDR correction was applied to 48 comparisons, resulting from 3 types of drivers × 8 signals × 2 transformed phase angles (sine and cosine). Type of trial showed a clear association with reaction time, with longer reaction times for incongruent trials (t = 17.489, p (FDR) < 0.001). Significant associations were also identified between reaction time and the sinetransformed phase angle of the response (R) in RSP (t = 2.864, p (FDR) = 0.038; shorter RT when closer to 270°), the cosinetransformed phase angle of S1 in HR (t = −2.654, p (FDR) = 0.048; shorter RT when closer to 0°), the sinetransformed phase angle of S4 in HR (t = −3.744, p (FDR) = 0.003; shorter RT when closer to 90°), the cosinetransformed phase angle of S1 in EEG (t = 3.042, p (FDR) = 0.028; shorter RT when closer to 180°), the sine and cosine transformed phase angle of S4 in EEG (t = 2.620, 4.716; p (FDR) = 0.040, <0.001, 2051817x, 2024, 19, Downloaded from https://physoc.onlinelibrary.wiley.com/doi/10.14814/phy2.70079 by Readcube (Labtiva Inc.), Wiley Online Library on [21/10/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 8 of 16 | MUÑOZCARACUEL etal. respectively; shorter RT when closer to 180° and 270°, respectively), the cosinetransformed phase angle of the response in EEG (t = 2.643, p (FDR) = 0.040; shorter RT when closer to 180°) and the sinetransformed phase angle of S1 in fNIRS- i (t = −2.663, p (FDR) = 0.040; shorter RT when closer to 90°). Random effects for participant indicated significant variability in the intercept of reaction time (participant: Var = 11,619, SD = 108, in milliseconds). Overall, these results indicate that longer reaction times were observed when the RSP phase angle was near 90° at response onset, approximating the mean phase angle of 63° across participants. Similarly, longer reaction times occurred when the fNIRS- i phase angle was near 270° at S1 onset, close to the group mean of 301°. In contrast, shorter reaction times were associated with an EEG phase angle around 180° at S4 onset, which approximates the group mean phase angle of 164°. Additionally, shorter reaction times were noted when the HR phase angle was near 90° at S4 onset, although this did not align with the mean phase angle of 347° across participants. Figure5 displays the average RT computed for all trials and participants at every phase angle (360°). Only signals that show significant associations between reaction time and sine/cosinetransformed phase angles for S1, S4, or R are displayed. 4 | DISCUSSION This study aimed to shed light on the complex dynamics of brain and body rhythms in response to rhythmic auditory stimulation and their behavioral implications, distinguishing between passive and active sensory encoding FIGURE 2 (a) PSD for the three experimental conditions (baseline, passive, and active) and the different data signals. The upper bounded line represents the accounted frequency window for statistical analysis. (b) Violin plots of the phase coherence between the different data signals and the auditory stimulation at 0.33 Hz for the three experimental conditions. For EEG, fNIRS- i, and fNIRS- e, the plotted PSD and phase coherence was the average of all channels. Vertical white lines inside the violins indicate the interquartile range, whereas the white point indicates the median value. Each colored point inside the violins represents a participant (n = 30). Significant differences between conditions are specified with an asterisk for p (FDR) < 0.05 and two asterisks for p (FDR) < 0.01. A, active condition; B, baseline condition; EEG, electroencephalographic activity; fNIRS- e, extracerebral HbO; fNIRS- i, intracerebral/cortical HbO; HR, heart rate; P, passive condition; PPG, peripheral pulse; PTT, pulse transit time; RSP, respiratory effort; SCL, skin conductane level. 2051817x, 2024, 19, Downloaded from https://physoc.onlinelibrary.wiley.com/doi/10.14814/phy2.70079 by Readcube (Labtiva Inc.), Wiley Online Library on [21/10/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 9 of 16 MUÑOZCARACUEL etal. processes. Our findings indicate that active engagement with rhythmic auditory stimuli holistically entrains various bodily oscillations to the external rhythm, offering insights into the intricate interplay between bodily and environmental rhythms and their implications for sensorimotor performance. 4.1 | Power spectral density The PSD results revealed a pronounced peak at 0.33 Hz across all central and autonomic nervous system signals recorded during the active condition, corresponding to the frequency of auditory stimulation. Statistical analysis confirmed significant differences in PSD between the active condition and both passive and baseline conditions, except for PPG, which did not exhibit significant differences between the active and passive conditions. These results align with existing literature on neural entrainment as a mechanism of selective attention (Lakatos etal.,2008, 2019; Obleser & Kayser, 2019), extending this phenomenon to cortical oxygenation levels and autonomic responses during attentional engagement in sensorimotor cognitive tasks. Additionally, the PSD of the EEG at 0.33 Hz was higher in the passive condition compared to the baseline condition, indicating neural entrainment to sensory stimulation even in the absence of conscious attentional effort, consistent with previous reports (Criscuolo, Schwartze, Henry, et al., 2023; Criscuolo, Schwartze, Prado, etal.,2023). However, no entrainment was observed for cortical oxygenation levels during the passive condition compared to the baseline condition, suggesting that passive entrainment may require lower metabolic brain resources. Concerning the cortical topography of the PSD results at the frequency of auditory stimulation, measured through fNIRS- i, the findings demonstrated a greater amplitude of cortical HbO at 0.33 Hz in the active condition compared to both passive and baseline conditions, particularly in the inferior frontal gyrus (IFG), superior temporal gyrus (STG), and inferior parietal lobule (IPL), which are cortical regions associated with keeping the meter in musical contexts (Vuust etal.,2006). These PSD results suggest a generalized entrainment across various bodily oscillations as a means of actively engaging with rhythmic auditory stimulation, while the entrainment of cortical oxygenation levels exhibits a localized effect in brain regions specialized for auditory processing. 4.2 | Phase coherence Consistent with the PSD results, the coherence analysis provided additional evidence of increased entrainment to auditory stimulation for all recorded signals during the active condition, in comparison to both passive and baseline conditions. In this case, alongside EEG, SCR also demonstrated heightened synchronization with incoming auditory inputs during the passive condition compared to the baseline. Interestingly, no significant differences were observed in SCR coherence with auditory inputs between the active and passive conditions. Given that phasic EEGEDA responses have been shown to represent saliencedependent reactions to incoming stimuli as part of the orienting or defensive reflex, irrespective of attentional engagement (Muñoz- Caracuel etal.,2021; Sokolov,1990), these findings could be interpreted as indicative of a subtle form of underlying communication between the CNS and ANS occurring in the context of rhythmic stimulation with variable sequence presentation. The cortical topography of coherence results with auditory stimulation, in contrast to PSD cortical topography, revealed greater coherence of HbO levels with the auditory stimulation across the entire recorded cortical surface during the active condition, rather than at more specialized TABLE 1 Statistical results from the PSD analysis around the frequency of interest (0.33 Hz). Signal df F p η2 SCR 1.236, 38.309 7.418 0.006 0.193 HR 1.400, 43.403 10.068 0.001 0.245 PPG 1.403, 43.487 3.586 0.051 0.104 RSP 1.279, 39.639 11.636 < 0.001 0.273 PTT 1.023, 31.710 7.659 0.009 0.198 fNIRS- e 1.426, 44.195 7.601 0.004 0.197 fNIRS- i 1.362, 42.234 14.172 < 0.001 0.314 EEG 1.088, 33.741 39.592 < 0.001 0.561 Note: Post hoc comparison analysis (FDR corrected for 24 comparisons: 3 conditions × 8 signals) revealed that the PSD of SCR, HR, RSP, EEG, fNIRS- i, fNIRS- e, and PTT around 0.33 Hz were higher for the active as compared with the passive condition [t (1, 29) = 2.546, 3.289, 3.149, 6.544, 3.509, 3.469, and 2.645, respectively; p (FDR) = 0.016, 0.003, 0.004, <0.001, 0.001, 0.002, and 0.013, respectively], as well as for the active as compared with the baseline condition [t (1, 29) = 3.093, 3.430, 3.872, 6.759, 4.308, 2.422, and 2.890, respectively; p (FDR) = 0.004, 0.002, <0.001, <0.001, <0.001, 0.021, and 0.007, respectively]. The PSD of PPG was greater for the active as compared with the baseline condition [t (1, 29) = 2.056; p (FDR) = 0.048]. The PSD of EEG and PTT was also higher for the passive as compared with the baseline condition [t (1, 29) = 3.333, 2.218, respectively; p (FDR) = 0.002, 0.034, respectively]. Altogether, these results indicate a significant amplitude peak at 0.33 Hz across all recorded signals during the active condition. Abbreviations: df, degrees of freedom; EEG, electroencephalographic activity; F, F- statistic; fNIRS- e, Extracerebral HbO; fNIRS- i, Intracerebral/ cortical HbO; HR, heart rate; p, pvalue; PPG, peripheral pulse; PTT, pulse transit time; RSP, respiratory effort; SCR, skin conductance response; η2, eta squared (effect size). 2051817x, 2024, 19, Downloaded from https://physoc.onlinelibrary.wiley.com/doi/10.14814/phy2.70079 by Readcube (Labtiva Inc.), Wiley Online Library on [21/10/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 16 of 16 | MUÑOZCARACUEL etal. Pezzulo, G., Rigoli, F., & Friston, K. (2015). Active inference, homeostatic regulation and adaptive behavioural control. Progress in Neurobiology, 134, 17–35. https:// doi. org/ 10. 1016/j. pneur obio. 2015. 09. 001 Pikovsky, A., Rosenblum, M., & Kurths, J. (2001). Synchronization: A universal concept in nonlinear sciences. Cambridge University Press. https:// doi. org/ 10. 1017/ CBO97 80511 755743 Rapp, P. E. (1987). Why are so many biological systems periodic? Progress in Neurobiology, 29, 261–273. https:// doi. org/ 10. 1016/ 0301- 0082(87) 90023 - 2 Ruiz- Martínez, F. J., Morales- Ortiz, M., & Gómez, C. M. (2022). Late N1 and postimperative negative variation analysis depending on the previous trial history in paradigms of increasing auditory complexity. Journal of Neurophysiology, 127(5), 1240–1252. https:// doi. org/ 10. 1152/ JN. 00313. 2021 Scholkmann, F., Spichtig, S., Muehlemann, T., & Wolf, M. (2010). How to detect and reduce movement artifacts in nearinfrared imaging using moving standard deviation and spline interpolation. Physiological Measurement, 31(5), 649–662. https:// doi. org/ 10. 1088/ 0967- 3334/ 31/5/ 004 Schroeder, C. E., & Lakatos, P. (2009). Lowfrequency neuronal oscillations as instruments of sensory selection. Trends in Neurosciences, 32(1), 9–18. https:// doi. org/ 10. 1016/j. tins. 2008. 09. 012 Seth, A. K., & Friston, K. J. (2016). Active interoceptive inference and the emotional brain. Philosophical Transactions of the Royal Society, B: Biological Sciences, 371(1708), 20160007. https:// doi. org/ 10. 1098/ RSTB. 2016. 0007 Sokolov, E. N. (1990). The orienting response, and future directions of its development. The Pavlovian Journal of Biological Science, 25(3), 142–150. https:// doi. org/ 10. 1007/ BF029 74268 Stefanics, G., Hangya, B., Hernádi, I., Winkler, I., Lakatos, P., & Ulbert, I. (2010). Phase entrainment of human delta oscillations can mediate the effects of expectation on reaction speed. Journal of Neuroscience, 30(41), 13578–13585. https:// doi. org/ 10. 1523/ JNEUR OSCI. 0703- 10. 2010 Vuust, P., Roepstorff, A., Wallentin, M., Mouridsen, K., & Østergaard, L. (2006). It don't mean a thing… Keeping the rhythm during polyrhythmic tension, activates language areas (BA47). NeuroImage, 31(2), 832–841. https:// doi. org/ 10. 1016/J. NEURO IMAGE. 2005. 12. 037 Vuust, P., & Witek, M. A. G. (2014). Rhythmic complexity and predictive coding: A novel approach to modeling rhythm and meter perception in music. Frontiers in Psychology, 5, 1–14. https:// doi. org/ 10. 3389/ fpsyg. 2014. 01111 Wilke, J. T., Lansing, R. W., & Rogers, C. A. (1975). Entrainment of respiration to repetitive finger tapping. Physiological Psychology, 3(4), 345–349. https:// doi. org/ 10. 3758/ BF033 26838/ METRICS Xiong, L., & Garfinkel, A. (2023). Are physiological oscillations physiological? Journal of Physiology, Advance online publication. https:// doi. org/ 10. 1113/ JP285015 Zelano, C., Jiang, H., Zhou, G., Arora, N., Schuele, S., Rosenow, J., & Gottfried, J. A. (2016). Nasal respiration entrains human limbic oscillations and modulates cognitive function. Journal of Neuroscience, 36(49), 12448–12467. https:// doi. org/ 10. 1523/ JNEUR OSCI. 2586- 16. 2016 Zion Golumbic, E. M., Ding, N., Bickel, S., Lakatos, P., Schevon, C. A., McKhann, G. M., Goodman, R. R., Emerson, R., Mehta, A. D., Simon, J. Z., Poeppel, D., & Schroeder, C. E. (2013). Mechanisms underlying selective neuronal tracking of attended speech at a “cocktail party”. Neuron, 77(5), 980–991. https:// doi. org/ 10. 1016/J. NEURON. 2012. 12. 037 How to cite this article: Muñoz- Caracuel, M., Muñoz, V., Ruiz- Martínez, F. J., Vázquez Morejón, A. J., & Gómez, C. M. (2024). Systemic neurophysiological entrainment to behaviorally relevant rhythmic stimuli. Physiological Reports, 12, e70079. https://doi.org/10.14814/phy2.70079 2051817x, 2024, 19, Downloaded from https://physoc.onlinelibrary.wiley.com/doi/10.14814/phy2.70079 by Readcube (Labtiva Inc.), Wiley Online Library on [21/10/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License