Cardiorespiratory Dynamic Response to Mental Stress: A Multivariate Time-Frequency Analysis

Abstract

Mental stress is a growing problem in our society. In order to deal with this, it is important to understand the underlying stress mechanisms. In this study, we aim to determine how the cardiorespiratory interactions are affected by mental arithmetic stress and attention. We conduct cross time-frequency (TF) analyses to assess the cardiorespiratory coupling. In addition, we introduce partial TF spectra to separate variations in the RR interval series that are linearly related to respiration from RR interval variations (RRV) that are not related to respiration. The performance of partial spectra is evaluated in two simulation studies. Time-varying parameters, such as instantaneous powers and frequencies, are derived from the computed spectra. Statistical analysis is carried out continuously in time to evaluate the dynamic response to mental stress and attention. The results show an increased heart and respiratory rate during stress and attention, compared to a resting condition. Also a fast reduction in vagal activity is noted. The partial TF analysis reveals a faster reduction of RRV power related to (3 s) than unrelated to (30 s) respiration, demonstrating that the autonomic response to mental stress is driven by mechanisms characterized by different temporal scales.

1. Introduction

Stress research has gained growing interest in the last decades. Results from the European Working Conditions Survey in 2000 suggested that 28% of all work-related health problems can be attributed to stress. Work absence associated with stress not only brings along highcosts; for example, in The Netherlands, they estimate a cost of 3 billion € per year [1], but it is also linked with serious health problems such as cardiovascular diseases [2, 3] and musculoskeletal disorders [1, 4, 5]. It is therefore important to identify the mechanisms underlying stress, such as physiological mechanisms.

The autonomic nervous system (ANS) dynamically coordinates, among others, cardiovascular variables (such as heart rate and contractility, blood pressure, and peripheral resistance), respiration, and complex interactions between them. The impact of stress on the cardiovascular system has been investigated extensively via the variability of the RR interval variation (RRV) series [6–10]. RRV analysis is widely used to assess the functioning of the ANS. Starting from the RR interval series, several RRV indices that quantify the activity of the ANS have been proposed [11]. In the power spectrum of RRV, a low-frequency (LF) band (0.04–0.15 Hz) and a high-frequency (HF) band (0.15–0.40 Hz) are defined. LF power is linked to both sympathetic and vagal activity, while HF power is only coupled to vagal outflow. All studies reported reduced vagal modulation and RRV and an increased sympathovagal balance during mental stress.

The response of the respiratory system to mental stress has been reported in [12–14], where it was shown that stress is associated with an increased respiratory rate. Vlemincx et al. also reported the effect on respiratory variability [14]; sustained nonstressful attention reduces the total respiratory variability, while mental load increases the total variability.

In this study, we will focus on the combined analysis of the cardiovascular and respiratory systems, which is motivated by the strong interaction between both systems. Respiratory sinus arrhythmia (RSA) is perhaps the best known cardiorespiratory interaction and is the phenomenon that the heart rate increases during inhalation and decreases during exhalation [15]. This influence is reflected in the HF spectral component of RRV, which is therefore often taken as a measure of RSA. Although many studies show that RSA is related to vagal control [16–18], other researches suggest that the magnitude of RSA changes with the depth of breathing (tidal volume) and the respiratory rate, independently of vagal activity [19–22]. Hence, questions arise regarding whether RSA is a true and valid measure of vagal activity. The lack of consensus about its interpretation limits its practicability.

Nonetheless, there is a strong influence of respiration on RRV, and there is a need to conduct a combined analysis during stress. This has been done in a few studies; Pattyn et al. investigated cardiorespiratory reactivity by means of RSA and separate cardiovascular and respiratory parameters and found an increase in heart rate and a decreased RSA during stress [23]. Zhang et al. reported the effects of mental tasks on cardiorespiratory synchronizations and found reduced synchronization epochs during mental arithmetic [24].

In this study we will focus on the changes in cardiorespiratory interactions during stress using time-frequency (TF) analyses. Spectral analysis has the advantage of having a clear link with physiology but has the limitation that it requires stationary signals, which is mostly not the case with physiological signals. The benefit of time domain analyses lies in the dynamic response to stress on cardiorespiratory coupling that is obtained. Time-frequency analyses combine the advantages of time and frequency domain analyses and thus can be used to analyze nonstationary signals and in particular to quantify the temporal variations of their spectral content [25, 26].

In addition, we aim to conduct partial TF analyses where influences of respiration on RRV are separated. Partial spectral analysis is a well-established, model-free, and fully data-driven technique used to modify the spectrum of a signal based on the information carried by another one [27]. A partial TF spectrum represents a modified version of the TF spectrum of signal 𝑥(𝑡), from which the magnitudes of the components linearly related to signal 𝑦(𝑡) have been reduced in a way which is proportional to the coherence between 𝑥(𝑡) and 𝑦(𝑡). This technique was previously used by Badra et al. to assess the partial coherence between RRV and systolic pressure by removing the influence that respiration exerts on both signals [28]. In this study, we will consider one partial spectrum that contains RRV related to respiration and one that is unrelated to respiration. It is important to note that we do not intend to provide a reliable estimate of vagal modulation. We intend to quantify changes in RRV power related and unrelated to respiration. The motivation to conduct this analysis originates from a previous study that showed that spectral features based on variations in the RR interval series unrelated to respiration yield an almost perfect classification (accuracy = 98%) between rest and stress, while traditional RRV analysis resulted in a classification accuracy of only 57% [29]. This result suggests that RRV unrelated to respiration contains important information about stress mechanisms.

In this paper, we will conduct (cross) TF analyses of RR interval series and respiration and assess their coherence during two similar mental stress tasks and a sustained attention task. We hypothesize finding a strong cardiorespiratory response during the first stress task and also a clear, but reduced, coupling during the second task due to habituation. The attention task is nonstressful and is hypothesized to only express a small reaction. As the interaction between RRV and respiration is vagally mediated and stress is linked to a reduced vagal control, the strength of the cardiorespiratory coupling is expected to decrease as a response to stress. The partial TF spectrum related to stress is hypothesized to mainly have HF power, while it is expected that the partial spectrum unrelated to respiration will primarily contain LF power and only little HF power.

2. Material and Methods

2.1. Experimental Setting

The database for this research consists of ECG (sampling frequency 𝑓_𝑠=1000 Hz) and respiration (𝑓_𝑠=50 Hz) measurements of 43 healthy volunteers (age: 18–22 years) that were recorded at the Department of Psychology and Educational Sciences of the KU Leuven (Leuven, Belgium). The respiration was measured using the LifeShirt System (Vivometrics, Inc., Ventura, CA, USA), which estimates the tidal volume, further used as respiratory signal, by means of respiratory inductive plethysmography (RIP) around the ribcage and the abdomen.

During the experimental protocol, the participants were instructed to conduct two types of tasks. The first task was a nonstressful attention task where the participants had to indicate the largest number on a computer using a mouse cursor. During the second task, the students had to perform a mental arithmetic task which induces stress. The whole protocol consists of an attention task (AT) and 2 mental stress tasks (MT1 and MT2), each followed by a recovery period. The order of the tasks was randomized. Prior to any task, a resting period was recorded during which the participants watched a relaxing documentary (RD). Each task and RD had a duration of 6 minutes. For this study, RD, AT, MT1, and MT2 were used. Due to missing data, the recordings of only 40 students were included in the study.

The experiment was approved by the Ethics Committees of the Department of Psychology and Educational Sciences and of the Faculty of Medical Sciences. The study was in accordance with the Declaration of Helsinki (2008).

2.2. Preprocessing

The RR interval series, the signal that contains the time between two heart beats, is composed by detection of the 𝑅 peaks in the ECG using the Pan-Tompkins algorithm [30]. All 𝑅 peak detections are automatically verified using the algorithm described in [31] and afterwards visually inspected. Next, the respiratory signal and the RR interval series are resampled at 4 Hz using cubic spline interpolation. Both signals are high-pass filtered at 0.003 Hz to remove very slow oscillations.

All processing steps of the data are performed in MATLAB R2012a (MathWorks, Natick, MA, USA).

2.3. Cross Time-Frequency Analysis

The cross time-frequency spectrum 𝑆_𝑥𝑦(𝑡,𝑓) of signals 𝑥(𝑡) and 𝑦(𝑡) is estimated using a time-frequency distribution (TFD) [32]: 𝑆_𝑥𝑦𝑡,𝑓=^+∞_−∞Φ(𝜏,𝜈)𝐴_𝑥𝑦(𝜏,𝜈)𝑒^{𝑗2𝜋(𝑡𝜈−𝜏𝑓)}𝐴𝑑𝜈𝑑𝜏,_𝑥𝑦(𝜏,𝜈)=^+∞_−∞𝑥𝜏𝑡+2𝑦^∗𝜏𝑡−2𝑒^{−𝑗2𝜋𝜈𝑡}𝑑𝑡, (1) where 𝐴_𝑥𝑦(𝜏,𝜈) is the cross-ambiguity function. Smoothing is performed by an exponential kernel, in the ambiguity domain, defined as Φ()𝜈𝜏,𝜈=exp−𝜋𝜈₀²+𝜏𝜏₀^22𝜆. (2) In this study, values of 𝜏₀, 𝜈₀, and 𝜆 are set to 0.050, 0.046, and 0.3, respectively, leading to a kernel function with a TF resolution of {Δ_𝑡,Δ_𝑓}={10.9s,0.039Hz}, where Δ_𝑡 and Δ_𝑓 quantify the spreading introduced by the kernel [25, 32].

Time-frequency coherence is an estimate of the strength of the local coupling between two signals and is determined by 𝛾_𝑥𝑦=𝑆𝑡,𝑓^𝑥𝑦𝑡,𝑓𝑆_𝑥𝑥𝑡,𝑓𝑆_𝑦𝑦𝑡,𝑓,𝛾_𝑥𝑦𝑡,𝑓∈[0,1]. (3)

Time-frequency phase difference (TFPD) is given by Θ_𝑥𝑦ℑ𝑆𝑡,𝑓=arctan^𝑥𝑦𝑡,𝑓ℜ𝑆_𝑥𝑦,Θ𝑡,𝑓_𝑥𝑦𝑡,𝑓∈[−𝜋,𝜋]. (4)

2.4. Partial Time-Frequency Analysis

The separation of respiratory influences from the RR interval series is conducted using partial TF spectra obtained by 𝑆_{𝑥𝑥/𝑦}𝑡,𝑓=𝑆_𝑥𝑥−𝑆𝑡,𝑓^𝑥𝑦𝑡,𝑓𝑆_𝑦𝑥𝑡,𝑓𝑆_𝑦𝑦=𝑡,𝑓1−𝛾²_𝑥𝑦𝑡,𝑓𝑆_𝑥𝑥𝑡,𝑓. (5) In this study, we focus on the partial spectrum of the RRV (𝑥=𝑅) from which the respiratory influences (𝑦=𝑟) are removed (𝑆_RR/𝑟(𝑡,𝑓)). Its complement, that is, the RRV which contains RR changes related to respiration (𝑆_RR,𝑟(𝑡,𝑓)), is defined as 𝑆_{𝑥𝑥,𝑦}𝑡,𝑓=𝛾²_𝑥𝑦𝑡,𝑓𝑆_𝑥𝑥𝑡,𝑓. (6)

In order to evaluate the performance of partial TF spectra, a simulation study is set up. Let 𝑥₁(𝑡) and 𝑥₂(𝑡) be two nonstationary signals, with each one being composed by one or two different complex exponentials showing both amplitude (𝐴_𝑘(𝑡)) and frequency (𝑓_𝑘(𝑡)) modulations: 𝑥_𝑖(𝑡)=^𝑀𝑖_𝑘=1𝐴_𝑘(𝑡)exp𝑗𝜙_𝑘(𝑡)+𝜃_𝑘 (7) with 𝑓_𝑘(𝑡)=(1/2𝜋)(𝑑𝜙_𝑘(𝑡)/𝑑𝑡), 𝑖={1,2}, and 𝑀_𝑖 is the number of spectral components. Let 𝑥(𝑡) be given as, 𝑥(𝑡)=𝑥₁(𝑡)+𝑥₂(𝑡)+𝜉(𝑡), (8) where 𝜉(𝑡) is complex zero-mean white Gaussian noise whose standard deviation is adjusted to obtain an SNR equal to 20 dB. In a first simulation study, 𝑀₁ and 𝑀₂ are, respectively, 2 and 1. The TF spectra 𝑆₁₁(𝑡,𝑓) and 𝑆₂₂(𝑡,𝑓) are given in Figure 1. In a second simulation study, an extra component is added to 𝑥₂(𝑡), yielding signals with TF spectra shown in Figure 2.

Figure 1 

Time-frequency spectra of signals 𝑥1(𝑡), 𝑥2(𝑡), and 𝑥(𝑡), see (7)-(8), of simulation study 1. Note that 𝑥1(𝑡) and 𝑥2(𝑡) do not share a spectral component.

Figure 2 

Time-frequency spectra of signals 𝑥1(𝑡), 𝑥2(𝑡), and 𝑥(𝑡), see (7)-(8), of simulation study 2. Note that spectral components of 𝑥1(𝑡) and 𝑥2(𝑡) intersect in the TF domain.

In the next step of the simulation study, the partial spectrum 𝑆_𝑥𝑥/2(𝑡,𝑓) is computed. The performance of partial TF spectra is evaluated based on the agreement between 𝑆_𝑥𝑥/2(𝑡,𝑓) and 𝑆₁₁(𝑡,𝑓). This agreement is quantified by the spectral distance 𝑑 between both TF spectra [33]: 𝑑=𝑆_𝑥𝑥/2𝑡,𝑓−𝑆₁₁𝑡,𝑓_𝑙1𝑆₁₁𝑡,𝑓_𝑙1 (9) with 𝑆_𝑥𝑥/2(𝑡,𝑓) being the mean partial TF spectrum estimated over all of the realizations of each simulation study and ‖𝐴‖_𝑙1 being the L1-norm of matrix 𝐴.

2.5. Derivation of Time-Varying Parameters

The time courses of several indices that characterize the interactions between RRV and respiration are derived. We compute the following time-varying parameters for TF spectrum 𝑆_𝛼(𝑡,𝑓)∈[𝑆_RR(𝑡,𝑓), 𝑆_RR/𝑟(𝑡,𝑓), 𝑆_RR,𝑟(𝑡,𝑓)] in specific frequency bands 𝛽, with 𝑓_𝛽 being the frequencies in 𝛽.(i)Instantaneous power is given as 𝑃^𝛽_𝛼(𝑡)=∑_𝛽𝑆_𝛼(𝑡,𝑓_𝛽)𝛿_𝑓 with 𝛿_𝑓 being the frequency step in the spectrum.(ii)Instantaneous frequency 𝐹^𝛽_𝛼(𝑡) is estimated as the frequency of the spectral peak in frequency band 𝛽. The considered bands 𝛽 are based on the traditional RRV frequency bands: (i)𝛽_LF=[0.04Hz,0.15Hz],(ii)𝛽_HF=[0.15Hz,0.40Hz],(iii)𝛽_TOT=[0.04Hz,0.40Hz].

In addition, two parameters that describe the local coupling between RRV and respiration in a specific frequency band 𝛽_𝑟 around respiration are computed as follows:(i)coherence 𝛾^𝛽𝑟_𝑅𝑟(𝑡)=mean_{𝑓∈𝛽𝑟}[𝛾_𝑅𝑟(𝑡,𝑓)],(ii)phase difference Θ^𝛽𝑟_𝑅𝑟(𝑡)=mean_{𝑓∈𝛽𝑟}[Θ_𝑅𝑟(𝑡,𝑓)],

where 𝛽_𝑟(𝑡) is a time-varying band defined as 𝐹_𝑟(𝑡)±(Δ_𝑓/2) with 𝐹_𝑟(𝑡), the respiratory frequency. This frequency band is centered around the breathing frequency and lies generally within the traditional HF band.

The last two time-varying parameters that are considered are the instantaneous respiratory frequency, 𝐹_𝑟(𝑡), and heart rate, HR(𝑡), expressed in beats per minute [bpm].

2.6. Statistical Analysis

The Wilcoxon signed rank test is used to assess statistical differences between the 4 conditions (AT, MT1, MT2, and RD). Because we are interested in the dynamic response to each task, statistical analysis was conducted sample by sample to track the 𝑃 values in time. Statistical significance is obtained when 𝑃<0.05.

In the case of the instantaneous frequency and power, we are only interested in the relative changes, regardless of the subject’s general condition or prior influences. Therefore, a correction was applied at the onset of each task prior to application of the Wilcoxon signed rank test to study only the relative response. The reference is taken as the mean value of each parameter in a window Δ_𝑡 around the onset of each task. No correction is applied for the coherence and phase difference.

3. Results

3.1. Simulation Study

Each simulation study is performed using 50 different realizations of noise and phases 𝜃_𝑘, derived from a uniform random distribution between −𝜋 and 𝜋. The presented results arise from averaging over these realizations. Figure 3 shows the averaged results of the simulation studies described in Section 2.4. In Figure 3(a), the TF spectrum 𝑆₁₁(𝑡,𝑓) is given. Figure 3(b) displays the partial TF spectrum 𝑆_𝑥𝑥/2(𝑡,𝑓) of the first simulation. A high similarity (𝑑=0.20) between 𝑆_𝑥𝑥/2(𝑡,𝑓) and 𝑆₁₁(𝑡,𝑓) is found in the first simulation.

(a)

(b)

(c)

Figure 3 

Averaged results of the simulation study: (a) TF spectrum 𝑆11(𝑡,𝑓); (b) partial TF spectrum 𝑆𝑥𝑥/2(𝑡,𝑓) of simulation study 1; (c) partial TF spectrum 𝑆𝑥𝑥/2(𝑡,𝑓) of simulation study 2. In the ideal case, the partial spectrum 𝑆𝑥𝑥/2(𝑡,𝑓) equals 𝑆11(𝑡,𝑓).

The result of the second simulation study is shown in Figure 3(c). We can observe “missing” parts around 𝑡=170 s and 𝑡=900 s, precisely for those portions of the TF domain in which spectral components of 𝑥₁(𝑡) and 𝑥₂(𝑡) intersect. A higher spectral distance 𝑑=0.22 is found, indicating a slightly lower agreement between 𝑆_𝑥𝑥/2(𝑡,𝑓), and 𝑆₁₁(𝑡,𝑓). Figure 4 shows the instantaneous power of the spectral components of 𝑥₁(𝑡) in simulation study 2 of 𝑆₁₁(𝑡,𝑓), 𝑆_𝑥𝑥(𝑡,𝑓) and 𝑆_𝑥𝑥/2(𝑡,𝑓), computed as the mean power in a window Δ_𝑓 around the frequency modulation 𝑓_𝑘(𝑡) from (7). Components 1 and 2 correspond, respectively, to the sinusoidal and linear components of 𝑥₁(𝑡) with a linear increasing and decreasing power (see 𝑃₁₁(𝑡)). Remark that the instantaneous power of 𝑆_𝑥𝑥(𝑡,𝑓) significantly increased around 𝑡=170 s and 𝑡=900 s, precisely when the instantaneous frequencies of 𝑥₁(𝑡) and 𝑥₂(𝑡) intersect. The instantaneous power of 𝑆_𝑥𝑥/2(𝑡,𝑓) decreased in correspondence to these intersections.

(a)

(b)

Figure 4 

Averaged instantaneous power of the two spectral components of 𝑥1(𝑡) in simulation study 2 of 𝑆11(𝑡,𝑓), 𝑆𝑥𝑥(𝑡,𝑓), and 𝑆𝑥𝑥/2(𝑡,𝑓). (a) component 1 of 𝑥1(𝑡), which corresponds to the sinusoidal component; (b) component 2 of 𝑥1(𝑡), which corresponds to the linear component.

3.2. Stress Monitoring

Figure 5 shows the RR interval series, respiratory signal, and TF spectra and coherence of one typical subject during the documentary watching (RD) and the first mental stress task. In the top 2 panels it is shown that both RR and RRV decrease when MT1 is compared with RD for this subject. An increase in the respiratory rate and an increased number of sighs are noted during mental stress. A high coherence around respiratory frequency is found during RD, while a reduction in coherence is noticed during stress for this subject.

(a)

(b)

Figure 5 

Example of RR interval series (RR) and respiratory (resp) signal and their TF spectra (𝑆RR(𝑡,𝑓),𝑆𝑟𝑟(𝑡,𝑓)) and coherence (𝛾𝑅𝑟(𝑡,𝑓)) during documentary watching and the first mental stress task. Vertical lines indicate the onset of RD (a) and MT1 (b).

Figure 6 shows the partial TF spectra for the same subject during RD and MT1. During RD, we observe that 𝑆_RR,𝑟(𝑡,𝑓) contains power in both the LF and HF bands. On the other hand, 𝑆_RR/𝑟(𝑡,𝑓) includes most power in the LF band. However, note that the HF band still comprises some power, though strongly reduced. 𝑃_RR,𝑟(𝑡) and 𝑃_RR/𝑟(𝑡) decreased during MT1 in comparison to RD.

(a)

(b)

Figure 6 

Example of partial TF spectra (𝑆RR,𝑟(𝑡,𝑓),𝑆RR/𝑟(𝑡,𝑓)) during documentary watching and the first mental stress task. Vertical lines indicate the onset of RD (a) and MT1 (b).

3.2.1. Time-Varying Parameters

Figure 7 shows the median instantaneous respiratory frequency (𝐹_𝑟(𝑡)) and heart rate (HR(𝑡)) in the top panels. In addition, the time instances of statistically significant differences between RD and the other tasks are indicated by bars below each subplot. Both 𝐹_𝑟(𝑡) and HR(𝑡) increase during AT, MT1, and MT2, while they slightly decrease during RD. Significant differences between RD and the other conditions are observed within 10 s after onset of each task, throughout the whole task for 𝐹_𝑟(𝑡). No differences are found between AT, MT1, and MT2. In the heart rate, significant differences are found a few seconds after onset of each task, between all conditions, except between AT and MT2. The heart rates differ maximally after 30 s from the onset of the task. After 100 s, no consistent differences between the 4 conditions are observed.

(a)

(b)

(c)

(d)

Figure 7 

Median instantaneous respiratory frequency (𝐹𝑟(𝑡)), heart rate (HR(𝑡)), coherence (𝛾𝛽𝑟𝑅𝑟(𝑡)), and phase difference (Θ𝛽𝑟𝑅𝑟(𝑡)) in the time-varying band 𝛽𝑟(𝑡). The standard error is shaded. The bars below each subplot indicate the time instances of significant differences between RD and the tasks.

The coherence (𝛾^𝛽𝑟_𝑅𝑟(𝑡)) and phase difference (Θ^𝛽𝑟_𝑅𝑟(𝑡)) are shown in the lower panels of Figure 7. Mental stress exhibits a reduction in cardiorespiratory coherence compared to RD and AT. This is only significant for MT1. No other differences in coherence between all conditions can be observed. The phase difference does not show any consistent difference as a result of mental stress or attention.

Figure 8 shows the median instantaneous power in the total frequency band. 𝑃^TOT_RR(𝑡) decreases during mental stress. This reduction is statistically significant already 10 s after onset of the mental task. Also a significant, but smaller, reduction is found during AT after 20 s. When considering the partial spectra, we observe a similar pattern in 𝑃^TOT_RR,𝑟(𝑡). 𝑃^TOT_RR/𝑟(𝑡) shows statistically significant differences between RD and MT1 from 40 s till 100 s, between RD and MT2 from 30 s till 120 s, and between RD and AT from 45 s till 80 s.

Figure 8 

Median instantaneous power for TF spectra 𝑆RR(𝑡,𝑓), 𝑆RR,𝑟(𝑡,𝑓), and 𝑆RR/𝑟(𝑡,𝑓) in the total frequency band. The standard error is shaded. The bars below each subplot indicate the time instances of significant differences between RD and the tasks.

Figure 9 shows the median instantaneous power for TF spectra 𝑆_RR(𝑡,𝑓), 𝑆_RR,𝑟(𝑡,𝑓), and 𝑆_RR/𝑟(𝑡,𝑓) in the LF, Figure 9(a), and HF, Figure 9(b), bands. Interestingly, 𝑃^LF_RR(𝑡) displays a significant reduction during MT2. This reduction is also found during MT1 and AT compared to RD, but to a lesser extent. A similar pattern is observed in 𝑃^LF_RR,𝑟(𝑡), suggesting that this reduction may be related to respiration; 𝑃^LF_RR/𝑟(𝑡) shows only a significant reduction during MT1 and AT, but not during MT2.

(a)

(b)

Figure 9 

Median instantaneous power for TF spectra 𝑆RR(𝑡,𝑓), 𝑆RR,𝑟(𝑡,𝑓), and 𝑆RR/𝑟(𝑡,𝑓) in the LF (a) and HF (b) bands. The standard error is shaded. The bars below each subplot indicate the time instances of significant differences between RD and the tasks.

In the HF band of 𝑆_RR(𝑡,𝑓), we observe significant differences between 20 s and 50 s when comparing MT1 and AT with RD. This reduction is also found in 𝑃^HF_RR,𝑟(𝑡) and to a lesser, but still significant, extent in 𝑃^HF_RR/𝑟(𝑡) (only for MT1), showing that the differences that appear in the HF band can be attributed to both RRV related and unrelated to respiration. In contrast with the findings in the LF band, MT2 shows no significant difference with any of the other conditions.

The median instantaneous frequencies for TF spectra 𝑆_RR(𝑡,𝑓), 𝑆_RR,𝑟(𝑡,𝑓), and 𝑆_RR/𝑟(𝑡,𝑓) in the HF band are given in Figure 10. The instantaneous frequencies in the LF band are not shown as they do not exhibit significant differences between the conditions. The top figure shows 𝐹^HF_RR(𝑡), which is expected to be highly coupled to respiration. This is confirmed when looking at 𝐹^HF_RR,𝑟(𝑡), which shows a highly similar pattern as 𝐹_𝑟(𝑡) in Figure 7. Also here, significant differences between AT, MT1, and MT2 and RD are found. 𝐹^HF_RR/𝑟(𝑡) shows no distinction between the tasks.

Figure 10 

Median instantaneous frequency for TF spectra 𝑆RR(𝑡,𝑓), 𝑆RR,𝑟(𝑡,𝑓), and 𝑆RR/𝑟(𝑡,𝑓) in the HF band. The standard error is shaded. The bars below each subplot indicate the time instances of significant differences between RD and the tasks.

4. Discussion

In this work, we aimed to characterize stress-related changes in cardiorespiratory interactions. We used time-frequency analyses to assess the cardiorespiratory coupling. In addition, partial TF spectra were introduced to distinguish between RR interval variations related and unrelated to respiration. This approach was also evaluated in two simulation studies. Next, several time-varying parameters were derived from the computed (partial) TF spectra and statistical analysis was conducted to assess whether different mental tasks provoked different cardiorespiratory responses compared to changes induced by a relaxing documentary watching task.

4.1. Cross and Partial Time-Frequency Analyses

Time-frequency analyses were conducted as they provide the time course of spectral indices which have an established physiological interpretation. The proposed TFD was used in previous research [25, 26, 32] and proved to provide a better TF resolution and localization than spectrogram and continuous wavelet transform of the local coupling between two signals.

The use of partial TF analyses, as described in Section 2.4, was evaluated in two simulation studies, where the performance was assessed by the similarity between 𝑆₁₁(𝑡,𝑓) and 𝑆_𝑥𝑥/2(𝑡,𝑓). In simulation study 1, none of the signal spectral components intersect; that is, they have no overlapping instantaneous frequencies. As shown in Figures 3(a) and 3(b), there is a high correspondence between 𝑆₁₁(𝑡,𝑓) and 𝑆_𝑥𝑥/2(𝑡,𝑓), which demonstrates that a good performance is obtained in the case that two signals do not share a spectral component.

In simulation study 2, the instantaneous frequencies of 𝑥₁(𝑡) and 𝑥₂(𝑡) overlap twice, around 𝑡=170 s and 𝑡=900 s, where 𝑃_𝑥𝑥(𝑡) increases and 𝑆_𝑥𝑥/2(𝑡,𝑓) is lower than 𝑆₁₁(𝑡,𝑓). A strong reduction in 𝑆_𝑥𝑥/2(𝑡,𝑓) is expected because in these portions of the TF domain, 𝛾²_𝑥𝑦(𝑡,𝑓)∼1, by definition, partialization is achieved by removing 𝛾²_𝑥𝑦(𝑡,𝑓)𝑆_𝑥𝑥(𝑡,𝑓) from 𝑆_𝑥𝑥(𝑡,𝑓). From this second simulation study, we can conclude that partial TF spectra cannot be used to separate the contribution of two spectral components which are simultaneously oscillating at the same instantaneous frequency. However, this limitation is expected to have little impact on the analysis of cardiorespiratory interactions for the following reasons: (1) RR oscillations locally coupled to and synchronous with respiration are assumed to come from respiration itself; (2) the hypothetical overlapping between RRV oscillations related to and not related to respiration is expected to be low, because the respiratory signal is relatively narrow banded; and (3) as long as the difference between the instantaneous frequencies of the spectral components of different signals is lower than the frequency resolution (Δ_𝑓), the partial spectrum correctly separates them. A critical evaluation on partial spectral analysis can be found in [34].

In addition, (5) shows that the amount of power in the partial TF spectrum 𝑆_{𝑥𝑥/𝑦}(𝑡,𝑓) which is removed from 𝑥(𝑡) around (𝑡₀,𝑓₀) is proportional to the coherence between 𝑥(𝑡) and 𝑦(𝑡) around (𝑡₀,𝑓₀). If at time 𝑡₀𝑥(𝑡) and 𝑦(𝑡) share a spectral component with the same instantaneous frequency 𝑓₀, their coherence at (𝑡₀,𝑓₀) will be close to 1 and thus the component will be removed, regardless of the amplitude modulation of both components.

Finally, it is important to note that the use of partial TF spectra is motivated by the fact that this approach is (1) nonparametric; that is, it does not rely on any model, and (2) it is fully data driven, in contrast to other techniques such as TF filtering: (1) model-based approaches, such as the IPFM model, MVAR models, or multivariate point process models, offer the opportunity to estimate hidden variables and the strength of directional couplings, but their outcomes strongly depend on the goodness-of-fit of the model. For this reason, we opted for a methodology which, in our view, is more robust to evaluate changes in RRV linearly related and unrelated to respiration. (2) TF filtering requires the use of a TF mask or a smoothing function, where it is necessary to decide on the geometry of the mask, its values, and, more importantly, its functioning. This makes the algorithm specific to a given set of signals in a given condition. Our approach, on the contrary, being data driven, can be generalized to any kind of signal in any possible situation and provides an easier interpretation of the results as the magnitudes of the partial spectra are proportional to the coherence between RRV and respiration.

4.2. Cardiorespiratory Response

The heart rate shows a clear increase due to mental stress and sustained attention. These findings are in agreement with those reported in [6, 9, 10, 23]. In line with the hypothesis, the largest increase is found during MT1, while a smaller increase is observed during MT2 and AT, suggesting a positive correlation between HR and mental load.

Also in line with our hypothesis is the decrease in the total power of RRV that is more pronounced during the mental tasks than during the attention task. Watching the documentary also gives rise to a brief reduction in total power immediately after onset of the task. These results suggest that mental stress causes vagal withdrawal, and sustained attention also a small reduction in vagal activity.

As previously reported in [9], the strongest reduction of HF power corresponds to MT1. Moreover, this study reveals that the reduction during AT is mainly due to the initial drop in HF between 20 and 50 s. The effect of MT1 is also strongest during this interval, but a small reduction persists throughout the whole task. In contrast to what was found in [9], no significant difference between MT2 and RD was observed, indicating that habituation might have occurred during MT2. The influence of the second mental stress task is only seen in the LF band, which corresponds to both sympathetic and vagal influences. A small reduction is found during AT and MT1. These results are in contrast with the findings by Taelman et al. using the original RR interval series [9]. They found a strong reduction in LF power during MT1 and AT, but MT2 did not differ significantly from RD. When conducting TF analyses, they did report significant differences between AT and MT2, but no differences between AT and MT1. We can observe a similar pattern; that is, the largest difference is between AT and MT2. However, these differences are not significant in this study. The differing results might arise from the applied correction to study only relative changes, thereby discarding possible effects from prior tasks. This was not implemented in [9].

The respiratory frequency increases during AT and MT, which was also reported in [14]. This is an important parameter to take into account as many authors suggest that an inverse relationship between RSA magnitude and respiratory rate exists, which is independent of vagal activity [19–22]. Seeing that HF power is often taken as a measure of RSA, it is of utmost importance to be able to distinguish changes in HF due to vagal control or due to changes in respiratory rate. The same reasoning holds for the depth of breathing, as it is also shown that this influences RSA magnitude independently of vagal outflow. In this TF framework, information on the depth of breathing can be obtained in 𝑃^𝛽𝑟_𝑟𝑟(𝑡). We found, however, that, in this study, the depth of breathing did not differ significantly among the tasks or compared to RD.

We observed that the coherence slightly decreases during stress. This effect as a result to mental stress was hypothesized, due to vagal withdrawal and sympathetic activation [17]. The phase difference remains constant during the different conditions and is thus not affected by stress or attention.

4.2.1. Partial TF Analysis

Figure 9 shows that there are similar time courses in the 4 conditions between 𝑃_RR(𝑡) and 𝑃_RR,𝑟(𝑡) in the LF and HF bands. This indicates that both the RR interval variations in the LF and HF bands are highly coupled to respiration and demonstrates the dominant effect of respiration on RRV, as can also be observed in the instantaneous high frequency in Figure 10. However, the question rises regarding what the physiological interpretation of this cardiorespiratory coupling in the LF band is. Note that during the computation of indices from the LF band we verified that the respiratory frequency did not fall within this band. If we did encounter slow respiratory rates, we did not compute the indices at those time instances. Possibly, the respiratory influences in the LF band can be explained by the baroreflex feedback theory or irregular breathing patterns, as was suggested by Yildiz and Ider [35].

Also in the total band, the effect of respiration is apparent. Moreover, the effects of stress in 𝑃^TOT_RR,𝑟(𝑡) are already pronounced after 3–7 s, while in 𝑃^TOT_RR(𝑡) a significant difference was only found after 10 s. Both in 𝑃^TOT_RR(𝑡) and 𝑃^TOT_RR,𝑟(𝑡) the effect of sustained attention is detected after 20 s.

As expected, most of the power in the HF band is related to respiration. However, there is still some HF power in 𝑆_RR/𝑟(𝑡,𝑓), which furthermore shows effects related to stress. Opposed to what we hypothesized, 𝑃^LF_RR/𝑟(𝑡) does not show a clear discrimination between the conditions. Only a slight reduction is found during MT1 and AT. The analysis of 𝑃^TOT_RR/𝑟(𝑡) exhibits differences between RD and the other tasks; the effects of sustained attention appear only for a short interval, from 45 to 80 s, while the influence of stress on physiological processes, other than respiration, is observed after 30 s and lasts till 100 s.

The results from partial TF analysis suggest that the effects of stress are mainly related to respiration. Although it was hypothesized that RRV unrelated to respiration would show enhanced stress influences, as was found during classification of rest and stress [29], we found only small differences between the active tasks and RD. The partial TF analyses revealed different temporal patterns as a result of stress; a fast response (within 7 s) is observed in RRV related to respiration, as shown in Figure 8, while RR interval variations unrelated to respiration are only apparent after 30 s.

Finally, it is important to note that the results not only display different responses to mental stress and attention, but also differences between MT1 and MT2 can be observed. Although both tasks are the same, they present a different response, in terms of magnitude and latency. We hypothesized that the response to stress would be faster and stronger during MT1, while during MT2 the effect of stress would be reduced due to habituation, as was also observed by Taelman et al. [9]. This hypothesis is confirmed in terms of heart rate, coherence, HF, and total power.

5. Conclusion

The goal of this study was to characterize the dynamic interactions in the cardiorespiratory regulation in response to mental stress. As a suitable approach, cross time-frequency analyses were conducted. In addition, partial TF spectra were computed to evaluate separately the response of RR interval variations linearly related to respiration and variations that are not linked to respiration. Although stress also influences the respiratory pattern [14], breathing is also under voluntary control and might not always be a suitable indicator of stress.

Sustained attention was adopted as a nonstressful task during which we observed an increased heart and respiratory rate and a slightly increased coherence. A transitory reduction of HF power suggests a vagal withdrawal during the first minute. Partial TF analyses showed that the response to sustained attention of RR interval variations related to respiration differs from that to documentary watching after 20 s, while RRV unrelated to respiration exhibits a significant response after 45 s. After 80 s, there is no effect detected of this attention task.

We found that mental stress causes an increase in heart rate and respiratory rate. Also a reduction in cardiorespiratory coherence, HF power, and LF power was found, indicating vagal withdrawal. The partial TF analysis revealed that the response to stress of RR interval variations related to respiration appears very fast (3 s), while the variations unrelated to respiration react with a slower temporal pattern (30 s). These results demonstrate that cross and partial time-frequency analyses carry valuable information on the cardiorespiratory stress mechanisms and suggest that this is a useful tool for biofeedback in stress-reducing therapies.

Acknowledgments

This research was supported by (i) Research Council KUL: GOA MaNet, PFV/10/002 (OPTEC), through several Ph.D./postdoc and fellow Grants; (ii) Flemish Government: (a) FWO: Postdoc Grants, Projects nos. G.0427.10N (integrated EEG-fMRI), G.0108.11 (compressed sensing), G.0869.12N (tumor imaging), G.0A5513N (deep brain stimulation); (b) IWT: Ph.D. Grants, Projects nos. TBM 070713-Accelero, TBM 080658-MRI, TBM 110697-NeoGuard; D. Widjaja is supported by an IWT Ph.D. Grant; (c) iMinds: SBO dotatie 2013, ICONs: NXT_Sleep, FallRisk; (d) Flanders Care: Demonstratie Project Tele-Rehab III (2012–2014); (iii) Belgian Federal Science Policy Office: IUAP P719 (DYSCO, 2012–2017); ESA AO-PGPF-01, PRODEX (CardioControl) C4000103224; (iv) EU: RECAP 209G within INTERREG IVB NWE Programme, EU HIP Trial FP7-HEALTH/2007–2013 (no. 260777), EU MC ITN TRANSACT 2012 (no. 16679), ERC Advanced Grant: BIOTENSORS (no. 39804), ERASMUS EQR: Community Service Engineer (no. 539642-LLP-1-2013). The scientific responsibility is assumed by the authors of this paper.