Nonlinear neural patterns are revealed in high frequency functional near infrared spectroscopy analysis

Functional Near Infrared Spectroscopy (fNIRS) is a useful tool for measuring hemoglobin concentration. Linear theory of the hemodynamic response function supports low frequency analysis ( < 0.2 Hz). However, we hypothesized that nonlinearities, arising from the complex neurovascular interactions sustaining vasomotor tone, may be revealed in higher frequency components of fNIRS signals. To test this hypothesis, we simulated nonlinear hemodynamic models to explore how blood flow autoregulation changes may alter evoked neuro-vascular signals in high frequencies. Next, we analyzed experimental fNIRS data to compare neural representations between fast (0.2 – 0.6 Hz) and slow ( < 0.2 Hz) waves, demonstrating that only nonlinear representations quantified by sample entropy are distinct between these frequency bands. Finally, we performed group-level distance correlation analysis to show that the cortical distribution of activity is independent only in the nonlinear analysis of fast and slow waves. Our study highlights the importance of analyzing nonlinear higher frequency effects seen in fNIRS for a comprehensive analysis of cortical neurovascular activity. Furthermore, it motivates further exploration of the nonlinear dynamics driving regional blood flow and hemoglobin concentrations


Introduction
Functional near-infrared spectroscopy (fNIRS) is a non-invasive method used to measure brain tissue hemodynamics as a proxy for neural activity (Strangman et al., 2002b).Slow hemodynamic response signals in fNIRS oscillations ( < 0.2 Hz) are well-suited for linear (i.e., spectral) analysis, particularly when the time intervals between stimuli are long (Glover, 1999).However, such a conservative cut-off frequency is aimed at avoiding overlap with frequency components in vascular dynamics that are independent of brain function, such as cardiac pulsatility or respiratory waves (Pinti et al., 2019).High frequency signals are often filtered out because they are assumed to correspond to systemic physiological or instrumentation noise (Huppert et al., 2009).However, disentangling fNIRS artifacts from systemic brain-heart and, more in general, brain-body interaction effects may be difficult (Vikner et al., 2021;Candia-Rivera et al., 2022a,b).
Recent studies suggest that high frequency ( > 0.2 Hz) signals in fNIRS may provide additional information about neurovascular activity (Yücel et al., 2021;Santosa et al., 2018;Ghouse et al., 2020).Indeed, the nonlinear nature of the autonomic nervous system (Goldberger et al., 2002;Marmarelis, 2004;Sunagawa et al., 1998;Barbieri et al., 2017) may impact the hemodynamic response in higher frequencies.Recent research by (Ghouse et al., 2020) assessed entropy estimates of fNIRS signals that contained an upper bound of 0.6 Hz and demonstrated complementary areas of activity when compared to neural correlates observed using linear analysis of fNIRS.This finding supports recent suggestions to use cutoff filters as high as 0.5 Hz (Yücel et al., 2021).Some studies even suggest minimal fNIRS preprocessing, avoiding filtering and instead opting for robust statistics (Santosa et al., 2018).These developments underscore the importance of investigating the potential contribution of higher frequency signals in fNIRS analysis to improve on our understanding of neural and neurovascular activity.
More specifically, to investigate the potential contribution of nonlinearities in higher frequency of fNIRS signals, this study assesses differences between neural representations of fNIRS in the traditional slow wave ( < 0.2 Hz) and the proposed fast wave (0.2-0.6 Hz) frequencies using nonlinear information-theoretic methods.We hypothesize that to fully characterize the cognitive phenomena reflected in fNIRS signals, we need to assess its full spectrum with additional nonlinear analysis.This study uses sample entropy (SampEn) (Richman and Moorman, 2000) to assess fNIRS signals irregularity (and so predictability) based on embeddings of the fNIRS signals (Sauer et al., 1991).
First, we simulated plausible modulations of autoregulatory feedback on blood flow control with stochastic dynamics in the hemodynamic model (Friston et al., 2000) to motivate fast wave nonlinearity analysis.Then, using a mental arithmetic paradigm and a control motor imagery paradigm (Berntson et al., 1996), we investigated the additional information that high frequency components may provide in experimental data when a subject is under cognitive stress.We performed a representational similarity analysis (RSA) (Kriegeskorte et al., 2008;Carlin et al., 2011) on data obtained from a public database (Shin et al., 2017), comparing fast and slow wave fNIRS modalities.Our hypothesis was that fast wave activity would add significant complementary regions, particularly in nonlinear analysis, and that the correlations between fast and slow wave fNIRS activity would be lower during mental arithmetic compared to motor imagery due to regional nonlinear interactions.We assessed multivariate correlation at each detector and hypothesized that mental arithmetic topoplots of activity would be less correlated between fast and slow wave fNIRS.

Simulations
For initial validation of plausible effects of autoregulation dynamics on hemoglobin, we first simulate a mechanistic balloon model of hemodynamic activity (Cui et al., 2010;Friston et al., 2000).First, a flow inducing signal (s), i.e. the response of the vascular system to a neural metabolic demand, is linearly described for the sustenance of incoming blood flow (f in ).Particularly: u is the control signal which is the neural signal generating the flow.ϵ is the efficacy by which the neural signal can sustain a flow-inducing signal, i.e. a response that may dilate vessels to increase blood flow into.This is related to the vascular resistance (r), whose dynamics are ṙ = − r 2 s (Friston et al., 2000).τ s is the time constant for this flow-inducing signal and τ f is the time constant describing the effects of autoregulatory feedback from blood flow.Considering τ f comprises information on autoregulatory effects of blood flow (the time constant for returning back to baseline), this is the parameter we later model to assess its nonlinear effects on hemoglobin concentrations.
From the sustaining blood flow generating signal, a so-called "balloon model" describes the dynamics of the volume of a blood vessel (the balloon), and the permeation of hemoglobin inside and outside the vessel (Buxton et al., 1998).Explicitly, the rate of decay of blood volume is related to blood flowing in and blood flowing out of a vessel: (2) α describes the capacity at which a balloon can expel water, having been distended by the surge of inflow.τ o is the time constant which governs the rate of change of the volume, which is similarly intertwined with the rate of change of deoxyhemoglobin in the venous compartment (q): (3) The function E describes the oxygen extraction coefficient, or efficacy of the tissue in extracting oxygen from the incoming blood, while E o is the resting oxygen extraction coefficient.In (Cui et al., 2010), an extension was proposed to relate the blood volume and total hemoglobin concentration as: Then, oxyhemoglobin is merely o = p − q.
We reiterate the observation in eq. ( 1) that the parameter τ f is particularly related to the autoregulatory feedback.Exactly how its value relates to the true autoregulatory changes is little understood, though literature states that responses to autoregulatory changes occur over a period of 1-2 min (or less than 0.02 Hz) (Lemkuil et al., 2013).To maintain this expected periodicity, while allowing random deviations due to uncertain dynamics, we designed a 2nd order stochastic dynamical model whose power spectral density (PSD) on average has a peak at the frequency of hypothesized autoregulatory changes.
The steady state value of τ f is at 0.8 (normalized units), while the steady state of its change is 0. A is a parameter that modulates how fast it returns to steady state, and the diffusion term is random variations outside the potential normal oscillatory behavior of τ f .When there is no autoregulatory activity, it quickly returns to steady state, with A= 10.When there is autoregulatory activity, it is more free to change with A = .1.black Eq. ( 5) was devised to model the realizations of τ f through simple, non-trivial stochastic differential equations.Such a generative model embeds interpretable 2nd order linear dynamics showing emergent nonlinar oscillatory properties in the observed variable.blackWe integrate the stochastic differential equations with the Euler-Maruyama method (Platen and Bruti-Liberati, 2010).
We numerically simulated 100 τ f time courses to illustrate expected spectral properties of the stochastic dynamics.This was achieved by obtaining an estimate of the power spectral density (PSD) using the Welch method (Welch, 1967).The units of the PSD for τ f is in arbitrary units, considering its the time constant of a normalized blood flow signal seen in eq. ( 1).We simulated a block design experiment using these generated τ f time courses, where each block is 40 s long with a minute long rest, and 5 stimulus blocks would either induce τ f modulations (i.e. where A = 0.1) or 5 stimulus blocks would not induce the modulations.The final simulated hemoglobin concentration time courses are corrupted with a signal-to-noise (SNR) ratio of 0 dB.

Experiment data
A publicly available dataset was used to obtain fNIRS signals with the desired experimental protocol for this study, as reported in (Shin et al., 2017).In summary, the experiment recruited twenty-nine healthy subjects (aged 28.5 ± 3.7), fifteen of which were female and fourteen male.Of these twenty-nine healthy subjects, all were right-handed expect for one.As according to Shin et al. (2017), all participants were free from neurological, psychiatric, or any brain-related disorders.They were fully informed about the experimental process, and written consent was obtained from each volunteer.After the experiment, they received financial compensation.The study adhered to the guidelines of the declaration of Helsinki and received approval from the Ethics Committee of the Institute of Psychology and Ergonomics, Technical University of Berlin (approval number: SH_01_20150330).Three trials were performed with ten repetitions of mental arithmetic and baseline events, and three trials were performed with ten repetitions of right-hand and A. Ghouse et al. left-hand motor imagery for each subject.We note that the motor imagery tasks and mental arithmetic tasks were done independently, not concurrently.Thirty-six fNIRS series were acquired for each subject with 10 Hz sampling rate.
The experiment design had sixty seconds of resting state to start data acquisition from a subject, after which an instruction was shown on the screen indicating which task was to be performed-either an arithmetic problem, a " − " for a baseline, or a left or right arrow for motor imagery.The subject performed the indicated task for ten seconds, with a subsequent fifteen second resting phase before the next instruction.After twenty repetitions of these instructions and tasks (ten repetitions for each task in an experiment run), a sixty second rest was performed.A total of three trials were performed, for a total of thirty repetitions per event.

fNIRS signals
Thirty six channels of optical densities (OD) were resolved from source detector pairs comprising 760 nm and 850 nm wavelengths at distances of 3 cm covering the frontal, lateral parietal and posterior cortical regions as seen in Fig. 1a.The modified Beer Lambert law was used to convert the ODs to deoxyhemoglobin (Hb) and oxyhemoglobin (HbO) (Strangman et al., 2002a).
Fig. 2 illustrates the preprocessing pipeline.After applying the modified Beer-Lambert law to resolve the 36 channels seen in Fig. 2b from the source-detector pairings in Fig. 2a, band-pass frequency filters were applied to extract traditional hemodynamic bands ( < 0.2 Hz) (Strangman et al., 2002a;Pinti et al., 2019) or the proposed increased hemodynamic band (0.2-0.6 Hz).A wavelet filtering approach using a Daubechies 5 wavelet, nine level decomposition was used to further reduce instrumentation noise such as movement in the oxy-and deoxyhemoglobin signals (Molavi and Dumont, 2012).Detrending and prewhitening with an AR(1) model was then performed to remove temporally structured noise in the signal (Huppert, 2016).The signals were separated into epochs of 30 s (comprising the 15 s task phase, and the 15 s rest phase), with each channel at each activity block being referenced to the mean of the previous 5 s.We exclude the 2 second instruction phase before the task-phase from the analysis to reduce potential confabulation of task-evoked responses from processing the instructions.Total hemoglobin was computed as the addition of both Hb and HbO.The total hemoglobin is important considering the integration of the concentrations may present unique temporal dynamics revealing distinct nonlinear temporal effects.For all three signals, entropies and mean estimates were extracted then passed into the 1st and 2nd level analysis.Before second level analysis, the results were spatially smoothed to increase the sensitivity of random effects from detector locations (Tak et al., 2016).A Gaussian smoothing kernel with a full-width-half-max of 1.5 cm (half the source-detector separation distance) was used as seen in Fig. 1b.

Entropy analysis
To analyze nonlinearity and regularity of the fNIRS signals, we exploited sample entropy (SampEn) (Richman and Moorman, 2000).SampEn is a method to calculate the entropy of a dynamical system in its phase space.In other words, it assesses how much information it takes to characterize the dynamics of a system.It shares similarities with methods like approximate entropy (ApEn) (Pincus, 1991), except the SampEn algorithm alleviates biases that are introduced in ApEn by not considering self-matches in its calculation of the correlation integral used in the definition of entropy in dynamical systems (Delgado-Bonal and Marshak, 2019b;Sinai, 1959).There are a whole slew of other entropy estimation methods that could have been used instead, however we have previously performed a study comparing a whole battery of entropy estimation methods for assessing fNIRS signals, with the conclusion that SampEn provides similar results as the other entropy estimation methods (Ghouse et al., 2020).
A delay-time τ and embedding dimension m are needed to reconstruct manifolds using delay-coordinates.τ was selected as the first zero of the autocorrelation, while m was found using the false nearest neighbors approach (Abarbanel et al., 1993).
For calculating SampEn, radius R = 0.2 × σ x was used as the threshold to determine whether states were neighbors, where σ x is the standard deviation of the fNIRS time series (Delgado-Bonal and Marshak, 2019a).The particular equation describing how to calculate SampEn is: X in this equation is a state space reconstructed using a time series.The superscript denotes the embedding dimension while the subscript denotes the state index, for which there exists N states.

Representational similarity analysis
Representational similarity analysis is a method used to compare multivariate data such as different brain data types (Kriegeskorte et al., 2008).This paper calculates similarities between tasks using the distance correlation in order to construct the representational similarity matrices for fast or slow wave (Székely et al., 2007;Geerligs et al., 2016).
In brief, given random vector X and random vector Y with dimensions R p and R q respectively, distance correlation evaluates independence by integrating the distance between the random vectors in conjunction with a weighting function w(t, s) = (c p c q ⃒ ⃒ ⃒t| Γ( 1+d2 ) corresponds to half the surface area of a unit sphere in the given dimensionality d.This leads to the following statistic: As the dimensions of the random vectors, p and q, approach infinity, the statistic converges to a Student's t-distribution and can be approximated as such for hypothesis testing (Székely and Rizzo, 2013).Given a ij = |x i − x j | and b ij = |y i − y j |, where i and j represent the ith and jth observations of x or y, the sample covariance can be estimated as: In our application, random vectors X and Y have equivalent dimensionality R 3 as we are looking at multivariate correlations between oxy, deoxy and total hemoglobin.Furthermore, due to there being 4 experimental conditions (baseline, mental arithmetic, left-and right-hand imagery), the resulting representational dissimilarity matrix (RDM) is a 4 × 4 matrix.Representational dissimilarity outputs were obtained for each detector location for each measure (entropies or mean value).Extracting the upper triangle, we have two 29 matrices for slow and fast wave fNIRS respectively to perform statistical analysis upon.

Fast vs. slow wave spatial analysis
For each subject, the data obtained from each detector had a shape of N repetitions × N concentrations.By calculating the median over the N repetitions for either fast or slow wave fNIRS, we obtain a vector of size N Concentrations x N detectors for each subject.This analysis enables us to determine whether the spatial fNIRS arrays for fast and slow waves are dependent, or whether the null hypothesis that they are independent can be rejected.We expect that the slow and fast wave analysis will become more independent when there is a task inducing changes in autoregulatory activity.We also perform a Wilcoxon paired analysis at the detector level to compare which medians are significantly different in the group of subjects for each concentration.

Statistical analysis
In order to perform statistical analysis of the data, we exploited the Scipy statistics package in Python (Virtanen et al., 2020).

Fig. 2.
Overview of the analysis pipeline used for each fNIRS signal in the dataset, which includes computation of hemoglobin concentrations, band pass filtering, motion artifact correction, removal of linear trends and serial correlations, total hemoglobin calculation, state space reconstructions using delay-coordinate embeddings, estimation of sample entropy or mean value for each trial, first level analysis using median value over trials, and spatial smoothing to improve sensitivity in group analysis results (Tak et al., 2016).
A. Ghouse et al.To compare RDMs between fast and slow wave fNIRS, a group analysis was conducted using a Friedman test, with a Bonferroni correction to account for multiple comparisons.A family-wise p-value of less than 0.05 was considered significant.Nonparametric Friedman statistical tests were used to avoid assumptions of normality in the dataset when performing intergroup analysis (Friedman, 1937).Then, the Bonferroni correction was applied to be as conservative as possible when correcting for family-wise error rates such that any inferences are robust to Type I errors (Tukey, 1953;Dunn, 1961).
For median analysis, rainclouds were generated from standardized Z transformations to compare the shapes of distributions between fast and slow wave analysis.A Kolmogorov-Smirnov test was applied to compare concentrations analyzed from fast and slow wave analysis, correcting for multiple comparisons.A paired sample Wilcoxon signed rank test was performed on the difference of median value for each concentration between mental arithmetic vs baseline or right-hand vs left-hand motor imagery.Detectors were retained for further analysis if any of the 3 concentrations returned significant with a false alarm rate of α = 0.05∕6 = 0.008.For remaining detectors, distance correlation was performed between fast and slow wave fNIRS.

Results
The first set of results are from simulations of neurovascular signals generated by a hemodynamic model, evoked from events that either induce or do not induce autoregulatory changes.The second set of results are from group statistics of representational similarity analysis between fast wave and slow wave cortical distribution of activity, and median analysis of cortical activity.Standard fNIRS analysis refers to time averaged value of fNIRS activity during a task, while nonlinear fNIRS analysis refers to SampEn during the time duration of the task.

Simulations
Fig. 3 illustrates the results of realizations of the hypothesized autoregulatory feedback changes as according to eq. ( 5).The median peak frequency was at 0.02 Hz, corresponding to periodicity of 50 s, with periodicity ranges from 12 ss to 100 s.Fig. 4 illustrates the effects that these autoregulatory changes may have on the hemodynamics, where HbO is the oxyhemoglobin and HbR is the deoxyhemoglobin.Over the 100 realizations of the block design experiment, the mean entropy of the hemodynamic response in the high frequency (0.2-0.6 Hz) regime without autoregulation modulations was at 0.212 (±0.015) bits compared to 0.163 (±0.011) bits when the modulations occur; through a paired Wilcoxon signed rank test, the SampEns were found to be significantly different (p ⋘ 0.001).On the other hand, the power of this high frequency band without autoregulation modulation was found to be 0.0368 (±0.00441) as compared to 0.0372 (±0.00385) during autoregulation modulation; the power bands were not significantly different according to the Wilcoxon signed rank test (p = 0.46).

Representational similarity analysis
Results from the Friedman analysis comparing the upper triangle of the RDMs constructed using either fast wave or slow wave fNIRS signals according to the methods described in section 2.5 can be seen in Fig. 5, for mean and SampEn respectively.For the mean value, no detector has significantly different representations between slow and fast wave fNIRS, whereas each quadrant of the sensor space on the cortex contained significantly different detectors for SampEn.Generally, SampEn contained higher Friedman test scores than the mean analysis.

Fast vs slow wave spatial analysis
Using the methods described in 2.6, we obtained group level median analysis results for fast and slow wave fNIRS.To reiterate, this entailed calculating the median value of the mean or SampEn signal over the repetition of the task performed by the participant, and then performing a contrast analysis between the effects of mental arithmetic (MA) and baselines (BL), or left-hand motor imagery (LH0 and right-hand motor imagery (RH).A distillation of the spatial results can be seen in the raincloud Fig. 6 to visualize distributions of contrasts between the task comparisons across subjects; Kolmogorov-Smirnov tests were performed to ascertain whether the distributions were significantly different.Supplementary fig.S1 shows the spatial cortical map of contrasts between mental arithmetic (MA) and baseline (BL), or right-(RH) or lefthand (LH) motor imagery tasks.
From the raincloud plots in Fig. 6, distributions for mean results consistently appear to be overlapped when comparing slow and fast wave different in medians, whereas SampEn demonstrates a flatter distribution for slow wave analysis as compared to fast wave analysis.Motor imagery also presents a flatter distribution in SampEn when using slow wave analysis compared to fast wave analysis, however this greater variance is similarly observed in mean value analysis.Performing a Kolmogorov-Smirnov two sample test, we found that the distribution of fast and slow wave analysis was significantly different between SampEn Fig. 3. Spectral analysis of realizations of τ f simulated from the stochastic differential equation in eq. ( 5).a) Represents the average power spectral density (PSD) of τ f as estimated by the Welch method over 100 realizations.The units of the PSD are arbitrary as τ f is the time constant of f in seen in eq. ( 1).b) on the other hand is a histogram of the peak frequency of the PSD of τ f over the 100 realizations.in oxy and deoxyhemoglobin, while in no mean analysis were they significantly different.
For mental arithmetic and baseline Wilcoxon paired analysis, mean estimates show significant detectors only in slow wave fNIRS, with a cluster in the right lateral cortical areas and detectors bilaterally in the frontal cortex.SampEn mainly determined significant detectors in the frontal right cortex and left parietal cortex in the slow wave analysis, while fast wave analysis complemented information in the right parietal cortical areas.
Left-and right-hand motor imagery Wilcoxon paired analysis demonstrated activity bilaterally in the parietal cortical areas with both SampEn and mean estimates.Both mean and SampEn demonstrated slow activity predominantly in the right parietal cortex, while fast activity was found in the left lateral parietal cortex.Mean also contained significant detectors in the medial left parietal cortex in slow wave analysis and in the occipital regions for fast wave analysis.
As a sanity check, we also compared the neural representations between slow frequency bands (0 Hz, 0.2 Hz) and cardiac frequency bands (0.8 Hz, 3 Hz).The distance correlation of topolots of fast wave vs slow wave was not significant, indicating that the fast wave fNIRS analysis in freqs (0.2 Hz, 0.6 Hz) was not a result of systemic physiological confounders (see also supplementary fig.S2).

Discussion and conclusion
In this study, we compared neural activity representations between  fast wave (0.2-0.6 Hz) fNIRS and the standard slow wave ( < 0.2 Hz) fNIRS signal using both linear measures of mean value and nonlinear measures derived from sample entropy.The fNIRS signals were obtained from a publicly available dataset of 29 subjects performing mental arithmetic and left-/right-hand motor imagery tasks (Shin et al., 2017).Our hypothesis was that mental arithmetic tasks, which induce cognitive stress, would provide more variation in areas of cortical activation than motor imagery when comparing slow and fast wave fNIRS using nonlinear analysis.This hypothesis was motivated by literature suggesting the nonlinear modulation of fNIRS dynamics sustained by various factors, including those induced by autonomic nervous system activity (Friston et al., 2000;Gianaros et al., 2012;Sheng and Zhu, 2018;Candia-Rivera et al., 2023).Hemodynamic models are also known to be nonlinear in nature (Friston et al., 2000;Buxton et al., 1998;Friston, 2001).
To investigate the potential effects of vasomotor property dynamics on autoregulation modulations, we first conducted a simulation study using equations ( 1) and ( 5).Our theoretical model for modulations of vasomotor tone closely matched the expected spectral profile of realworld vasomotor tone variations at frequencies below 0.2 Hz, as shown in Fig. 3. Additionally, we found that autoregulation variations in flow signals that propagate hemoglobin concentration dynamics (as depicted in Fig. 4) reduce entropy in the observed signal even with an SNR of 0 dB.This suggests that with proper statistical power, these autoregulatory dynamics could potentially be observed in real-world data using nonlinear analysis methods.It is important to note that our empirical data analysis should not be interpreted as evidence that high frequency signal entropy is a measure of autoregulation.Rather, we hypothesize that the nonlinear activity observed in high frequency fNIRS provides complementary information that is crucial for a comprehensive characterization of cognitive states, such as those associated with stress processing during mental arithmetic tasks.
To explore potential distinctions in neural signal representations between fast and slow waves in fNIRS data, we conducted a Friedman test to compare results on the tasks of baseline, mental arithmetic, and motor imagery, as shown in Fig. 5.The mean value analysis did not reveal any significant differences in task combinations, as expected since the mean corresponds to the DC value of the signal and high pass filtering only affects attenuation of the DC value.In contrast, nonlinar SampEn analysis demonstrated significant differences between fast and slow wave fNIRS detectors in all four quadrants of the sensor space covering cortex, with a higher Friedman score across the cortex than the mean analysis.While it is unclear how filtering affects state space regularity, high pass filtering may result in lower mean entropy and higher variance than low-pass filtering (Borges et al., 2020).This effect may correspond to the significantly different SampEn representational Upon examining motor imagery results we observed changes in laterality in SampEn analysis depending on the frequency band of focus, with fast wave analysis being more sensitive to left hemisphere activity and slow wave analysis being more sensitive to right hemisphere activity (see Supplementary fig.S1).A similar phenomenon was seen in the mean analysis, although an additional medial left hemispheric detector was observed for slow waves.We expect this activity to correspond to the somatosensory cortex, located in the central cortical regions, where the hands of the cortical homunculus are located between the lateral and medial portions of the cortex (Grodd et al., 2001).By considering both frequency bands using SampEn analysis, we may have discovered its bilateral effects.This bilateral difference between left-and right-hand motor imagery was only found using slow wave analysis with just the mean analysis.It is important to mention that the handedness of a participant plays a critical role in the effect seen in motor imagery.Right-handed individuals appear to have a stronger lateralization effect specific to the side of the hand imagined; left-handed individuals on the other hand have bilateral responses (Crotti et al., 2022).Considering the set of participants used in this study were primarily right-handed, the results presented in this study may be specific to this group.Nonetheless, considering we do paired analysis of the same participants performing left-and right-hand motor imagery, the resulting differences of cortical signal effects observed when comparing fast and slow wave analyses appear to indeed be a neurovascular result.This suggests that to fully describe the complementary areas of activity that entropy provides, we must consider the full frequency profile of potential hemodynamic activity.
We also conducted a distance correlation analysis of the difference of median maps and found that all mean maps were significantly correlated, as expected.The cortical representations reflected in median difference maps between tasks should be similar if we scale the signal by merely a DC value, as previously discussed.However, our hypothesis was that mental arithmetic should not be significantly correlated in entropy analysis while motor imagery should be significant.Instead, we found that neither was significant, although mental arithmetic vs. baseline had a lower correlation than the motor imagery distinction.This may imply that the nonlinear irregular effects induced by vasomotor mechanical dynamics are always occurring, although the degree to which they affect the resulting neural signal varies depending on the task.It's known that a significant portion of cerebral vascular resistance comes from vasomotor control in arterioles (Mandeville et al., 1999).The nonlinear transform of the constant dynamics of vasomotor control to hemoglobin concentrations may always be present, yet their scale is modulated by a stress task such as mental arithmetic, making their high frequency contributions to hemoglobin concentrations more irregular in the case of mental arithmetic.We finally remark that the proposed fast wave analysis within (0.2 Hz, 0.6 Hz) is not affected by non-cognitive phenomena associated with oscillations in the (0.8 Hz, 3 Hz) band (Hoshi, 2016).
We recognize that this study may have its limitations.First, we did not perform simultaneous fNIRS and autonomic signals analysis.Respiration and heart rate, typically occurring within the fast-wave band, could be confounding factors in our data analysis since a multitude of regulatory processes may change during cognitive tasks.However, the heterogeneity in the response across the cortex suggests that such oscillations may not systematically affect fNIRS series all over the scalp.Literature suggests possible differential respiratory effects on the scalp, indicating the need for future studies to take actual respiratory and heart rate signals into account (Candia-Rivera et al., 2022c;Dubois et al., 2016).Additionally, signals corresponding to blood pressure in Mayer waves contaminate fNIRS signals within the low frequency bands (Pinti et al., 2015;Kirilina et al., 2013), and nonlinearity evoked in blood pressure variations could similarly reflect in high frequency fNIRS.Hence, future studies should also monitor end-tidal CO2 as partial pressure of CO2 has been shown to affect vascular dynamics (Mas et al., 2000).Despite these limitations may introduce some uncertainty regarding the origin of nonlinear fNIRS dynamics in the fast wave band, our study's outcomes remain robust and reliable.This is further substantiated by the corroborative evidence provided through our comprehensive simulation-based analysis.
In conclusion, our study shows that nonlinear analysis can detect distinct neural representations in fast (0.2-0.6 Hz) waves compared to conventional slow (<0.2 Hz) wave frequencies in fNIRS data.Traditional methods like mean estimates could not reveal such distinct representations.While the origin of neural hemodynamic activity from a behavioral stimulus may be associated with oscillations at frequencies <0.2 Hz, nonlinear neurovascular interactions may generate fNIRS oscillations at higher frequencies.Thus, when assessing effects on fNIRS, comprehensive characterizations should also consider the nonlinear properties of high-frequency band oscillations.

Fig. 1 .
Fig. 1.Illustration of spatial analysis of fNIRS signals.(a) shows the position of the optodes used to derive the 36 fNIRS channels seen in (b).(b) shows an example of the full-width-half-max of the smoothing Gaussian kernel prior to second level group analysis (Tak et al., 2016) in the color green.

Fig. 4 .
Fig. 4. (a) demonstrates a single realization of a band passed oxyhemoglobin power spectrum density when either autoregulatory activity is occurring with a task, or when there is no autoregulatory activity with a task.(b) displays a representative simulation of the signals in the model and the events that influence them (whether a task induces autoregulatory behavior or not).

Fig. 5 .
Fig. 5. Results of a group level Friedman test comparison between the upper triangle of the RDM for fast (0.2-0.6 Hz) and slow wave ( < 0.2 Hz) fNIRS on the group level, where significance indicates at least one element in the upper triangle had a significantly different median between fast and slow wave fNIRS.Channels that are significant are marked by a thick outline on the marker.

Fig. 6 .
Fig. 6.Rainclouds illustrating the standardized Z transformed group level distributions of the absolute value of the difference between medians for mental arithmetic (a) and motor imagery (b) for either fast wave (0.2-0.6 Hz) or slow wave (0-0.2Hz) fNIRS signals.The "*" represents significant difference between the fast and slow wave distribution for the given measure.