Quantifying neurophysiological mechanism through heart rate variability: the case of cognitive stress

Prolonged exposure to daily stress results in a continuous and destructive effect on the human mind. Cognitive dysfunctions and depressions are the various effects of the stress, which may also lead to death (Saxena et al., 2002). It has been reported that 50% of the employees in the world face several work-related stresses that lead to depression as well (Engin, 2004). Several studies show that prolonged exposure to stress may lead to decreased quality of life and acquired acute depression, which affect mental health (Arzeno et al., 2006; Israel et al., 2005).

Nowadays, individuals are experiencing significantly more mental anxieties than in the past. While various reasons contribute to the current state, the fundamental driver would be the inaccessibility of unwinding time and the consistent need for difficult mental work. Cognitive tasks, which are relatively harder, cause frustration due to underestimation of their time of completion, which ultimately results in anxiety and stress (Badre, 2021). Persistent stress is a significant contributor to morbidity and mortality (Dossett et al., 2020).

In this study, we will be focusing on the effect of cognitive stress on the autonomous nervous system (ANS). The ANS helps us in realizing instantaneous motor functions based on external stimuli (Israel et al., 2005; Saxena et al., 2002). Thus, it will be considered an essential factor to measure acute stress in a non-invasive way, which is known to be fundamental and reliable for measuring the activity of ANS (Bansal et al., 2007). Control of this activity is done by sympathetic and parasympathetic responses of ANS and can be determined by the responses of heart rate variability (HRV) signals (Dobbs et al., 1984). HRV can be used to measure the mental anxiety of normal healthy individuals by comparing them at rest and at stressful conditions. It was reported that even a subtle change in cognitive activity can be detected through HRV (Choi, 2017); no wonder HRV is regarded as a clinical index for assessing psychophysiological resilience (Mourtakos et al., 2021).

The sinoatrial node (SA node) refers to a bundle of specialized cells in the heart responsible for initiating heartbeats. HRV refers to the alterations of the time intervals between subsequent heartbeats. The variations are the result of different complex physiological functions. Heart–brain interaction and autonomic nervous system dynamics are probably the most dominating physiological function behind HRV. The heart–brain connection pathways are explained in Figure 1.

Figure 1:

Neurophysiological mechanism of HRV under cognitive stress. HRV, heart rate variability.

There are two sets of nerves, sympathetic and parasympathetic (vagus), which, in combination, maintain the equilibrium of the heart rate (HR) by slowing down or speeding up the rate, respectively. During routine tasks, rest, or sleep, HR is lower than the intrinsic rate, and parasympathetic activity predominates. When a stressful situation emerges, characterized by an HR >100 bpm, sympathetic activity takes over.

The descending sympathetic nerve targets the SA node via the underlying cardiac nervous system and much of the myocardium. Action potentials conducted by these motor neurons trigger the release of noradrenaline and adrenaline, increasing HR and enhancing atrial and ventricular contractility. After the onset of sympathetic stimulation, there is a delay of up to 5 s before the heart rate gradually increases with stimulation, reaching a steady level within 20–30 s of continuous stimulation. Even brief sympathetic stimulation can affect HR rate and HRV rhythm for 5–10 s. The relatively slow response to sympathetic stimulation contrasts with near-instantaneous vagal stimulation. Therefore, an increase or decrease in heart rate, or an abrupt change between heart rates, is mediated primarily by parasympathetic nerves.

Although the sympathetic nervous system (SNS) is dominant during stress, including cognitive one, there is yet another emerging perspective, especially when the stress is related to a cognitive task. This is the vagal reactivity and the recuperation measures that have been connected to intellectual execution (Capuana et al., 2014). Vagal reactivity addresses the change among standards and a particular occasion, such as finishing an undertaking, and it is fundamental to chronicle it to assess the person’s flexibility to the circumstance (Laborde et al., 2018). Recuperation is typically seen as a cycle of rebuilding; it alludes to the change between the occasion and a period point after the occasion when the vagal action must be like the pattern. Equivalent to vagal reactivity, vagal recuperation assumes a significant part in the versatility of the creature (Sinha et al., 2015). These two perspectives are ineffectively investigated with intellectual working. In any case, as indicated by the vagal tank hypothesis (Laborde et al., 2018), considering the vagal action and the vagal recuperation during various psychological assignments could be intriguing. This kind of study could permit to see better what heart vagal control means for a few key self-administrative parts of conduct (Delmastro et al., 2020; Engert et al., 2017; Ghiasi et al., 2020; Holzman & Bridgett, 2017; Ishaque et al., 2021; Pykett et al., 2020; Rachakonda et al., 2020; Rudovic et al., 2018). Furthermore, the study is indicative of the evident contrasts between benchmark, task execution, and recuperation, which are identified with psychological debilitation (Betti et al., 2018; Capdevila et al., 2021; Chen et al., 2020; Hong et al., 2020; Koren et al., 2020; Masood & Alghamdi, 2019; Ni et al., 2021; Saini & Gupta, 2019; Wu et al., 2019).

The meta-analysis includes a series of papers that are listed in the zone of time-domain of HRV analysis. Some of the measuring criteria are sorted out according to their relevance of use in different articles. Such zonal criteria are the mean of RR intervals (MeanRR), standard deviation of the RR intervals (RRStdd), square differences between root means square of the two successive RR intervals (RRMmsds), and proportion of RR intervals that differ more than a time period of 50 ms (NN50p), as depicted in Table 1.

The rigorous study that we have conducted leads to the fact that the MeanRR, NN50p, and RRMmsds get an enormous hike during the stressed condition, although some observations did not notice any difference in these features. By contrast, a reduced RRStdd was majorly observed in most of the studies (Chandola et al., 2010; Gabella, 2012; Grossmann et al., 2016; Keller et al., 2011; Taelman et al., 2011; Tharion et al., 2009), involving time-domain analysis. Exceptions include the study by Schubert et al. (2009), which came up with an observation completely contradictory. Our general observation shows a subjugated incremental behavior of the RRStdd function (Goldberg & Grandey, 2007; Grossman & Taylor, 2007; Grossmann et al., 2016). Justification of this contradictory theory is also given by them, where the authors mentioned that a slow respiratory schematic and relative decrease in the ventilating mechanism was the main cause for achieving this result during the activity of public speech delivery.

There are several papers that are shortlisted for quantification in frequency-domain analysis. The analytical feature of these quantifications is done on the basis of three regional distributions of the frequency-dominated signal. First, the signals are separated on the basis of low-frequency power area (LowF) and then segregated with respect to high-frequency power area (HighF), and finally, as a method of comparison, by the ratio between the HighF and the LowF.

Table 2 depicts the observation under such criteria. Among the papers listed, very few stated that there is a considerable amount of increase in LowF with a statistical significance, whereas the other observation was completely contradictory with this fact. Tharion et al. (2009) observed perfectly normal general trends for the situations involved. Although, the fact has to be considered that the assessment of the readings done is taken on the day of the examination, the readings during the cause examination were completely ignored.

In accordance with that, Hjortskov et al. (2004) presented us with a completely different histology where there is an introduction session continued with the rest and stress situations. In the end, a noticeable dominant decrease in the low-frequency spectrum is vividly seen during the stress situation. To our knowledge, this situation arises due to the fulfillment of different protocols to simulate stress within the subject. Taelman et al. (2011) also support the factors behind these considerations.

In the HighF power domain, several studies agree with the fact that there will be a reduction in measure during acute mental activity or stress. An exception to this fact has also been registered in the study (Vuksanović, 2007).

Studies registered and analyzed on the basis of the HighF to LowF ratio criterion agree with the same fact unanimously that the signal measures increased during stress conditions (Tharion et al., 2009; Vuksanović, 2007), whereas the rest of the studies showed a different statistical analysis, which is not considered in conventional statistical analysis.

A combined time-frequency metric is sometimes preferred to capture the temporal variation of the frequency components. These include short-time Fourier transform (STFT) and wavelet-based metrics.

In the STFT method, the segmented signal of the mathematical expression in the continuous time-domain is given by Y(t_C,t) = y(t)w(t_C − t), where t_c is in the middle of the real-valued, symmetric window function. The mathematical equation after using the Fourier transform on the segmented signal to restore its frequency properties is shown in the expression Ytc,wc=∫−∞∞Ytc,te−jwctdt Y\left( {{t_c},{w_c}} \right) = \mathop \int \nolimits_{ - \infty }^\infty Y\left( {{t_c},t} \right){e^{ - j{w_c}t}}dt .

For QRS complex detection, a wavelet-based approach is generally preferred. Some preprocessing such as noise removal and removal of baseline wandering is required before the wavelet transform can be applied. Window detection by the method of determining QRS complex is given by Wm,n=∫−∞∞f(y)φm,n*φm,nydy W\left( {m,n} \right) = \mathop \int \nolimits_{ - \infty }^\infty f(y)\varphi _{m,n}^*{\varphi _{m,n}}\left( y \right)dy , where the function of RR interval is given as Wt=∑i=1nδt−∑n=1iytn W\left( t \right) = \sum\limits_{i = 1}^n {\delta \left[ {t - \sum\limits_{n = 1}^i {yt\left( n \right)} } \right]} .

The non-linear analysis involves the application of non-linear system theory, such as chaos theory, to quantify non-stationary HRV signals. The HRV dynamics can be characterized in the non-linear domain, maintaining clinical utility through Poincare plots, short-term fractal scaling exponent measured with detrended fluctuation analysis (DFA), the Lyapunov exponent, and different entropy measures such as Shannon entropy (ShnEn) and sample entropy.

a.

Poincare plot and recurrence plot

Standard deviation in the case of the Poincare plotting system and the recurrence method are linear in nature and defined by SDstdd and SDrec, respectively. More than three state space trajectories can be reduced down to only two dimensions by this method of recurrence plot. The methodology is depicted as R_i,j = 0,w_i − w_j > ɛ;R_i,j = 1,w_i − w_j ≤ ɛ. Here in the formulation, only two states i and j are shown although it can be done for the higher number of states. The method needs a threshold value ɛ that is user defined according to the situation. Lengths of the distance between two trajectories are calculated, and if the value is less than the threshold it will be taken as 1 while if it is greater than ɛ then it will be taken as 0. Table 3 explains different measures of recurrence plot and their definitions.

The determination of the recurrence plot like Figure 4 is given by RECdet (Hautala et al., 2009; Lehrer et al., 2006) along with that the recurrence rate is also defined as RECrate in an approved condition, where the length of lines (Lmax) can be considered as a defining criterion, while the lines’ mean length (Lmean) plays a vital role in that system.

Figure 4:

Recurrence plot of the data analyzed.

The quantification of two-dimensional (2D) and three-dimensional (3D) Poincare plotting are given as e+=Ejk=1+τN,τ=1,⋯,N−1 {e^ + } = \left| {E\left( j \right)} \right|_{k = 1 + \tau }^N,\tau = 1, \cdots ,\left( {N - 1} \right) , where Ejj=1N \left| {E\left( j \right)} \right|_{j = 1}^N is a representation of a discrete signal, and e⁺ and e⁻ are subgroups of the discrete signal with time delay τ. The system of coordinates is replaced with a two-dimensional rotation with an angle of π/4 with respect to X-axis, i.e., emen=cosπ/4−sinπ/4sinπ/4cosπ/4e+e− \left( {\matrix{ {{e_m}} \cr {{e_n}} \cr } } \right) = \left( {\matrix{ {\cos \left( {\pi /4} \right)} & { - \sin \left( {\pi /4} \right)} \cr {\sin \left( {\pi /4} \right)} & {\cos \left( {\pi /4} \right)} \cr } } \right)\left( {\matrix{ {{e^ + }} \cr {{e^ - }} \cr } } \right) .

On the other hand, 3D Poincare plotting has been done by the method of ellipsoid fit, shown in Figure 5.

Figure 5:

Ellipsoid fit applied on the dense region.

In the new coordinate used in the system, Em′ = Mean Em, En′ = Mean (En), SD1 = [Var(Em)]^1/2, SD2 = [Var(En)]^1/2, where Ejj=1N \left| {E\left( j \right)} \right|_{j = 1}^N is a representation of discrete signal and e⁺, e⁻, and e⁻⁻ are subgroups of the discrete signal with time delay τ. Further, e+=Ejj=1N−2τ {e^ + } = \left| {E\left( j \right)} \right|_{j = 1}^{N - 2\tau } , e−=Ejj=1+τN−τ {e^ - } = \left| {E\left( j \right)} \right|_{j = 1 + \tau }^{N - \tau } , and e−−=Ejj=1+2τN {e^{ - - }} = \left| {E\left( j \right)} \right|_{j = 1 + 2\tau }^N , with τ = 1,⋯,(N − 1).

b.

Entropy

A mathematical formulation of ShnEn can be obtained as follows. The definition of energy for discrete-time ECG signal, x[n], may be represented in the following manner (Beyramienanlou & Lotfivand, 2017): Ex=∑−∞∞xn⋅xn*=∑−∞∞xn2 {E_x} = \sum\limits_{ - \infty }^\infty {\left[ {{x_n} \cdot x_n^*} \right]} = {\sum\limits_{ - \infty }^\infty {\left[ {{x_n}} \right]} ^2} . To find ShnEn, zero-phase filter and Shannon energy approximate are playing a vital role. Considering a[n] as the filtered and normalized version of x[n], the energy of the local spectrum for an individual sample is ShnEn = −|a[n]| log(|a[n]|). The process output, s[n] = −a²[n] log(a²[n]), is further normalized, s[n] = (s[n] − μ)/δ, where μ is the process normalization vector. Next, a rectangular filter of length h is used to diminish the spikes. The filtering operation is m[n] = filter(h, j, S), where m[n] defines moving average. The diagnosed peaks are linked together by the difference equation, d[n] = f[n] − f[n − 1], n = 2, 3, 4, …, where f[n] is the input to the Shannon energy calculation process.

To detect QRS complex in non-linear domain, the threshold selection has to be made by the following relations: Threshold = |kμ(1 − σ²)| if σ < μ, or Threshold = |kσ(1 − μ²)| if σ > μ, where k is a constant. The largest Lyapunov exponent is LowEmax, where the dimension correlation is D2. The fluctuation correlation can be done in short- and long-term modes.

The approximate entropy, En(rchon), was devised by Chon et al. (Chon et al., 2009), where the tolerance threshold (r) has been set to 0.03 NNstdd to obtain the maximum entropy value. These calculations by the embedding dimension methodology are indicated in Figure 6. The methodological technique is as follows. In this case, we initialize from a scalar time series placed in a m-dimensional phase space using a time-delay procedure, where the number of points placed in the m-phase space is N = N_c − (m − 1)h, where h is the time-delay represented as multiple of time interval. The correlation between them is expressed by: Cmr=limn→∞2NN−1∑i=1N∑j=i+1NΘr−yi−yj {C_m}\left( r \right) = \mathop {\lim }\limits_{n \to \infty } {2 \over {N\left( {N - 1} \right)}}\sum\limits_{i = 1}^N {\sum\limits_{j = i + 1}^N {\Theta \left( {r - \left| {{y_i} - {y_j}} \right|} \right)} } , where r is the correlation length. In most of the conditions, C_m(r) grows as a power of r. Thus, a plot of log(c_m(r)) vs log(r) shows the scaling region, log(C_m (r)) ≈ D2log(r). By this methodology, we can define the number of false neighbors in the system to its embedding dimension (D2; Casaleggio et al., 1995).

Figure 6:

False nearest neighbor test.

The two features, namely, entropy at the maximum threshold, En(rmax), and the max Entropy reported by Chon et al., En(rchon), are explored in various studies (Melillo et al., 2011; Papousek et al., 2010), which reported a decrease in Stdd under stress. Melillo et al. (2011) and Schubert et al. (2009) have taken their signals while performing an academic examination and arithmetic task, respectively. In a different paper, the same author (Melillo et al., 2013) reported a completely different trend, and the study by Vuksanović et al. (2007) confirmed the opposite trend. Different characteristics of non-linear parameters are expressed in Table 4.

Our current sutdy presents a schematic literature review of investigated articles in a meta-analytic way. We present how HRV signals are manipulated under artificial mental stress.

In the time-domain, the report measures four HRV metrics, namely RRStdd, RRMmsds, NN50p, and D2, which depict a unique behavior during stress conditions (Chandola et al., 2010; Li et al., 2009; Papousek et al., 2010; Schubert et al., 2009; Taelman et al., 2011; Tharion et al., 2009). The sample selection of the study coincides with the fact that there is a major decrease in the normal properties of RR signals during the zone of stress. Specifically, RRStdd (Schubert et al., 2009) was the measure where there is a significant decrease shown during stress.

The same characteristics have also been noticed in the HighF frequency-domain analysis (Chandola et al., 2010; Hjortskov et al., 2004; Li et al., 2009; Papousek et al., 2010; Taelman et al., 2011; Tharion et al., 2009), while the HighF/LowF ratio and the LowF itself have shown a significant behavior, during stress, which is increasing in nature (Chandola et al., 2010; Hjortskov et al., 2004; Kofman et al., 2006; Lackner et al., 2011; Li et al., 2009; Papousek et al., 2010; Traina et al., 2011; Vuksanović, 2007). HRV measures are mostly taken for 5 min (Schubert et al., 2009), while in this case, it has been taken for 3 min (Chandola et al., 2010; Taelman et al., 2011; Tharion et al., 2009) by the respective researchers. Coming to the domain of LowF, certain experts’ studies pointed out a major lowering in the behavior during stressed conditions (Hjortskov et al., 2004; Taelman et al., 2011; Tharion et al., 2009), while others reported a different scenario (Lackner et al., 2011; Li et al., 2009; Papousek et al., 2010; Traina et al., 2011; Vuksanović, 2007). Although their mode of study is completely different from the others, the first and the second ones used physical stress to mimic stress conditions, while the others used arithmetic tests to simulate stress. In the first work, Hjortskov et al. (2004), the volunteers were asked to perform a typing task with their dominant hand on the non-dominant section of a keyboard at a relatively fast pace under a time constraint. On the contrary, to simulate stress, participants have to choose a correct answer from the alternatives given to a question by clicking the mouse but not with their dominant hand. This activity seems to activate a different article that had not been activated when they were asked to use their usual working hand, as depicted by Taelman et al. (2011). An important observation is that no physical activity has been involved in this process. A majority of the studies asked their subjects to perform a mathematical test using only the keyboard. This experiment has been performed by Yu and Zhang (2012). On the other hand, a single study reported a decreased HRV during stress with a LowF value less than the usual stress condition. In this situation, the performers were the students of the university, and their HRVs are jotted on the day of their examination by Tharion et al. (2009) and Acharya et al. (2006). This experiment focused on the examination day but failed to report the conditions of the pupil during the examination (Balda et al., 1977; Melillo et al., 2011). These properties of the layering of stress-building procedures were delivered in Bracha (2004). Each of these layerings has a different and unique impact on HRV modulation when studied on ANS. The values of HRV measures are projected in different aspects in two different unique papers (Melillo et al., 2011; Vuksanović, 2007). Both of these papers have significant statistical significance. But the clarity of such heterogeneous results was not identified in this paper. The measures that are tabulated can be considered as a checkpoint for further researchers in this area (Bando et al., 2015; Bland & Altman, 1996; Blaug et al., 2007; Castaldo et al., 2015, 2016; Crowley et al., 2011; Jelinek et al., 2011; Melillo et al., 2012; Pan & Li, 2007; Raaijmakers et al., 2013; Satya, 2009; Siegel & Castellan, Jr., 1988; Sinha et al., 2016).

HRV measures like RRStdd, RRMmsds, NN50p, D2, HighF, LowF, and HighF/LowF are the seven significant metrics through which meta-analysis of the pooled values has been done. Thus, the values pooled together can be considered more reliable rather than picking on the values from the individual paper. Although the values received by the HighF/LowF measures are statistically significant, the differences in measurement were too small when compared to the measures of RRStdd during rest as well as stress (Singh et al., 2013). The values determined by the LowF measures can be considered effective if it is applied to different system designs (Hjortskov et al., 2004; Taelman et al., 2011).

From the above-measuring aspects, the non-linear analysis was found to be the most robust and apt system, as every analytical procedure in the non-linear domain shows similar results for all measures. From the study, it can be concluded that an analytical study of the cardiac signals can prove to be an elementary tool for diagnosis (Bishop, 2006; Eckberg, 1983; Garcia-Mancilla & Gonzalez, 2015; Hernández-Vicente et al., 2021; Liu & Ulrich, 2014; Salahuddin et al., 2007; Sano & Picard, 2013; Sinha et al., 2015; Tabar & Halici, 2016; Timothy et al., 2016). Moreover, the effect of human stress on the chaotic biological signals can be effectively measured with the aid of non-linear-domain analysis of the same (Barreto et al., 2007; Bousefsaf et al., 2013; Healey, 2015; Kelsey et al., 2004; Laird et al., 2005; McDuff et al., 2014; Nunan et al., 2010; Quintana & Heathers, 2014; Thayer et al., 2012). The accuracy of the non-linear metrics can significantly benefit the medical community.

We would like to conclude that the papers investigated in this study show a unanimous report on the non-linear measures of HRV, and they all proved that there is a significant deviation in it during acute mental stress. In the case of non-linear studies, where the same stress has been incorporated to manipulate the situation, the conclusions achieved are much more robust and give better results than the frequency-domain analysis. The situation can be simulated, or the stress condition may occur naturally, but the results in the non-linear domain are more accurate than that of the frequency-domain without any changes in the condition. The situations support the prior observations in this area, which proves that a shift has always been incorporated in ANS, maintaining a balanced equilibrium in the parasympathetic activation while maintaining avoidance while the subject is under stress. As we have seen, these non-linear studies can efficiently point out those shifts. Although various case studies studying the psychological effects of human stress have been incorporated in this study, a few more methodologies have not been used in the field of HRV analysis. Thus, this study provides a rigorous benchmark for future researchers to work on this field of non-linear activities of HRV and its characteristics during mental stress.

Although we have given a brief idea about the techniques used for HRV analysis, we have not developed any predictive model or analytical technique to support a correlation between the processes involved. To defeat these impediments, in future psychophysiological experiments, it will be helpful to use the arising rules for detailing HRV boundaries (Laborde et al., 2017), which can improve the nature of information, permit more straightforward announcing, and lead to more analyzable information in the quantitative investigation (e.g., meta-examination). Further examination should plan to build the investigations on the connection between HRV and some psychological areas, like consideration, language, preparing speed, and visuospatial abilities that are dismissed by the examinations up to this point.