What Does Sleeping Brain Tell About Stress? A Pilot Functional Near-Infrared Spectroscopy Study Into Stress-Related Cortical Hemodynamic Features During Sleep

1 Introduction

There is abundant evidence that mental stress is often linked to reduced sleep quality, suggesting abnormalities in brain activity during sleep when people are stressed (Buysse et al., 2011). While our understanding into how stress affects brain activity when we are awake (and are engaged in lab-based stress induction tasks) has been greatly advanced in recent years (Alonso et al., 2015; Kramer et al., 2017; Chang and Yu, 2018; Rosenbaum et al., 2018; Rampino et al., 2019; Schaal et al., 2019; Rosenbaum et al., 2021), no study has looked at how stress modulates brain activity during sleep. Attempts to study stress during sleep face several challenges. First, traditional neuroimaging techniques for studying stress in daytime are not suited for in-sleep measurement due to various methodological restraints imposed by these techniques. Hemodynamic imaging methods such as functional magnetic resonance imaging (fMRI) can generate hemodynamic profiles at high spatial resolution, but they are invasive as people could hardly fall asleep in noisy fMRI scanners. Functional near-infrared spectroscopy (fNIRS) achieves better trade-off between convenience and spatial resolution, but traditional fNIRS systems still use many cables which make them unsuited for measurement during sleep. Electrophysiological neuroimaging techniques such as electroencephalography (EEG) is widely used to measure brain activity during sleep, but their spatial resolution is limited, and they are sensitive to motion artifacts. Second, laboratory-induced stress response is often temporary, and the effect could barely sustain until and throughout nocturnal sleep (Rosenbaum et al., 2021). Established methods for inducing social stress (e.g., the Trier Social Stress Test (Chang and Yu, 2018)), emotional stress (e.g., viewing scary pictures (Rampino et al., 2019)), and physical stress (e.g., sleep deprivation (Alonso et al., 2015)) may fail to mirror natural stress responses, as a laboratory setting often does not represent the typical conditions under which stress occurs in real life (Wolfram et al., 2013). Laboratory-induced stress responses often fade out in an hour (Rosenbaum et al., 2021; Rosenbaum et al., 2018), while real-life stress responses could last hours to days or even longer after the onset of the stressors (Joëls and Baram, 2009). Another pitfall of lab-based stress induction protocols is that they are unsuited for longitudinal repeated measurement from individual subjects as they are likely to cause response habituation especially in the hypothalamic–pituitary–adrenal (HPA) axis as indicated by the cortisol secretion level (Schommer et al., 2003; Kudielka et al., 2006; Jönsson et al., 2010; Gianferante et al., 2014). In addition, stress has been routinely treated as a dichotomous variable (i.e., stress is either present or absent) in many research studies. In real life, however, people may experience various levels of stress with different temporal profiles (Joëls and Baram, 2009). The dichotomous perspective of stress is also unnatural when biomarkers of stress responses such as cortisol is used as an indicator because it is difficult to set a universal cutoff line that accommodates interpersonal variability.

This pilot study is the first to look at how bedtime stress associates to brain activity during sleep. We aimed to explore which cortical areas demonstrate stress-related blood flow patterns during the first sleep cycle. Especially, we focused on answering the following two research questions:

• What hemodynamic features are significantly associated to each stress indicator?

• Which sub-regions in the PFC are significantly associated to each stress indicator?

The study design included addressing the limitations of the existing research paradigm. Table 1 highlights the originality of the present study in comparison with previous studies. This study adopted the N-of-1 approach which is an idiographic research methodology that overcomes the pitfalls of the widely adopted large-sample approach. The large-sample approach requires stringent conditions such as cohort homogeneity—a condition difficult to meet no matter how large the sample size is. When the within-subject variability is much larger than the inter-subject variability, which is common in psychology and physiology studies, the large-sample approach often fails to provide us with findings that generalize well to individuals (Molenaar, 2004; Barlow and Nock, 2009; Mehl and Conner, 2012; van Ockenburg et al., 2015; Burg et al., 2017; Fishera et al., 2018; Piccirillo et al., 2019). In contrast, the N-of-1 approach embraces longitudinal repeated measurement on a single subject to generate the most relevant and reliable information for the specific person, which represents a true scientific undertaking (Barlow and Nock, 2009).

TABLE 1. Research paradigm comparison between the current study and previous studies.

With respect to the experiment settings, we performed the measurement at the subject’s home using a cutting-edge wearable fNIRS system together with non-invasive wearable and mobile devices to achieve the highest level of ecological validity. We also did not rely on lab-based stress induction protocols, as they require the subjects to be actively engaged in cognitive tasks and often fail to induce stress responses that sustain until bedtime. Instead, we tracked the stress indicators on a daily basis, and utilized the days of lower levels of measured stress as natural controls to the days of higher levels of measured stress. Stress responses in human may manifest in multiple physiological systems with varied temporal profile (Joëls and Baram, 2009). In this study, stress was quantified using both objective and subjective indicators. The objective indicators included two widely used stress biomarkers that reflect the hormonal and immunological responses to stress: salivary cortisol and secretory immunoglobulin A (sIgA). The rise of salivary cortisol reveals the stress-related changes in the hypothalamic–pituitary–adrenal (HPA) axis. Meanwhile, stress-associated immunological response could occur more rapidly compared to the HPA axis, characterized by a quick and temporal rise and then decrease in sIgA (Engeland et al., 2016). sIgA may also be a valuable indicator for differentiating between positive and negative stress effects or between successful and unsuccessful adaptation or coping with situational demands (Zeier et al., 1996). The subjective perception of stress can be measured using psychometric instruments ranging from as simple as a Likert scale to as complex as the 30-item Perceived Stress Questionnaire (PSQ) (Levenstein et al., 1993). In this study, the perceived stress level was rated on a 1–10 Likert scale that was implemented using a mobile application. Another significant difference of this study is that the measurement was mostly performed when the subject was sleeping, while in previous stress neuroimaging studies the subjects were all awake and were engaged in cognitive tasks.

In this study, stress response was treated as a continuous phenomenon in contrast to the traditional dichotomous perspective. The data analysis focused on finding significant associations between cortical hemodynamic features and each individual stress indicator. Previous studies mostly rely on one feature type—the mean of the concentration changes in oxyhemoglobin (ΔO₂Hb) and deoxyhemoglobin (ΔHHb). While this feature type has the merit of easy interpretation, it fails to fully capture the characteristics of the cortical hemodynamic signals. In search for the most useful stress-association features, we derived a wide range of time-domain, frequency-domain, and nonlinear features from the cortical hemodynamic signals. We also proposed an original ensemble feature ranking algorithm that leverages six different statistical and machine learning techniques to generate an average importance score for each feature.

This pilot study does not intend to generate conclusive findings, but rather serves as a proof-of-concept for a new research paradigm that can be implemented to study stress in unexplored settings (e.g., during sleep). Understanding the neurophysiological mechanism that underlies the relation between stress and sleep has the significance of giving hint to the development of brain activity markers of stress, which can be readily measured and monitored using wearable brain imaging technologies. Despite of being a small-scale pilot study, the data collection and analysis protocols are readily applicable to large-scale studies. The observations from this study serve as a foundation for future research to elucidate where the brain processes stress during sleep, based on which new stress indicators or stress coping strategies may be developed.

2 Data Collection

2.1 Measuring Stress

In this study, we quantified stress using both objective (i.e., cortisol and sIgA) and subjective indicators (i.e., perceived stress rating). Salivary cortisol and sIgA were measured using the SOMA Dual Analyte LFD test kits. These kits can be used for real-time measurement in a naturalistic setting. Saliva samples were collected using oral fluid collector (OFC) swabs and were incubated for 15 min in OFC buffers before being read. The participant was instructed not to eat, drink, or brush teeth 30 min prior to providing saliva samples. The calibration range of cortisol and sIgA were 1.25–40 nmol/L and 25–800 μg/ml, respectively (Dunbar et al., 2015). The validity of the SOMA kits has been examined in previous studies (Mitsuishi et al., 2019). The measured salivary cortisol and sIgA data were manually logged in a CSV file. Perceived stress was rated on a 1–10 Likert scale (1 = not stressed at all; 10 = extremely stressed) which was implemented using a mobile application named HealthLog.

2.2 Measuring Prefrontal Hemodynamics

A wearable functional near-infrared spectroscopy (fNIRS) (Brite 24; Artinis Medical Systems Co., Netherlands) was used to measure the concentration changes in oxyhemoglobin (ΔO₂Hb) and deoxyhemoglobin (ΔHHb) in the PFC. The fNIRS is a non-invasive brain imaging technique that strikes a good trade-off between temporal and spatial resolution (Tak and Ye, 2014). The advantage of the Brite 24 system—which weighs only 300 g—is that it permits the monitoring of ΔO₂Hb and ΔHHb without imposing constraints on the posture and movement of the subject, and thus is suited for studying cortical hemodynamics during sleep. In this study, the Brite 24 consists of 10 transmitters (Tx) and 8 receivers (Rx). The Txs take turns to emit light at wavelengths of 760 nm (dominantly absorbed by HHb) and 850 nm (dominantly absorbed by O₂Hb). They were fixed on a soft neoprene head cap, which ensures the alignment of optode placement across different measurements. The optodes were placed at an interoptode distance of 3 cm to achieve the maximum penetration depth of 1.5 cm and were configured into 27 channels as shown in the Template DAQ state at the bottom of Figure 1. All optodes were placed between the FpZ–F3–Cz–F4–FpZ regions in the PFC according to the international 10–20 EEG system. The sampling rate was set to 50 Hz. The Brite 24 device has a battery life of up to 2.5 h when it is used for continuous online measurement. While the battery life could be extended by connecting the device to an external power supply, we decided not to use that strategy out of safety concern for the subject. The PFC was selected as the region of interest in this pilot study because previous studies performed during wake time have shed light on the role of the PFC in responding to acute and chronic stress (Cerqueira et al., 2007; Hains and Arnsten, 2008; Dedovic et al., 2009; Yuen et al., 2009; Arnsten et al., 2015; Nejati et al., 2021).

FIGURE 1. Screenshot of the OxySoft.

The Brite 24 system consists of companion software named OxySoft. The software allows the real-time inspection of the signal quality of each channel when the fNIRS device is paired up via Bluetooth connection. Channels with poor quality are marked by red dots, as shown at the bottom of Figure 1. The OxySoft supports several visualization formats of the ΔO₂Hb and ΔHHb signals, including time series plots, 2D heatmap, and 3D heatmap in a glass head. The data recorded by the Brite 24 device were synchronized with the software at regular time interval and were stored in temporal files. When a measurement was stopped, OxySoft processed the temporal files to generate a complete data file that contained raw optical density (OD) data.

2.3 Measuring Complementary Physiological Data

In addition to the Brite 24 system, we also used a Fitbit Sense together with the companion Fitbit app to collect complementary data of sleep, heart rate, and breath rate. These data were utilized in the data preprocessing pipeline to remove physiological artifacts, which is described in detail in the next section. Fitbit is well-suited to this study as it supports the collection of multiple streams of physiological signals without imposing additional burden to the subject. Despite that Fitbit devices may not offer medical-grade measurements, numerous validation studies have demonstrated that Fitbit devices can achieve reasonable accuracy and a better trade-off between accuracy and ecological validity (Menghini et al., 2020; Liang and Chapa-Martell, 2019; Liang and Chapa-Martell, 2018).

2.4 Data Collection Procedure

Grounded on the N-of-1 approach, a longitudinal data collection experiment was conducted with a healthy subject (male, 30 years). The principle of the N-of-1 method allows the exclusion of confounding factors pertaining interpersonal differences in health and physiological conditions. In comparison, the traditional large-sample approach requires stringent conditions such as cohort homogeneity and the findings often do not generalize well to individuals (Molenaar, 2004; Barlow and Nock, 2009; van Ockenburg et al., 2015; Fishera et al., 2018). The large-sample approach becomes especially problematic when the variability within subject is much larger than the variability across subjects (Mehl and Conner, 2012; Fishera et al., 2018). There has been increasing evidence that the large-sample approach may not provide us with information that generalizes well to individuals (Molenaar and Campbell, 2009; Burg et al., 2017; Piccirillo et al., 2019), and hence should not be deemed as more scientific than other approaches that explicitly address within-person variability (Mehl and Conner, 2012). On the other hand, the N-of-1 approach has been well-recognized to provide the highest reliability at the individual level (Molenaar, 2004; Molenaar and Campbell, 2009; Mehl and Conner, 2012; van Ockenburg et al., 2015). The subject was recruited through personal connections. The inclusion criteria were 1) healthy subjects aged 18–65 years without chronic diseases, sleep disorders, and mental disorders, 2) has a smartphone, and 3) understand the contents of the informed consent. This study was approved by the Ethics Committee of the Kyoto University of Advanced Science. Written informed consent was obtained from the subject before the data collection experiment started.

The data collection procedure is illustrated in Figure 2. Saliva samples were collected before bedtime at night. We ensured that saliva sample collection was always done during a fixed time period 22:00–23:00 to control the confounding effect of the circadian hormonal rhythm (Oster et al., 2017). The subject was asked to rate how stressful he felt on the HealthLog app after a saliva sample was collected. The Brite 24 and Fitbit Sense were put on the subject when he was ready for sleep. The Brite 24 head cap was placed symmetrically on the subject’s head, and the Fitbit Sense was worn on the non-dominant wrist. To reset the brain to a common baseline, the subject first went through a wake rest phase where he simply sat quietly for 2 min while staying awake. The wake rest phase was followed immediately by the sleep phase. The Brite 24 was left on until it ran out of battery. The subject was instructed to remove and stop the Brite 24 (simply by pressing the main button) when he needed to go to the restroom early morning or when he woke up, whichever happened first. The subject was asked to synchronize the Fitbit Sense with the companion mobile application after waking up.

FIGURE 2. Data collection procedure.

3 Data Analysis

The objective of the data analysis was to identify the channel-wise features derived from the hemodynamic signals that are significantly associated to stress indicators. We first processed the raw OD signals to yield cleaned high quality ΔO₂Hb and ΔHHb signals, and then derived features from the cleaned signals at each channel. The data analysis pipeline was implemented using Python 3.8.8.

3.1 Data Preprocessing

We exported data from all the devices and instruments for preprocessing. Stress data and Fitbit data were aggregated at a 1-day resolution (i.e., one data point for each day during the experiment period). The perceived stress data were exported from the HealthLog app into a CSV file with a premium account subscription. Fitbit data of sleep, heart rate, and breath rate were exported using a web app that we developed in our previous study (Liang et al., 2016). All these data were then merged by matching date stamps.

The fNIRS data were collected at a high sampling rate of 50 Hz and required more complex preprocessing. The total raw signals measured by the Brite 24 consist of several components, and the ΔO₂Hb and ΔHHb related to neural activity is only a small portion. Noisy components are those related to breath, heartbeat, and movement, which need to be removed (Tak and Ye, 2014). The fNIRS data preprocessing pipeline is illustrated in Figure 3.

1) Export raw OD signals. Using the OxySoft, we exported raw OD signals in EDF format so that they were semi-compatible with the data formats supported by the MNE-NIRS Python library (Luke et al., 2021). Although the OxySoft also allows the export of ΔO₂Hb and ΔHHb signals that have been converted from raw OD signals, unfortunately the Artinis format of these signals is not supported by the MNE-NIRS library at the time of this study. It is worth noting that the MNE-NIRS library read in the EDF data as EEG signals by default; hence, additional processing was needed to convert the signal type to fNIRS after loading EDF files using the MNE-NIRS library.

2) Trim the signals. Since we were only interested in the first sleep cycle, we discarded the signal segments before the sleep start time and after the first sleep cycle. The sleep start time as recorded by the Fitbit Sense was used as the start time (T_s) of the effective data. The end time (T_e) of the effective data were set to T_e = T_s + 90 min as the average sleep cycle of healthy adults is 90 min (Feinberg and Floyd, 1979).

3) Remove channels with poor signal quality. The quality of the OD signals could be compromised by many factors during the measurement. While the OxySoft allows real-time inspection of signal quality, it does not provide computational tools to remove channels with poor signals. In our data preprocessing pipeline, we removed channels with poor signal quality using the Scalp Coupling Index (SCI) method (Pollonini et al., 2014). We first performed channel-wise filtering on the OD signals at both wavelengths using a band-pass filter (0.7–1.5 Hz) to preserve only the heartbeat components. The resulting signals were normalized to balance any difference between their amplitudes. The zero-lag cross-correlation between the resulting signals of the same channel—defined as the SCI—was computed and used as a quantitative measure of the signal-to-noise ratio of the channel. Channels with an SCI-value below 0.75 were regarded as poor channels and were removed from the subsequent analysis.

4) Transform OD to ΔO₂Hb and ΔHHb. The modified Beer–Lambert law (MBLL) was applied to convert the OD signals to ΔO₂Hb and ΔHHb signals (Delpy et al., 1988). As shown in Eq. 1, ${ϵ}_{{\text{O}}_2\text{HB}}({\lambda }_i)$ and ϵ_HHB(λ_i) are the extinction coefficients of O₂Hb and HHb at wavelength of λ_i, respectively. L denotes the interoptode distance. PPF(λ_i) denotes the partial pathlength factor, which represents the sensitivity of the measured optical density to the hemoglobin concentration change in a focal region (Steinbrink et al., 2001). In this study, L and PPF(λ_i) were set to 0.03 and 0.1 (Strangman et al., 2014).

5) Filter out physiological systemic responses. The hemodynamic response due to neural activity has frequency content predominantly below 0.5 Hz (in many cases around 0.1 Hz). The ΔO₂Hb and ΔHHb signals were band-pass filtered to remove cardiac and respiratory noise. According to the Fitbit data, the subject typically had a breath and heart rate between 11–13 bmp and 55–80 bmp, respectively, during sleep. Hence, the cutoff frequency of the band-pass filter was set to 0.02–0.18 Hz.

6) Remove motion artifacts. Although fNIRS is considered more resilient to motion artifacts than EEG, abrupt head motion such as tossing and turning in sleep may still induce spikes that contaminate the true cortical hemodynamic signals. We used the correlation based signal improvement (CBSI) method (Cui et al., 2010) for motion artifacts removal. This method is based on the observation that the ΔO₂Hb and ΔHHb signals, which are typically strongly negatively correlated, will become more positively correlated when contaminated with motion artifacts. Correspondingly, the CBSI method removes motion artifact through recovering the negative correlation between the ΔO₂Hb and ΔHHb signals.

7) Compute epoch-wise average. The effective data of each measurement trial spanned over 90 min (generating 270,000 data points each night). The duration was significantly longer than that of traditional fNIRS studies where a measurement is usually at the scale of several minutes. To efficiently analyze such huge amount of data, we averaged the cleaned ΔO₂Hb and ΔHHb signals epoch-by-epoch at a 30-s interval. Each epoch contains 1,500 data points. This step was compliant with the standard procedure for sleep analysis (Iber et al., 2017). The output time series signals are denoted as {X_n: n = 1, … , N} where N = 180.

FIGURE 3. Pipeline for fNIRS data preprocessing.

3.2 Feature Construction

We derived 36 features from the ΔO₂Hb signal and 36 features from the ΔHHb at each channel. These features fall into three groups. The first group contains 11 time-domain features. These features were directly extracted from the cleaned signals.

• Descriptive statistics: mean $(\bar{X})$, standard deviation (σ), maximum (X_max), and minimum (X_min).

• Skewness (skew): a normalized measure of the asymmetry of the probability distribution of a signal.

• Kurtosis (kurt): a normalized measure of the relative importance of tails versus shoulders in causing dispersion of a signal.

• The 5th-order moment (mmt₅): a measure of the relative importance of tails versus center in causing skew of a signal.

• Mean absolute value (MAV): the average of the absolute value of the signal amplitude.

• Root mean square (RMS): a measure of the average power of a signal.

• Zero crossing (ZC): the number of times the signal changes value from positive to negative and vice versa. It can be interpreted as a measure of the noisiness of a signal.

• Crest factor (CF): an indicator of how extreme the peaks are in a signal.

The second group contains two most typical frequency-domain features. Fast Fourier transform (FFT) was applied to convert {X_n: n = 1, 2, 3, … , N} from time domain to frequency domain to extract the following two features.

• Total power (totalSpec): the sum of the spectral components of a signal.

• Maximal power (maxSpec): the maximum amplitude of the spectral components of a signal.

The third group contains 23 features that characterize the nonlinear characteristics of the cortical hemodynamic signals. Human physiological systems are dynamical systems that often exhibit nonlinear characteristics (Goldberger and West, 1992; Cheffer et al., 2021). Previous studies found that nonlinearities are particularly present in the brain (Toyoda et al., 2008; Ma et al., 2018). To extract nonlinear features, we first used Takens’ time-delay embedding to construct a phase space representation of the system as:

where τ is the time delay and d the embedding dimension. The optimal value of τ and d were decided by minimizing the time-delayed mutual information and by the false nearest neighbors method (Kantz and Schreiber, 2003), respectively. The search range was set to [1, 10] for τ and [2, 6] for d at an increment of 1. The signals were then embedded using the optimal τ_opt and d_opt. The maximal Lyapunov exponent (MLE), Hurst exponent (HE), and correlation dimension (CD) were computed from the embedded ΔO₂Hb and ΔHHb signals. Several nonlinear analysis techniques were also applied to derive features, including recurrence quantitative analysis (RQA), Poincaré plots (PP), and detrended fluctuation analysis (DFA). Different measures of entropy were also calculated.

The RQA computes several quantitative metrics from a recurrence plot (RP). A RP is a visualization of the recurrence behavior of the phase space trajectory $\vec{u}(i)$ of a dynamical system. Each element in the RP is calculated by the following equation:

where $\mathrm{\Theta }:\mathbb{R}\to (\mathrm{0,1})$ is the Heaviside step function, ϵ is a cutoff distance, and ∥▪∥ is the Euclidean norm. In this study, ϵ was set to 0.85. The metrics derived from a RP quantify the recurrence behavior of a dynamic system.

Poincaré plot (PP) is a special type of RP used to quantify self-similarity of a dynamical system. It is a scatter plot of each pair of consecutive data points in a time series signal (technogram), which is often in a shape of ellipse. The minor axis (or width) of the ellipse, denoted as SD₁, reflects the level of short-term instantaneous variability. The major axis (or length) of the ellipse, denoted as SD₂, reflects the long-term variability. PP has been widely used in ECG analysis to help diagnose cardio abnormalities (Hoshi et al., 2013).

The DFA method is often used to quantify the fractal scaling properties and is useful for revealing the statistical self-similarity of a signal (Peng et al., 1994). It has been proven particularly useful in neurology studies (Peng et al., 1994; Hardstone et al., 2012). The DFA first converts a signal to mean-centered cumulative sum. The output signal is then split into epochs, detrended, and the RMS is computed. This process is repeated over a range of epoch sizes n at different scale. A linear trend line is then fit to the log (RMS) − log(n) plot. The slope of the fitted trend line, denoted by α, is called scaling exponent.

The derived nonlinear features are summarized below.

• Optimal delay (τ_opt).

• Optimal embedding dimension (d_opt).

• Maximal Lyapunov exponent (MLE): a measure of separation rate of a signal’s trajectories in the phase space. It indicates the predictability of a dynamic system. A positive maximum Lyapunov exponent is an indicator of the presence of chaos (Eckmann and Ruelle, 1985).

• Hurst exponent (HE): a measure of long-term memory (or long-range dependency) of a signal. A value of HE in the range 0–0.5 indicates long-term negative autocorrelation, while a value in the range 0.5–1.0 indicates long-term positive autocorrelation. A value of 0.5 can indicate a completely uncorrelated signal.

• Correlation dimension (CD): an indicator used to distinguish deterministic chaos from stochastic processes.

• Recurrence rate (RR): the number of black dots in a RP excluding the main diagonal line. It is a measure of the relative density of recurrence points in the entire RP.

• Percent determinism (DET): the fraction of recurrence points that form diagonal lines. It reports the percentage of recurrent points in diagonal structures.

• Maximal diagonal line length (D_max): the length of the single longest line in the diagonal direction within an entire RP. The smaller the D_max, the more divergent the trajectories.

• Average diagonal line length (D_avg): the average time that two segments of the phase space trajectory are close to each other.

• Entropy of diagonal lines lengths (ENT_D): the Shannon entropy of the frequency distribution of the diagonal line lengths. It reflects the complexity of the deterministic structure in the system.

• Laminarity (LAM): the histogram of lengths of vertical lines in a RP. It reports the percentage of recurrent points in vertical structures.

• Trapping time (TT): the average length of the vertical lines. It indicates the mean time the system will abide at a specific state.

• Longest vertical line length (V_max): the maximal length of the vertical lines in the entire RP.

• Entropy of vertical lines lengths (ENT_V): the Shannon entropy of the frequency distribution of the vertical line lengths.

• Standard deviation of the minor axis of a PP (SD₁).

• Standard deviation of the major axis of a PP (SD₂).

• Ratio of SD₁ and SD₂ (SD_ratio); computed as SD₁/SD₂.

• Area of the fitted ellipse (S_e): the area of the ellipse fitted into the PP. It is computed as S_e = π × SD₁ × SD₂.

• Scaling exponent (α): the slope of the fitted trend line in DFA, where each epoch has no overlap.

• Scaling exponent with overlap (α_OL): the slope of the fitted trend line in DFA, where each epoch has 50% overlap.

• Sample entropy (sampEn): a measure of the negative natural logarithm of the probability that if two sets of data points of length m have Euclidean distance D [X_m (n₁), X_m (n₂)] < r (n₁ ≠ n₂) then two sets of data points of m + 1 also have Euclidean distance D [X_m+1 (n₁), X_m+1 (n₂)] < r. In this study, r was set to 0.2σ. A lower value for the sample entropy corresponds to a higher probability indicating more self-similarity and less noise in the signal (Richman and Moorman, 2000).

• Permutation entropy (perEn): a complexity measure that captures the order relations between the values of a signal. Signals with smaller perEn are more regular and deterministic, and those with higher perEn are noisier and more random.

• Statistical complexity (SC): the product of the normalized permutation and a normalized version of the Jensen–Shannon divergence between the ordinal distribution and the uniform distribution (López-Ruiz et al., 1995).

3.3 Feature Ranking

We proposed an original ensemble approach to rank channel-wise hemodynamic features for each stress indicator. As outlined in Figure 4, this approach utilized six feature selection statistical and machine learning techniques to generate an average importance score for each feature $\bar{\zeta }$, calculated the correlation coefficient cor and the corresponding p-value between each feature and a target stress indicator, and performed feature pruning based on the specified criteria. Feature ranking was performed on ΔO₂Hb and ΔHHb features separately. The six feature selection techniques included F-test, mutual information, multivariate linear regression, least absolute shrinkage and selection operator (Lasso) regression, Ridge regression, and recursive feature elimination (RFE).

FIGURE 4. Proposed ensemble method for selecting important hemodynamic features.

The F-test and mutual information (MI) are univariate methods that consider the relationship between each feature and a target stress indicator individually. The F-test method performs a hypothesis testing between a model created by just a constant and another model created by a constant and a hemodynamic feature, and hence reveals the significance of each feature in improving the model. The calculated F-statistic was used as the importance scores of features. While the F-statistic only reflects the linear relationship between a feature and a target stress indicator, the MI method analyzes nonlinear relationships by calculating information gain (Estevez et al., 2009; Ross, 2014). The MI between a hemodynamic feature and a target stress indicator reveals the reduction in uncertainty for the stress indicator given the known value of the feature (Ross, 2014), which was used as the importance score of features. The multivariate linear regression (LR) method fits a linear model with coefficients to minimize the residual between the observed values and the predicted values of the stress indicator. The Lasso regression and Ridge regression methods address the over-fitting problem of linear regression through L1-regularization and L2-regularization, respectively. In addition, the Lasso regression is considered a very useful technique in selecting a strong subset of features as it aggressively produces coefficients of 0 for some features. On the other hand, the Ridge regression is suited for data interpretation because useful features tend to have non-zero coefficients. The α parameters were all set to 0.5 for Lasso and Ridge regression. For the three linear regression methods, the estimated coefficients were used as the importance scores of the features. The RFE method selects features by recursively removing the weakest features until only five features are left. The ranking position of RFE was used as the importance score.

The cortical hemodynamic features were all scaled between [0, 1] before each feature selection technique was performed. For feature i (i ∈ S_feature), the importance score generated by method k, denoted by ${\zeta }_i^k(k\in {\mathbf{S}}_{\text{method}})$, were also scaled between [0, 1]. For a stress indicator $j\left(j\right.$ is either cortisol, sIgA, or perceived stress), the scores of all feature ranking methods for feature i were then averaged to produce an average score (denoted as ${\bar{\zeta }}_{i,j}$) of feature i. The computation of the ${\bar{\zeta }}_{i,j}$ is explained in Eq. 4. We also defined the support of feature i with respect to stress indicator j (denoted as support_i,j) as the number of feature selection methods that yielded an importance score above 0.50 for feature i. Since the $\bar{\zeta }$ only indicates the relative ranking of a feature, we calculated the linear correlation between individual feature and a target stress indicator (denoted as cor_i,j) to better interpret the quantitative relationships. The Pearson’s correlation analysis was performed when cortisol and sIgA were used as the stress indicator, and Spearman’s correlation analysis was performed when perceived stress was used as the stress indicator. The p-values (denoted as p_i,j) were calculated to indicate the significance of the correlation coefficients at a significance level of 0.05. The features that satisfied the following four criteria were selected as important features: 1) ${\bar{\zeta }}_{i,j}>0.50$, 2) support_i,j ≥ 3, and 3) cor_i,j > 0.50, and 4) p_i,j < 0.05. Feature ranking was performed in a channel-wise manner and for ΔO₂Hb and ΔHHb signals separately.

4 Results

4.1 Descriptive Statistics

In total 15 days of data were collected from the subject. As shown in Table 2, the average level of salivary cortisol and sIgA were 3.4 nmol and 295.3 μg/ml, respectively. These values were within the normal ranges for healthy adults (Zeier et al., 1996; Oster et al., 2017). Perceived stress ranges from 2.0 to 8.0 with an average score of 3.8. Only 4 out of the 15 days were rated above 5, indicating that the subject did not perceive constant chronic stress during the data collection experiment. A correlation analysis found no significant linear relationship among the three stress indicators.

TABLE 2. Descriptive statistics of stress indicators.

4.2 PFC Hemodynamic Features Associated to Stress

Channel-wise hemodynamic features associated to stress indicators are summarized in Table 3–8. The aggregated frequency of each feature type is illustrated in Figure 5. Visualization of the channels associated to each stress indicator is provided in Figures 6, 7.

TABLE 3. Channel-wise ΔO₂Hb features associated to salivary cortisol.

TABLE 4. Channel-wise ΔHHb features associated to salivary cortisol.

TABLE 5. Channel-wise ΔO₂Hb features associated to salivary sIgA.

TABLE 6. Channel-wise ΔHHb features associated to salivary sIgA.

TABLE 7. Channel-wise ΔO₂Hb features associated to perceived stress.

TABLE 8. Channel-wise ΔHHb features associated to perceived stress.

FIGURE 5. Frequency of associated feature types summed across all channels.

FIGURE 6. Visualization of channels associated to ΔO₂Hb on the optode template. Stress indicator is (A) cortisol, (B) sIgA, and (C) perceived stress, respectively.

FIGURE 7. Visualization of channels associated to ΔHHb on the optode template. Stress indicator is (A) cortisol, (B) sIgA, and (C) perceived stress, respectively.

Stress was associated to features in both time and frequency domains as well as to nonlinear features. For cortisol, the top three associated cortical hemodynamic features of both the ΔO₂Hb and ΔHHb signals are the $\bar{X}$ of channel 16, the sampEn of channel 26, and the SD_ratio of channel 12. These three features were supported by all six feature selection techniques of the ensemble feature ranking algorithm. Higher cortisol level was strongly associated to increased mean of ΔO₂Hb but decreased mean of ΔHHb at channel 26, strongly associated to increased sample entropy of both the ΔO₂Hb and ΔHHb signals at channel 16, and moderately associated to decreased SD_ratio of both the ΔO₂Hb and ΔHHb signals at channel 12. With all channels counted in, the most frequently referred feature type was the time-domain feature $\bar{X}$. Five feature types were found to be associated only to cortisol but not to the other two stress indicators: mmt₅, sampEn, PE, DET, and maxSpec.

Less channel-wise features were found to associate to sIgA. Only two ΔO₂Hb features (i.e., HE of channel 21 and SD_ratio of channel 26) and one ΔHHb feature (i.e., SD_ratio of channel 26) were supported by all six feature selection technique. A lower sIgA level was moderately associated to higher values of the Hurst exponent of the ΔO₂Hb signal at channel 21 as well as increased SD_ratio of both the ΔO₂Hb and ΔHHb signals at channel 26. The most frequently referred feature type for sIgA were HE, SD_ratio, τ_opt, and ZC. In addition, ZC was associated only to sIgA but not to the other two stress indicators.

Two time-domain features associated to perceived stress were supported by all six feature selection methods. These two features were also the only features that had non-zero coefficient when the Lasso method was applied. Higher perceived stress was strongly associated to reduced skewness of the ΔO₂Hb signal and increased crest factor of the ΔHHb signal at channel 16. The most frequently referred feature types were HE, MLE, and α_OL. In the meantime, α, α_OL, and X_min were the feature types specific to perceived stress.

Figures 6, 7 demonstrated that channel 3 (optode pair Tx3-Rx1), 16 (Tx6-Rx5), 20 (Tx9-Rx6), 21 (Tx5-Rx7), and 26 (Tx9-Rx8) were the most relevant channels. The features of both the ΔO₂Hb and the ΔHHb signals at channel 16 were strongly associated to all three stress indicators. Channel 20 had similar relevance but with only moderate associations. The features of the ΔO₂Hb signals at channels 3 and 21 were moderately associated to all three stress indicators. In addition, the features of both the ΔO₂Hb and the ΔHHb signals at channel 26 were moderately associated to both cortisol and sIgA.

5 Discussion

This pilot study demonstrated the feasibility of investigating stressed brain during sleep by configuring a digital ecosystem with wearable/portable devices and mobile applications. The sleep data collected with Fitbit Sense provided information on sleep start time, and the heart and breath rate data facilitated personalized filtering of the fNIRS signals. The use of mobile applications such as the HealthLog app can also help reduce the burden of manual log, and thus is likely to improve subjects’ adherence to the study protocol. We do not intend to draw conclusions due to the pilot nature of the study; instead, we discuss several observations that may inspire future studies in the same direction.

Stress can be measured along multiple dimensions using several indicators. In this study, we measured salivary cortisol, sIgA, and collected subjective ratings on perceived stress. The correlation analysis revealed that while some features derived from the cortical hemodynamic signals associated with all stress indicators, others may be specific to only one or two stress indicators. The analysis showed that the brain activity may be characterized using various features derived from the hemodynamic signals. Time-domain features, frequency-domain features, and nonlinear features all showed promise. Taken together, the top five most frequently referred feature types were Hurst exponent, mean, the ratio of the major/minor axis standard deviation of a Poincaré plot, statistical complexity, and crest factor. Breaking down into individual stress indicator, the mean of the cortical hemodynamic signals is the most frequently referred feature type for cortisol. This coincides with the fact that most studies that rely on cortisol as the stress indicator solely characterize the cortical hemodynamic patterns using the mean of the ΔO₂Hb and the ΔHHb signals. On the other hand, nonlinear features of the hemodynamic signals could be more useful when sIgA and perceived stress are used as stress indicators. Four feature types had the same highest frequency for sIgA: Hurst exponent, major/minor axis standard deviation ratio of the Poincaré plot, optimal delay, and zero crossing. For perceived stress, Hurst exponent, maximal Lyapunov exponent, and over lapped α in DFA were the most frequently referred feature types. It is also found that the time-domain features (e.g., mean) derived from the ΔO₂Hb signals and those from the ΔHHb signals demonstrated opposite correlation directions to stress indicators, whereas the nonlinear features (e.g., Hurst exponent, correlation dimension, and statistical complexity) derived from the two cortical hemodynamic signals demonstrate the same correlation direction. While the two hemodynamic signals share some common important features, each one also contributed unique features. This suggests the necessity of consider both signals when investigating braining activity using fNIRS, which provides support to the argument made in previous studies (Tachtsidis and Scholkmann, 2016). It is also worth mentioning that the time-domain and frequency-domain features have the merit of their interpretability, whereas some nonlinear features such as the Hurst exponent and the optimal delay may hinder straightforward interpretation.

This study also generated preliminary observations on the sub-regions in the PFC associated to stress. The left rostral prefrontal cortex (channel 16; optode pair Tx6-Rx5), or RLPFC for short, is undoubtedly the most relevant sub-region. Significantly strong correlations were found between the hemodynamic features derived at this sub-region and all three stress indicators. The specific features selected at channel 16 varied depending on the target stress indicator. Higher cortisol level (indicating stronger stress response) was associated to increased mean ΔO₂Hb and decreased mean ΔHHb. Lower sIgA level (indicating stronger stress response) was associated to higher levels of noisiness in the ΔO₂Hb and ΔHHb signals characterized by increased zero crossing. Higher perceived stress level was associated to increased symmetry in the ΔO₂Hb signal and higher peak in the ΔHHb signal. The dorsolateral prefrontal cortex (DLPFC) was also a relevant cortical area. The area of mid-DLPFC (channel 3; optode pair Tx3-Rx1) and caudal-DLPFC (channel 20; optode pair Tx9-Rx6) both demonstrated significant and moderate associations to all three stress indicators. In the case of the caudal-DLPFC (channel 20), consistent relationships were found between stress and a same feature: the statistical complexity of the ΔHHb signal. In contrast, a ΔO₂Hb feature at the mid-DLPFC (channel 3) correlates differently to different stress indicators. To be more specific, lower crest factor of the ΔO₂Hb signal at channel 3 was associated to higher cortisol level (indicating stronger hormonal stress response) but at the same time higher sIgA level (indicating weaker immunological stress response). Rather than viewing this as contradictory, an alternative interpretation could be that different stress indicators characterize different facets of the human stress response, suggesting the necessity of using multiple indicators in stress studies. The relevance of the DLPFC in stress response of healthy subjects during wake time has been documented in previous stress studies using near-infrared spectroscopy (NIRS) technique (Yang et al., 2007; Yanagisawa et al., 2011; Rosenbaum et al., 2018; Schaal et al., 2019). Table 9 summarizes the experiment protocols and the main findings of these studies, which were all conducted when subjects were awake and were engaged in cognitive tasks. Our findings provide supplementary support to the within-person role of the DLPFC in processing stress during sleep. In addition, the preliminary result shed light on the possible role of the RLPFC, especially the left RLPFC, in processing stress during sleep. On the other hand, no relevance was found in the ventrolateral prefrontal cortex (VLPFC). While (Yanagisawa et al., 2011) found negative association between the activity in the VLPFC and the subjective ratings on social pain (which was induced by a feeling of social isolation), our finding suggests that this region may not be involved in processing daily life stress during sleep.

TABLE 9. Comparison of Experiment Protocols and Main Findings of Previous Daytime fNIRS Stress Studies and the Current In-sleep fNIRS Stress Study.

The preliminary findings should be interpreted with caution due to the strong limitations of the present study. First, while the idiographic N-of-1 approach adopted in this study represents a true scientific undertaking (Barlow and Nock, 2009), the preliminary findings solely hold for this specific subject. Future studies may conduct the longitudinal measurement on more subjects to identify possible common patterns across subjects. Second, the data analysis protocol in this study did not count in the interplay among different channels nor the confounding effect of different sleep stages. Analyzing the orchestration of the cortical hemodynamic signals from a dynamic network perspective may leads to new insights into how the brain responds to stress during sleep. Furthermore, this study only focused on the activity in the PFC area in the first sleep cycle. The potential role of other cortical areas in stress response during a full course of sleep demands further studies.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics Statement

The studies involving human participants were reviewed and approved by the Ethics Committee of Kyoto University of Advanced Science. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

ZL conceived and designed the study; ZL performed the data collection experiments, retrieved the data, and performed data preprocessing and analysis; ZL drafted the manuscript and made revision.

Funding

This work was supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI (Grant Number 16H07469, 19K20141, and 21K17670).

Conflict of Interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References