Mind-body relationships in elite apnea divers during breath holding: a study of autonomic responses to acute hypoxemia.

The mental control of ventilation with all associated phenomena, from relaxation to modulation of emotions, from cardiovascular to metabolic adaptations, constitutes a psychophysiological condition characterizing voluntary breath-holding (BH). BH induces several autonomic responses, involving both autonomic cardiovascular and cutaneous pathways, whose characterization is the main aim of this study. Electrocardiogram and skin conductance (SC) recordings were collected from 14 elite divers during three conditions: free breathing (FB), normoxic phase of BH (NPBH) and hypoxic phase of BH (HPBH). Thus, we compared a set of features describing signal dynamics between the three experimental conditions: from heart rate variability (HRV) features (in time and frequency-domains and by using nonlinear methods) to rate and shape of spontaneous SC responses (SCRs). The main result of the study rises by applying a Factor Analysis to the subset of features significantly changed in the two BH phases. Indeed, the Factor Analysis allowed to uncover the structure of latent factors which modeled the autonomic response: a factor describing the autonomic balance (AB), one the information increase rate (IIR), and a latter the central nervous system driver (CNSD). The BH did not disrupt the FB factorial structure, and only few features moved among factors. Factor Analysis indicates that during BH (1) only the SC described the emotional output, (2) the sympathetic tone on heart did not change, (3) the dynamics of interbeats intervals showed an increase of long-range correlation that anticipates the HPBH, followed by a drop to a random behavior. In conclusion, data show that the autonomic control on heart rate and SC are differentially modulated during BH, which could be related to a more pronounced effect on emotional control induced by the mental training to BH.

1. Introduction

The voluntary apnea is a psychophysiological condition characterized by the cognitive control of ventilation. As competitive sportsmen, breath-hold divers aim at prolonging the duration of their performance as long as possible. In general, a number of factors may influence breath-holding(BH) time, such as hypoxia and hypercapnia (chemical factors), lack of respiratory movements, and lung capacity (mechanical factors). In fact, practicing apnea induces physiological and psychological adaptations responsible for the long-lasting of the BH. The physiological ones include a desensibilization to CO2 increases (Ferretti, 2001; Andersson et al., 2004) as well as a minimization of the involuntary contractions of respiratory muscles provoked by the absence of respiratory movements (Delapille et al., 2002); on the psychological side, an adaptation to tolerate the need of breathing has been reported (Hentsch and Ulmer, 1984).

Apnea is commonly divided into two phases: an initial normoxia lasts till a physiological breaking point, when the O2 saturation has a drop, then a hypoxic condition takes place, during which the subject feels a growing urge to breathe. In the latter phase, activity of the diaphragm shows, after a period of sporadic activity, progressive involuntary movements followed by major contractions reaching their maximum when breath-hold becomes unbearable and the subject reaches his limit, and starts breathing again. These respiratory muscle contractions are observed both in untrained and in elite divers; however, trained subjects learn to tolerate the involuntary respiratory muscle contractions which allows them to postpone the limit point so as to enhance their BH time (Lin et al., 1983). In other words, during BH, chemical stimulation sends information to the respiratory center, but it appears that the greater effect on the respiratory center and respiratory muscle is indirectly sent by nonchemical stimuli; this components appear to be weaker in experienced subjects, permitting the increase of BH time (Delapille et al., 2001).

As a regular, continuous supply of oxygen is essential to sustain life, human beings possess innate defense mechanisms to raise tolerance to acute and/or chronic hypoxemia. In the particular case of apnea divers, the most important physiological adaptation to the oxygen lack is the so-called “diving reflex” (Ferretti, 2001; Andersson et al., 2004; Lindholm and Gennser, 2005). This reflex makes the breath-hold divers tolerate extreme hypoxic/hypercapnic condition thanks to ventilatory, cardiovascular, and cerebrovascular adaptations, minimizing the tissue oxygen uptake and to improve the dive duration. Cardiovascular modifications include peripheral vasoconstriction through sympathetic activation, connected with early hypertension and bradycardia caused by an increased vagal activity associated with a reduction of the cardiac output (Fagius and Sundlf, 1986; Gooden, 1994; Lindholm and Gennser, 2005). Before the physiological breaking point cardiac output is, therefore, preferentially directed to vital organs, and its reduction is predominantly due to an increase of intrathoracic pressure which delays the venous return and decreases the stroke volume (Ferrigno et al., 1986). Several factors influence these cardiovascular responses initiated by apnea and augmented by face immersion, such as arterial hypoxia perceived by chemoreceptors, the related hypercapnia, cessation of respiratory movements, and facial cold receptors. Also, recent studies propose that an active contraction of the spleen is also a feature of the human diving response and may play an O2-conserving role (Hurford et al., 1990; Schagatay et al., 2001). In absence of diving reflex, for “dry” apnea, some studies have found neither bradycardia nor a decrease in the cardiac output in response to BH (Perini et al., 2008; Pingitore et al., 2008).

In addition, in order to improve a O2-preserving mechanism, breath-hold divers have developed metabolic adaptations, as a reduced basal metabolic rate (Ferretti, 2001), and increased lactate production, implying a greater reliance on anaerobic metabolism (Andersson et al., 2004). Yet a different adaptive mechanism to hypoxia is the ability to induce growth of new blood vessels into ischemic areas (Lee et al., 2000; Eckle et al., 2008). In particular, a cardioprotective role of hypoxia-inducible factor 1 (HIF-1) levels, essential in the maintenance of oxygen homeostasis (Hagen et al., 2003), has been demonstrated.

Finally, the introduction of meditation techniques such as the Pranayama, performed during the training program, has led to a significant average increase in apnea time. Indeed, meditation enhances the control of breathing (Cysarz and Büssing, 2005) and modulates the emotional processing (Sobolewski et al., 2011), thus contributing to a better air hunger control. The air hunger is considered as a dominant and primitive sensation with a strong emotional component which alerts the organism for a risk to survival (Brannan et al., 2001; Liotti et al., 2001).

Many studies focus on the physiological effects of meditation techniques to gain insight into the physiological requisites responsible for the improvement of health. In particular, the cardiorespiratory interaction seems to play an important role (Bernardi et al., 2001; Cysarz et al., 2004), making use of a slow paced breathing pattern and producing low frequency breathing oscillations. Regardless of the different techniques, meditation causes an alteration in respiration rate, heart rate, metabolic changes and skin conductance (SC) response (Wenger and Bagchi, 1961; Dillbeck and Orme-Johnson, 1987; Jevning et al., 1992; Critchley et al., 2001; Infante et al., 2001; Jerath et al., 2006). All these works have suggested that physiological measures of heart rate variability (HRV), SC response, respiratory amplitude and rate, could be considered as reliable biomarkers of autonomic control during meditation. Meditation was also reported to induce positive variations in immune function (Davidson et al., 2003), and reduced activation of the endocrine system, with subsequent changes in hormonal levels, both during meditation and after (MacLean et al., 1997). Generally, voluntary breathing control may be used to preserve/improve health and seems to be beneficial for recovery processes, reducing anxiety levels, and as part of prevention programs (Bernardi et al., 2001). Hypertensive subjects practicing meditation, for example, show a decrease in systemic and diastolic pressure (Wenneberg et al., 1997).

For these reasons, it is reasonable to hypothesize that the use of meditation, combined with physical training, can modulate psychophysiological responses during BH, extending the duration of BH time and making people less sensitive to air hunger. Moreover, although many studies have investigated the physiological adaptive responses that occur during BH, few have characterized the normoxic phase.

Given this evidence, we made use of experienced breath-hold divers, who practice meditation technique during their training program, to obtain long-lasting electrocardiogram recordings on the one side, and SC responses, as marker of emotional arousal, on the other side. This experimental approach allows us to analyze the cardiovascular time courses during extreme physiological situations, also allowing us to characterize the relative components of the autonomic nervous system corresponding to the regulation of sweat gland function. There is a large volume of published studies describing the role of SC response as a function of the emotional arousal, for a detailed description view (Sequeira et al., 2009).

2. Materials and methods

2.1. Participants and experimental protocol

Data were collected on 14 healthy male diving athletes, mean age 28 (SD 6) years. They had been practicing competitive apnea for at least 2 years without interruption, and meditation techniques during the training program. Inclusion criteria were: nonsmokers, light caffeine, and alcohol consumers, drug free for at least 2 months, and no history of cardiovascular, pulmonary, or neurological/psychiatric diseases.

Each athlete underwent two successive recording sessions, both with closed eyes, while seating in a semi-dark and sound-attenuated room (T = 22°C). The first session was a 5 min basal condition recording, (free breathing condition, FB), the second one a prolonged breath-holding recording (BH condition), whose duration was dependent on the ability of each athlete. Each athlete performed BH in sub-maximal inspiration.

In the following analyses, BH was split in two phases: the normoxic phase of BH (NPBH) and the hypoxic phase of BH (HPBH). The point of separation between NPBH and HPBH was identified by the knee in the SaO2 time course and corresponded to a SaO2 drop rate exceeding 15%/min. Figure 1 shows the time course of signals as a function of BH; the vertical red line indicates the breaking point between NPBH and HPBH.

Figure 1

A typical set of autonomic signals recorded during BH. From top to bottom: electrocardiogram (ECG), respirogram (resp), skin conductance (SC) and oxygen saturation level (SaO2). The red line indicates the breaking point between the normoxic breath-holding phase (NPBH) and the hypoxic breath-holding phase (HPBH).

We based the psychophysiological characterization of subjects throughout the conditions on autonomic-related signals sampled at 500 Hz: DII-lead electrocardiogram (ECG), respirogram (Grass Technologies, US-RI), and SC (Contact Precision Instruments, US-MA) were recorded in parallel with the oxygen saturation (SaO2) level (finger clip pulse oximeter, Nonin Medical, US-MN). The SC signal was recorded on thenar/hypothenar of the non-dominant hand.

All experimental procedures were performed according to the ethical guidelines of the Local Ethical Committees and all athletes gave their informed consent.

2.2. Data analysis

In order to study the autonomic adaptations occurring during BH, a set of features was estimated from the recorded signals. In particular, we considered feature changes occurring both during NPBH and HPBH, separately.

Twenty features were estimated from the signals recorded during FB, NPBH, and HPBH. In particular, from ECG we derived the series of interbeat interval durations (RR series) that constituted the basis for the estimation of features describing the HRV. The autonomic control on cardiac function was characterized by means of HRV features. In the following we provide a definition of these features that we organize into the classes of time-domain measures, frequency-domain measures and non-linear measures. From SC we characterized the spontaneous SC responses (SCRs).

In order to describe the autonomic adaptations, we analyzed the changes of each feature as a function of the conditions (Univariate Statistical Analysis) and the changes of the latent factors model clustering the features in the different conditions (Factor Analysis).

HRV features From the RR series we derived five features in the time-domain, four in the frequency-domain (Camm et al., 1996) and six by using non-linear methods (Pincus, 1991; Kantz and Schreiber, 2004). All the extracted HRV features were signal-length independent. For the time-domain features we considered: mean value of RR intervals (RR), standard deviation of RR intervals (SDNN), standard deviation of successive RR interval differences (SDSD), root mean square of successive RR interval differences (RMSSD), and the coefficient of variation of the RR series (SDNN/RR). Both SDSD and RMSSD are measures of the HRV short-term variability (Brennan et al., 2001), while SDNN and SDNN/RR represent the overall (both short-term and long-term) variability within the RR series (Camm et al., 1996). For the frequency-domain features: power in very low frequency (VLF: 0.003–0.04 Hz), low frequency (LF: 0.04–0.15 Hz), and high frequency (HF: 0.15–0.4 Hz) bands and the ratio between the low and high frequencies bands (LF/HF). The estimation of HRV power within each band was nonparametric since based on the Hamming-windowed FFT.

The parasympathetic and sympathetic nervous activities influence the HRV at different frequency bands; in general the heart rate oscillations in the heart rate band are correlated to the parasympathetic tone, while those in the LF band are correlated both to the parasympathetic and sympathetic activities, depending on the experimental conditions. Finally, the LF/HF ratio was used to analyze the shift of sympathovagal balance toward a sympathetic predominance and/or a reduced vagal tone (Camm et al., 1996).

As non-linear features we considered information-related, dispersion, autosimilarity, and correlation measures. The information content was measured by means of the Shannon entropy (H); the dispersion measures corresponded to the main axes (minor axis—SD1—and major axes—SD2—) in the Poincaré plot; the autosimilarity properties were evaluated by the Approximate entropy (ApEn) and the Correlation Sums (CSs) with embedding dimensions m=2 and m=3; the short-term correlation scaling within each signal was measured by the α coefficient of Detrended Fluctuation Analysis (αDFA).

Nonlinear features of HRV provide additional information about the cardiac dynamics. Both autosimilarity (ApEn and CS) and information-related (H) features were used to quantify the unpredictability and the information creation rate in the RR series. For example, both ApEn and H are dependent on heart failure severity (Beckers et al., 2001). At variance, the α index from the DFA is a measure of the RR autocorrelation and for example it resulted discriminant in patients with acute coronary syndromes: high αDFA values are associated to unresolved cardiac stress (Hsin et al., 2010). The SD1 and SD2 of the Poincaré plot complete the description of RR dispersion and estimate the short and long-term HRV variability, respectively.

SCR features Spontaneous SCRs were characterized by four features estimated in a time-domain according the model of the sudomotor nerve activity underlying SCRs (Alexander et al., 2005). Beside measuring the SCR mean rate (SCRr), this approach enables separating out temporally overlapped SCRs, thus allowing the single SCR measure of response amplitude, (SCRa), of rise latency (SCRl), and of decay time to the baseline (SCRd). The present analysis took into account responses with amplitude equal or greater than 0.05 μ Siemens. The average value of each SCR morphological feature evaluated over all responses during the condition was used to characterize the condition itself.

Both rate and amplitude of SCR provide estimates of the tone of the cutaneous pathway of the sympathetic branch (Nishiyama et al., 2001). In addition, this pathway is tightly modulated by the emotional arousal, since the central nucleus of amygdala projects to the lateral nucleus of hypothalamus, which is one of the structures that regulates the SC activity (Davis and Whalen, 2001).

2.3. Univariate statistical analysis

In order to identify feature changes occurring early and late during BH compared to FB, we performed multiple comparisons between each pair of phases. Since the distribution of some features were not normally distributed (Lilliefors normality test), the comparisons have been carried out via the Kruskal–Wallis one-way test. The level of significance was set at p-value ≤ 0.05. All statistical analyses were performed using the SPSS Statistics 18 (IBM, Armonk, NY, USA).

2.4. Factor analysis

A factor analysis (estimated by a maximum likelihood maximization) was applied to uncover the latent structure underlying the extracted features. Factor analysis models each observed measure as a linear combination of factors, weighed by coefficients (loadings), plus an error term. The eigenvalue spectrum of the correlation matrix gives indications on the number of latent factors of the model: usually, a frequent choice is given by setting the number of factors equal to the number of eigenvalues preceding the scree point in the eigenvalue spectrum (scree test).

Herein, we estimated a model for each condition since we were interested both in identifying the factors describing patterns or responses of latent mechanisms of adaptation, and in highlighting model changes, i.e., significant modifications of the features clustering within the factors. We chose to apply the factor analysis only on the features that were significantly discriminant between conditions and we set the number of factors according to the scree test. Features were standardized in order to deal with their different range of variability. To facilitate factor interpretation on the basis of its loadings, the varimax rotation method has been employed. Using the varimax rotation, each factor tended to have either large or small loadings coefficients. This behavior helps in pointing out the sub-items related to each factor, namely in identifying the groups of inter-related sub-items. The threshold for the absolute values of the loadings was set to 0.7.

3. Results

Fourteen volunteers underwent 5 min of eye-closed relaxation, followed by sub-maximal dry apnea while their ECG, SC, and respirogram were recorded. The subjects were trained athletes and performed the apnea without any preceding hyperventilation procedure. ECG and SC provided two different displays of the autonomic output, while the respirogram was used as a control for the BH.

Despite the low levels of SaO2 they reached at the end of BH (88% ± 5.9%), nobody had any medical and/or neurological complication as a consequence of the BH. On average, the lengths of NPBH and HPBH were 113 s (SD 37 s) and 69 s (SD 31 s), respectively. No arrhythmic events have been detected during BH.

3.1. Univariate statistical analysis

In order to study the autonomic changes from the FB condition throughout BH, we estimated a number of features from ECG and SC signals.

Figure 2 shows the qualitative trends of the features that presented significant changes between conditions. Significant trends between FB and NPBH separate the features in two groups: a rising trend (SCRr, SCRa, αDFA, and CS with m=2 and m=3) and a descending trend (SDSD, RMSSD, HF, ApEn, H, and SD1). From NPBH and HPBH each group of features underwent an additional split. Moreover, LF/HF was the only feature with a significant change from NPBH and HPBH, without a significant one from FB to NPBH.

Figure 2

The qualitative trends of the features between the three conditions (FB, NPBH, and HPBH). The significant changes between conditions are shown as lines, the non-significant changes are not reported. Ascending (descending) lines indicate an increase (decrease) of the related features. For each feature, the order of the y-axis levels indicates how mean values are ordered with respect to different conditions. The FB mean values (FBv) of each measure are considered as the reference levels for the values in the NPBH and the HPBH. The label of each feature has the same color of related trend. All the abbreviations are described in the text.

Table 1 reported the descriptive statistics, separately for the three conditions, and the results of the univariate statistical analysis of the HRV features calculated in the time-domain. On average, the mean RR values showed a slight, not statistically significant, decrease in BH respect to FB. At variance, NPBH induced a significant decrease of SDSD and RMSSD mean values as respect to FB, while the HPBH was accompanied by a significant increase of both measures compared with NPBH. The BH did not have any significant effects on SDNN and SDNN/RR features.

Table 1

Descriptive statistics—mean (SD)—of time-domain HRV features and related Univariate Statistical Analysis results (.

Measures	Descriptive statistics			Condition effect	Pairwise comparisons
	FB	NPBH	HPBH		FB vs. NPBH	FB vs. HPBH	NPBH vs. HPBH
RR	960 (130)	880 (90)	902 (150)	NS	NS	NS	NS
SDNN	80 (20)	70 (30)	103 (50)	NS	NS	NS	NS
SDSD	50 (30)	30 (10)	83 (50)	0.00	0.01	NS	0.00
RMSSD	52 (30)	30 (10)	76 (50)	0.00	0.01	NS	0.00
SDNN/RR	0.080 (0.030)	0.080 (0.030)	0.110 (0.040)	NS	NS	NS	NS

All the time values are expressed in ms. NS stands for not significant at p-value ≥0.05. All the abbreviations are described in the text.

Table 2 reports the descriptive statistics separately for the three conditions, and the results of the univariate statistical analysis of the HRV features calculated in the frequency-domain. Both VLF and LF did not exhibit any between-conditions significant changes; whereas the HF power firstly had a significant drop in the NPBH respect to FB, then a significant increase during HPBH with respect to NPBH. For the LF/HF ratio we did not identify any significant variation between NPBH and FB, while it decreased in HPBH respect to NPBH.

Table 2

Descriptive statistics—mean (SD)—of frequency-domain HRV features and related Univariate Statistical Analysis results (.

Measures	Descriptive statistics			Condition effect	Pairwise comparisons
	FB	NPBH	HPBH		FB vs. NPBH	FB vs. HPBH	NPBH vs. HPBH
VLF	2400 (2100)	1700 (1900)	2500 (5100)	NS	NS	NS	NS
LF	3000 (2300)	1500 (1600)	3300 (3000)	NS	NS	NS	NS
HF	650 (500)	220 (200)	2100 (3400)	0.00	0.01	NS	0.00
LF/HF	5.9 (6.3)	10.2 (9.2)	4.8 (5.9)	0.03	NS	NS	0.01

All the frequency power values are expressed in mV2/Hz. NS stands for not significant at p-value ≥0.05. All the abbreviations are described in the text.

Table 3 completes the description of the HRV features and shows the descriptive statistics, separately for the three conditions, and the results of the univariate statistical analysis of the non-linear HRV features.

Table 3

Descriptive statistics—mean (SD)—of non-linear HRV features and related Univariate Statistical Analysis results (.

Measures	Descriptive statistics			Condition effect	Pairwise comparisons
	FB	NPBH	HPBH		FB vs. NPBH	FB vs. HPBH	NPBH vs. HPBH
ApEn	0.97 (0.16)	0.64 (0.12)	0.59 (0.20)	0.00	0.00	0.00	NS
H	1.36 (0.20)	1.00 (0.19)	1.05 (0.24)	0.00	0.00	0.00	NS
CS (m = 2)	0.04 (0.02)	0.06 (0.01)	0.05 (0.04)	0.00	0.00	0.02	NS
CS (m = 3)	0.02 (0.01)	0.03 (0.01)	0.03 (0.03)	0.00	0.00	0.00	NS
αDFA	0.91 (0.14)	1.17 (0.30)	0.77 (0.14)	0.00	0.01	0.04	0.00
SD1	27.7 (20.8)	15.0 (7.5)	43.9 (32.6)	0.00	0.03	NS	0.00
SD2	63.0 (17.7)	63.0 (24.8)	77.9 (38.6)	NS	NS	NS	NS

NS stands for not significant at p-value ≥0.05. All the abbreviations are described in the text.

Compared to FB, ApEn, and H exhibited a significant decrease in both BH conditions; conversely CS (with m=2 as well as with m=3) significantly increased from FB to both BH conditions. The highest average value of αDFA was estimated during NPBH and it was significantly greater than those estimated during both FB and HPBH. Moreover, the average value of αDFA during HPBH was significantly lower than during FB. Considering the features related to the Poincaré plot, the SD1 was significantly decreased from FB to NPBH and increased in HPBH respect to NPBH, whereas SD2 not exhibited any between-conditions significant variations.

Finally, Table 4 summarizes the descriptive statistics, separately for the three conditions and the results of the univariate statistical analysis of the SC features. HPBH induced significant increments in the rate and the amplitude of SC responses (SCRr and SCRa) compared to FB; the (SCRa) also increased in NPBH (with respect to FB). On the contrary both the time to peak (SCRl) and the decay time SCRd of the individual responses did not change across conditions.

Table 4

Descriptive statistics—mean (SD)—of SC features and related Univariate Statistical Analysis results (.

Measures	Descriptive statistics			Condition effect	Pairwise comparisons
	FB	NPBH	HPBH		FB vs. NPBH	FB vs. HPBH	NPBH vs. HPBH
SCRr	0.10 (0.07)	0.15 (0.09)	0.20 (0.10)	0.01	NS	0.00	NS
SCRl	3.85 (0.80)	4.10 (0.80)	3.30 (0.90)	NS	NS	NS	NS
SCRa	0.02 (0.02)	0.05 (0.05)	0.10 (0.10)	0.01	0.03	0.01	NS
SCRd	5.20 (0.10)	5.30 (0.25)	5.20 (0.20)	NS	NS	NS	NS

The time and amplitude values are expressed in seconds and μ Siemens, respectively. NS stands for not significant at p-value ≥ 0.05. All the abbreviations are described in the text.

3.2. Factor analysis

As suggested by the visual analysis of trends in Figure 2, most of the features have common direction of changes as a function of the conditions, thus suggesting the existence of latent factors underlying the features. To this aim, a Factor Analysis was applied to the features having significant changes between FB and at least one phase of BH. In order to identify the number of latent factors we used the scree test. Figure 3 shows that the scree of the eigenvalues spectra was located at the third eigenvalue for all the conditions. As shown in Figure 4, each factor clustered distinct features, thus representing distinct psychophysiological mechanism that we called: autonomic balance (AB), information increase rate (IIR), and central nervous system driver (CNSD).

Figure 3

The eigenvalue spectra for the three conditions (FB, NPBH, and HPBH). The scree test applied to the eigenvalues spectra was used to estimate the number of factors. In figure the first nine eigenvalues of each spectrum were reported and the scree point, corresponding the third eigenvalue for all the three conditions, is indicated.

Figure 4

Factor analysis loading maps for the three condition (FB, NPBH, and HPBH). For each condition is shown the loading set associated to the three factors (AB, IIR, and CNSD) and related to each feature. The threshold for the absolute values of the loadings was set to 0.7, the loading values below 0.7 was in gray. All the abbreviations are described in the text.

In the FB condition the AB factor was composed by all the short-term HRV features (SDSD, RMSSD, HF, LF/HF, SD1) with the SC measures (SCRr and SCRa) and the Shannon entropy (H). The IIR factor was constituted by the autosymilarity measures of ApEn and CS (m=2 and m=3), and the CNSD factor was composed by variance scaling index αDFA.

It is worth noting that BH did not disrupt the FB factorial structure, and only few features moved among factors. In particular, from FB to NPBH the αDFA shifted from the CNSD factor to the IIR one and both the skin measures (SCRr and SCRa) moved to CNSD factor. The skin measures belonged to the CNSD factor also during the HPBH, whereas H shifted in IIR factor and the αDFA in AB one.

4. Discussion

The present study shows how the autonomic control on heart rate and SC is modulated by means of the psychophysiological conditions occurring along with BH. Now we discuss in some detail the results reported in the previous section, before discussing the conjecture that will emerge from our analysis, namely that trained subjects are able to actively modify some of the driving mechanisms of the central nervous system to the autonomic one.

First of all, as previously reported in (Pingitore et al., 2008), the heart rate did not change during both BH phases, compared to the FB condition. This finding contrasts with many reports that studied the apnea with face or full immersion (Andersson et al., 2004; Heusser et al., 2009). However, the immersion is related to a reflex, induced by a trigeminal stimulation such as the contact of face with cold water, named diving reflex. Indeed, diving reflex is considered to be the major contributor to the initial bradycardia observed in those studies (Andersson et al., 2004). In addition, our divers practiced BH in sub-maximal inspiration, which is not associated with the occurring of bradycardia as observed during the initial phase of the Valsalva maneuver. Thus, the lack of these elements may justify the non-significant decrement of HR (RR intervals increase) across conditions.

In this study many features describing the HRV measure the variability at high frequency, i.e., in the frequency band of the respiratory sinus arrhythmia. These variables are from time (SDSD, RMSSD), frequency (HF) and non-linear (SD1) domains. All these features strongly decreased from FB to NPBH, likely due to the ventilation suspension. Indeed, ventilation produces a modulation of the atrial sinus pacemaker rhythm and of the preload, both coupled with the respiratory activity. In agreement with this line of interpretation, the occurrence of involuntary respiratory movements during the HPBH is accompanied by a general increase of these features compared to NPBH, without reaching the FB levels.

In this experimental protocol no direct information can be derived about the vagal tone, due to the confounding effect of the ventilation. However, insights about the sympathetic tone were derived both from the long-term variability measures of the heart rate (LF, VLF, and SD2) and from the features concerning the SC (SCRr, SCRa, SCRl, and SCRd). Summarizing, the cardiovascular sympathetic pathway as estimated from the long-term components of HRV did not indicate any significant modification. Thus, our study does not highlight sympathetic tone increases as a function of BH. This is an interesting finding since it apparently contrasts with the increase demonstrated in the sympathetic skeletar vascular pathway (Heusser et al., 2009). Both effects are compatible in the perspective of a trophotropic principle: the lack of increase of sympathetic tone in the cardiac pathway is probably aimed at controlling the metabolic need of heart whereas the increase in the vascular pathway is aimed at increasing peripheral resistances, thus diverting the blood flow to the vital organs (Heusser et al., 2010).

As far as the SC features are concerned, they corresponded to the activity of sweat glands which are under the central control via the cutaneous sympathetic pathway. In particular, the rate of discharge of the sudomotor nerve is proportional to the SCR rate, and the amplitude of each SCR is linearly related to the amplitude of the sudomotor burst (Nishiyama et al., 2001). More in details, skin districts linked to electrodermal response are innervated by postganglionic sudomotor neurons located in ganglia of the lateral sympathetic chains (Jänig and Kümmel, 1981). This neurons are influenced by spinal sudomotor neurons in the spinal cord (Cabot, 1990), which are under excitatory or inhibitory input from reticular system in the brain stem. In particular, both reticular activating system and reticular inhibiting system send impulses to electrodermal activity through the descending neural pathways (Sequeira et al., 2009). Besides, hypoxia can induce modifications of the SC features, since the reduction of oxygen supply to skin could influence the sweat secretion (Nishiyama et al., 2001). Some works have described temporal feature variations as markers of particular emotional or cognitive function (Sequeira et al., 2009). Thus, during BH two opposite effects are expected to act: the sympathetic activation and the hypoxia effects. As a result, we identified significant SCR rate and amplitude increases, but we did not detect any significant changes on the temporal features from FB, across the BH phases.

Also, we made use of a set of nonlinear features, as they have proved to efficiently discriminate health and abnormal heart dynamics. This effectiveness has, during the last three decades, corroborated the hypothesis that the autonomic system controlling heart rate is guided by a low-dimensional chaotic dynamics, which, in normal conditions, works in the regime of weak chaos, namely with long-range correlations, hence 1/f power spectrum. Our results on FB confirm in full what has already been reported in the literature (Perkiomaki et al., 2000). The nonlinear features adopted herein have different meaning, looking at different aspects of the nonlinear dynamics under study. We discuss them in some detail in the following, together with our experimental results.

The Shannon Entropy, originally a measure of information, is indeed a measure of variability, as maximal information is accomplished when all the values of the variable of interest are visited with a flat probability. Shannon entropy thus provides information about the shape of the probability density function of the variable of interest (in our case RR spectral frequencies). This function is a priori well defined if the system is stationary and the time sequence is long enough to provide information on the system “invariant distribution,” or in other terms, the system has reached equilibrium. It is however efficient also in discriminating stationarity breakdowns, as the density function develops “holes,” yielding a steep decrease in the Shannon functional. A decrease, however, is also expected when the system becomes more regular, developing a single, major peak in the distribution, in our case, in the spectrum. A steep decrease is in fact apparent when subjects switch from FB to NPBH. Therefore, we are led to suppose that the absence of the stationary mechanical drive of respiratory rhythms may make the system departs from its dynamical equilibrium, in a non-stationary way. Accordingly, the involuntary respiratory movements of HPBH are probably responsible for the partial recovery in the values of H, in a way that is similar to what earlier discussed for the HRV features.

Now we discuss the features stemming from Taken's theorem (Strogatz, 1994), namely the correlations sums CS (m=2, m=3) and the Approximate Entropy ApEn, i.e., the logarithmic ratio of the two CS. These measures assume that a low-dimensional attractor exists with dimension m. The attractor reconstruction can be accomplished via the method of the “embedding dimensions,” namely building a fictitious “phase space” using always the same time course of a scalar function, but with different delays. The celebrated Taken's theorem assures that under general hypotheses of smooth trajectories, the topology of the dynamics in the fictitious (embedding) space is isomorphic to the original one. Although these hypothesis are almost never fulfilled, the analysis proved to be useful in more general cases. In the specific case of CS, one measures the amount of short-term predictability of the trajectory. ApEn measures how the predictability increases with increasing the embedding dimensions, that also correspond to a selection of time steps in the past where the system dynamics were “similar” (below a threshold). In this sense it measures the increase rate of information. Our results show a significant increase in these short-term predictability when passing from FB to BH. This is obviously associated with a decrease on unpredictability, which is normally described in terms of the dynamic entropy increase rate of the system, herein measured through the feature ApEn. This behavior mirrors the earlier described behaviors of H and of HRV features. At variance, we see a significant further decrease of ApEn from NPBH to HPBH. We are forced to induce that involuntary movements play a different role in this case. In fact, these movements tend to be paroxystically regular (periodic), therefore, increasing the predictability of the system.

The Detrended Fluctuation Analysis (DFA) conveys yet different information, since it is basically a measure of scale-free properties in the autocorrelation function of a time series (Goldberger et al., 2002). In this sense it is equivalent to measuring the exponent of a “1/f” power spectrum (Perkiomaki et al., 2000). In HRV literature it is usual to split the DFA analysis into two regimes, a short- and a long-range one. Herein we were forced to only study the short-term, because of the limited duration of our series. However, it must be stressed that “short-term” in DFA is in fact longer than all other measures, as the scale-free behavior was fit to data up to 30 heart beats. The FB indexes were in the region α <1, namely a stationary regime. The reported values are typical of healthy subjects. It has to be noticed that very similar indexes has been reported in basal EEG measurement, if one considers only the signal stemming from global events (abrupt changes in EEG amplitudes or frequency, involving many electrodes). The striking similarity between the two indexes (heart and brain) was already pointed out by (Allegrini et al., 2011), as an example of the driving mechanism taking place in the brain-heart axis. When passing from FB to NPBH, we notice a passage of αDFA from values smaller than one to values higher than one. This is in line with both hypothesis mentioned earlier, namely with a non-stationary drift, due to lack of ventilation and with the presence of a very slow periodic component, possibly related with the low respiration frequency during meditation and maintained also during BH. Also in this case the involuntary diaphragm movements are able to explain how the index changes when switching from NPBH to HPBH, since the regularity of that component of the signal is associated with a very low value of αDFA. The overall index, going to a value not far from 0.5 (the index of white noise) shows that long-range correlations are completely destroyed during hypoxia.

The factor analysis puts the above interpretation of results in a global and clear perspective. The first panel of Figure 4 refers to the FB condition and shows how the different features cluster into three factors that we named AB, IIR, and CNSD. The first factor nicely groups features known to describe the vagal-sympathetic balance. These are all HRV features as well as the SC ones. The second factor refers to entropy increase, and the corresponding features are grouped there. The third factor is represented by αDFA, nicely describing the driving mechanism of the nervous system in a medium-range time. The BH did not disrupt the factorial structure: the factor meanings remain unaltered and only a few of features move from a factor to another. With the passage to NPBH we see an interesting effect, namely that αDFA switched to the second factor and SC features switched to the third. This is probably the most intriguing result of the paper. It suggests that during BH performances, even before the physiological SaO2 break point, trained athletes are able to disconnect heart rate dynamics from the control of nervous system, probably because this disconnection is beneficial for a good performance. However, the SC features, which do not imply a relevant oxygen consumption, are, as expected, driven by central nervous system emotional drive. In the perspective of BH as an emotional laden condition, the separation between HRV and SCRs is at odds with the common findings of an integrated autonomic response (Carrillob et al., 2001). This view is supported by the hypothesis that electro-dermal activity is linked to reticular, limbic, and cortical areas assigned to the control of emotional processing (Sequeira et al., 1995), and it also accords with other studies (Grassi et al., 1998; Critchley et al., 2011).

Finally, during hypoxia, long-range correlations are completely destroyed, and αDFA is not even associated with entropy increase while the SC appears to be the only active emotional output of the increasing urge to breathe sensation.

We conclude stressing the fact that, although we cannot disentangle the effects of prolonged BH, i.e., hypoxia and absence of respiration, from the counterbalance induced by meditation techniques, we provided evidence of significant changes of autonomic features. The factor analysis, in particular, shows changes that (1) cannot be attributed to the mental states in the hypoxic phase, so that elite apnea divers can be taken as a model of general population in this phase; (2) cannot be solely attributed to absence of respiration during the normoxic phase, making the subjects under study peculiar in this phase. This latter fact may likely be due to their particular training, including meditation techniques.

Further studies should be carried out in trained subjects during long-lasting BH exposed to emotional stimulations (i.e., using the International Affective Picture System, IAPS), in order to well describe the effect of mental control during BH on cardiovascular, autonomic, and nervous system.

Conflict of interest statement

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