Effects of Supplemental Oxygen on Cardiovascular and Respiratory Interactions by Extended Partial Directed Coherence in Idiopathic Pulmonary Fibrosis

Idiopathic pulmonary fibrosis (IPF) is a chronic and restrictive disease characterized by fibrosis and inflammatory changes in lung tissue producing a reduction in diffusion capacity and leading to exertional chronic arterial hypoxemia and dyspnea. Furthermore, clinically, supplemental oxygen (SupplO2) has been prescribed to IPF patients to improve symptoms. However, the evidence about the benefits or disadvantages of oxygen supplementation is not conclusive. In addition, the impact of SupplO2 on the autonomic nervous system (ANS) regulation in respiratory diseases needs to be evaluated. In this study the interactions between cardiovascular and respiratory systems in IPF patients, during ambient air (AA) and SupplO2 breathing, are compared to those from a matched healthy group. Interactions were estimated by time series of successive beat-to-beat intervals (BBI), respiratory amplitude (RESP) at BBI onset, arterial systolic (SYS) and diastolic (DIA) blood pressures. The paper explores the Granger causality (GC) between systems in the frequency domain by the extended partial directed coherence (ePDC), considering instantaneous effects. Also, traditional linear and nonlinear markers as power in low (LF) and high frequency (HF) bands, symbolic dynamic indices as well as arterial baroreflex, were calculated. The results showed that for IPF during AA phase: 1) mean BBI and power of BBI-HF band, as well as mean respiratory frequency were significantly lower (p < 0.05) and higher (p < 0.001), respectively, indicating a strong sympathetic influence, and 2) the RESP → SYS interaction was characterized by Mayer waves and diminished RESP → BBI, i.e., decreased respiratory sinus arrhythmia. In contrast, during short-term SupplO2 phase: 1) oxygen might produce a negative influence on the systolic blood pressure variability, 2) the arterial baroreflex reduced significantly (p < 0.01) and 3) reduction of RSA reflected by RESP → BBI with simultaneous increase of Traube-Hering waves in RESP → SYS (p < 0.001), reflected increased sympathetic modulation to the vessels. The results gathered in this study may be helpful in the management of the administration of SupplO2.


INTRODUCTION
Idiopathic pulmonary fibrosis (IPF) is a chronic, restrictive, and progressive disease characterized by fibrosis and inflammatory changes in lung tissue. IPF also affects lung vasculature and prevents an adequate gas exchange, thus leading to exertional chronic arterial hypoxemia and dyspnea. The disease has a bad prognosis since after diagnosis the 50% of the patients die within 3-5 years (King et al., 2019). IPF etiology, diagnosis, treatment, and influences on quality of life, among other, have been investigated, but the focus has always been lungs performance. There is evidence that IPF patients also manifest comorbidities such as cardiovascular diseases, lung cancer, or pulmonary hypertension (Caminati et al., 2019). Furthermore, arterial hypertension has been related to dysfunctional autonomic cardiovascular control (Mancia and Grassi, 2014) and considerable impact on IPF disease progression and patient mortality (Buendía-Roldán et al., 2017).
Recently, IPF has been hypothesized as a systemic disease, but its influence on the autonomic nervous system (ANS) regulation has not been assessed as in other pulmonary disorders. For example, chronic obstructive pulmonary disease (COPD) negatively affects the cardiovascular system and the ANS regulation; the autonomic dysfunction is an important factor in the underlying pathophysiological mechanisms of the disease. The former conclusion has been mainly stated based on the analysis of heart rate variability (Van Gestel and Steier, 2010;Mohammed et al., 2015). However, for IPF disease studies about ANS regulation are scarce. Furthermore, supplemental oxygen (SupplO 2 ) has been prescribed to IPF patients to improve clinical symptoms, but its impact on ANS regulation of cardiovascular and respiratory systems has not been evaluated in these patients.
The effects of oxygen supplementation have been reviewed on healthy subjects and patients with cardiovascular diseases, but on respiratory patients just for COPD. Until now the evidence about the benefits or drawbacks of oxygen supplementation is not conclusive. On the one hand, some authors point out that this type of clinical intervention does not help to increase oxygen delivery, i.e., the rate at which oxygen is transported from lung to microcirculation depending on cardiac output and arterial oxygen contents (Smit et al., 2018a). Furthermore, it is plausible that in critically ill patients, it could be associated with increased hospital mortality (You et al., 2018). Conversely, oxygen therapy in COPD patients is clinically well-accepted and it is recommended to enhance exercise capacity (Stoller et al., 2010). However, some trials have not found enough evidence to support that long-term oxygen therapy improves COPD patients' mortality rate (Khor et al., 2019). It is important to point out that few studies have addressed the efficiency of oxygen administration in COPD and Interstitial Lung Disease (ILD) patients (Khor et al., 2019) and that studies about the ANS regulation in COPD patients were based on classic heart rate variability (HRV) parameters (Mohammed et al., 2017). Regarding IPF disorder, the authors of the present contribution compared the hemodynamic response to SupplO 2 between ill and healthy subjects, showing potential detrimental effects of SupplO 2 on IPF hemodynamics, particularly on total peripheral resistance (TPR) and cardiac output (Santiago-Fuentes et al., 2021). Different methods have been employed to assess ANS regulation by the analysis of heart rate variability. In the case of HRV, it reflects the variations in the beat-to-beat interval and the corresponding variability time series is built up by diverse oscillatory modes. To extract the information from the time series, linear and nonlinear schemes have been proposed. For instance, the short-term HRV spectral density representation has been broadly used to analyze the sympathetic and parasympathetic modulation by the power in the low frequency (LF) and high frequency (HF) bands, respectively. Another way to assess ANS regulation of cardiovascular system is the nonlinear approach that allows to incorporate the analysis of complexity of underlying physiological mechanisms. Nonlinear indices by symbolic dynamics (SD) and detrended fluctuation analysis (DFA) have enabled the exploration of the cardiovascular system adaptability, among other, for example in elderly population (Beckers et al., 2006;Voss et al., 2015). Furthermore, the analysis of physiological variability time series has advanced from the univariate to bivariate and, to multivariate type to discover the complex interactions between human body subsystems. Nowadays, an open research area is to assess the cause-effect relationship to elucidate the complex picture of the autonomic control of the cardiovascular system and the complex interplay between cardiovascular and the respiratory systems.
Recently, the use of multivariate autoregressive models and the assessment of the directional interactions among a set of physiological variables (i.e., the so-called Granger causality) has been proposed for the analysis of ANS regulation under pathological and non-pathological conditions Charleston-Villalobos et al., 2019). Particularly, the extended partial directed coherence (ePDC) has gained interest as a tool to estimate the causality in the frequency domain in presence of instantaneous interactions , i.e., as the effect from respiration to systolic blood pressure and beat to beat interval as well as the systolic blood pressure to beat to beat interval. Consequently, this study aimed to analysis cardiovascular and respiratory times series of variability by linear and nonlinear indices as well as Granger causalities (GC), via a multivariate autoregressive model including instantaneous effects, in IPF patients in comparison with healthy subjects under the effect of short-term SupplO 2 .

Subjects, Acquisition Protocol, and Preprocessing
This study includes 19 (8 women and 11 men) healthy subjects (CON) and 20 IPF patients (9 women and 11 men) with 67.79 ± 5.00 and 65.8 ± 6.48 years old, respectively. All subjects were medically evaluated at the National Institute of Respiratory Diseases in Mexico City after they accepted the invitation to participate and signed an informed consent according to the Declaration of Helsinki. Table 1 depicts different parameters related to clinical measures and respiratory functional tests of both groups. Signals acquisition was performed via a Biopac MP150 system during morning hours including ECG, continuous noninvasive arterial blood pressure, and peripheral blood oxygen saturation. Also, a thoracic belt was used to acquire the respiratory signal and all signals were sampled at 1,000 Hz. Raw signals were acquired in supine position continuously during 10 min with the subjects breathing spontaneously ambient air (AA) and an additional 10 min breathing SupplO 2 at 3 L/min to ensure an arterial oxygen saturation above 94% (Santiago-Fuentes et al., 2021). Time series of successive beat-tobeat intervals (BBI), respiratory amplitude (RESP) at BBI onset as well as systolic (SYS) and diastolic (DIA) blood pressure (BP) were extracted from recorded signals; all extracted time series were manually reviewed and corrected. For GC analysis, the variability time series were resampled at 2 Hz using spline interpolation and normalized to zero mean and unit variance. For dynamic data analysis, consecutive windows of 5 min shifted by 30 s (90% overlap) were used. Therefore, the influence of SupplO 2 on IPF was studied using 31 windows including three phases labelled as ambient-air (AA), transition phase (TPH) and steady supplemental oxygen (SupplO 2 ), as indicated in Figure 1. The transition phase (TPH) is characterized by windows sharing AA and starting SupplO 2 conditions. In each window, univariate, and bivariate indices, as well as the ePDC were estimated.

Linear Analysis
Time and frequency domain linear indices were extracted as the mean value, the root mean square of successive differences (rmssd), and power in the very low-frequency range (VLF, 0.003-0.04 Hz), low (LF, 0.04-0.15 Hz), and high frequency (HF, frequencies >0.15 Hz) bands in agreement with the standardization proposed by the Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology (AuthorAnonymous, 1996). The index rmssd is considered the more precise marker of chronotropic cardiac vagal influence (Minarini, 2020), for example, it provided a good differentiation between healthy subjects and vasovagal syncope patients for systolic and diastolic blood pressure variability (Reulecke et al., 2016). For estimating the power of LF and HF bands a parametric autoregressive model was used in each temporal window along of the variability time series. The model parameters and optimal order were estimated using the Burg method and the Akaike Information Criterion (AIC), respectively.

Nonlinear Analysis
To explore the nonlinear time series properties two techniques were applied: Detrended Fluctuation Analysis (DFA) to quantify the fractal scaling properties of a time series, and Symbolic Dynamics (SD), which is a coarse-grained method based on symbols (Reulecke et al., 2016). For DFA analysis, the short-term fractal scaling exponent (α 1 ) was calculated over equal and non-overlapping segments with length between 4 and 16 beats while the long-term exponent (α 2 ) for segments with length of 16-64 beats (Peng et al., 1995). Furthermore, by SD technique, nonlinear indices related to the low or high variability of the signal were determined. An alphabet, consisting only of symbols "0" or "1", was used to create words of length six and for BBI time series, thresholds of 2, 5, 10 and 20 ms were established (Voss et al., 2015). In this sense, a low variability index associated with the word "000000" and a high variability index associated to the word "111111", were counted (Schulz and Voss, 2017). In the case of blood pressure variability (BPV), the 6-length words were created by selecting thresholds of 1, 2, 3 and 4 mmHg (plvar 3 and 4) (Reulecke et al., 2016).

Bivariate Analysis
The Dual Sequence Method (DSM) is a linear technique to estimate the arterial baroreflex sensitivity (BRS) by the analysis of spontaneous fluctuations in systolic BP and BBI time series. Traditionally, a pattern of three consecutive increments or decrements in SYS and BBI are tagged as bradycardic (bslope, increase in SYS that causes an increase in BBI) or tachycardic (tslope, a decrease in SYS that causes a decrease in BBI) sequences, respectively. However, we avoided using the classical thresholds for SYS and BBI due to the criticisms about them (Gouveia et al., 2007) and in consequence, different pattern length was used including one up to three samples (N1-N3). Furthermore, the synchronous responses in the same beat interval (T0) and delayed BBI responses shifted by one up to three beats (T1-T3) were also explored, as suggested by (Malberg et al., 2002;Gouveia et al., 2007).

Multivariate Autoregressive Modeling, Extended Partial Coherence and Time-Frequency Analysis of Interactions
Establishing the cause and effect, or the driver-response relationship, between physiological systems has been of great interest in diverse biomedical applications. Fundamental to the driver-response relationship is the concept of Granger causality which states that if a signal improves the prediction of a second signal, above and beyond, its prediction in terms of its own past, then the first signal causes the second one. GC has been formalized in terms of multivariate time series analysis and its treatment in the frequency domain leading to the concept of partial directed coherence. Particularly, the extended partial where the model coefficients B(k), k = 0, . . . ,q, are related to instantaneous and strictly causal effects, q is the model order and W(n) is the innovation process formed by white and uncorrelated noises with diagonal covariance matrix Λ = diag (λ i 2 ). The identification of the extended MVAR model can be achieved from a strictly causal MVAR model with coefficient matrix Charleston-Villalobos et al., 2019). Furthermore, to calculate B(0) is necessary to perform a Cholesky decomposition of the input covariance matrix Σ LΛL T , with Σ cov(U(n)), of the strictly causal model to obtain the matrix L (I − B(0)) −1 . From this decomposition L is a lower triangular matrix as well as B(0) with null diagonal, but it is required to order the times series in Y(n) to set the direction, but not the strength, of physiological influence among them. Consequently, in the present study, to achieve the former constrain, two models were proposed: i) MVAR model 1 (MVAR 1 ) where Y(n) was built up as y 1 = RESP, y 2 = SYS and y 3 = BBI and ii) MVAR model 2 (MVAR 2 ) with y 1 = RESP, y 2 = DIA, and y 3 = BBI. The model order q was obtained by the minimum of the function AIC(q) Nlog(detΣ) + M 2 q, where Σ is the input covariance matrix of the strictly causal model. From the frequency domain representation of Eq. 1 is possible to obtain the spectral density of Y(f) and its inverse, from which the partial coherence (PC) function Π ij (f) between y i and y j is defined. However, PC cannot provide information about causality due to its symmetrical nature and, consequently, a factorization of Π ij (f) is necessary to produce ePDC, named χ ij (f), as: with T equal to the sampling period and λ i is a diagonal value of the noise covariance matrix of the W process in Eq. 1. Therefore, ePDC is a directional frequency domain measure of connectivity quantifying the influence of the process y j on the process y i , removing the influence of other processes. The ePDC is normalized with respect to the structure that sends the signal taking values between 0 (absence of causal coupling) and 1 (full causal coupling) at frequency f. In this study, the identification of eMVAR model was performed by the standard least-squares method and the estimated ePDC was corrected by statistical hypothesis testing based on setting a threshold for significance using causal Fourier transform surrogates Charleston-Villalobos et al., 2019). In a few words, once the eMVAR model is estimated, if there is a direct causality from the process y j to y i , the corresponding coefficients b ij (k, w l ), k 0, . . . , q, at each temporal window (w l ) are set to zero, and from there the surrogates are obtained; causality is selectively destroyed only over the direction under study, details of the procedure can be found in  To accept or reject the estimated ePDC, its magnitude at each frequency was compared with the threshold obtained from 100 surrogates (95% confidence interval). If the magnitude of ePDC was below the threshold, the null hypothesis was accepted, reflecting the absence of interaction, and the estimated value was set to zero. In the present study, the time series of variability were segmented producing a time-frequency representation of interactions (TFR i ). In TFRi, the x-axis is related to the 31 temporal windows of 5 min and the y-axis depicts the frequency in Hz (0.0-0.50 Hz) using 512 bins. Furthermore, the z-axis represents the magnitude of the ePDC (0-0.4) by a color palette from blue to red, respectively. To count with the same number of samples in each window a resampling was performed at 2 Hz. It is worthy to note that, at least for this research protocol, the TFRi without or with resampling are close to each other and the dynamic changes produced by IPF, and supplemental oxygen remain.

Statistical Analysis
To evaluate the effect of supplemental oxygen, in this study the statistical analysis was performed with two statistical test, withingroups and between-groups. The within-groups test was carried out on univariate and bivariate indices values comparing window 5 vs. windows 12 to 31 (TPH and steady SupplO 2 ) by the Wilcoxon Sign Rank-sum. All p-values were corrected using the Benjamini and Hochberg correction (Benjamini and Hochberg, 1995). In the case of between-groups comparison, the analysis was carried out on the univariate and bivariate indices values, and ePDC magnitude at each tile, in each temporal window by the nonparametric Mann-Whitney-Utest. In both statistical analyses, the significance was set at three levels for descriptive purposes: slightly significant for p < 0.05 (green), moderately significant for p < 0.01 (yellow), and highly significant for p < 0.001 (orange).

RESULTS AND DISCUSSION
The section is organized presenting first the results for the AA phase followed by the findings in the TPH and steady SupplO 2 phases. The SupplO 2 effect on healthy subjects and IPF patients is obtained looking for the time course of linear and nonlinear indices as well as by the dynamics of cardiovascular and respiratory interactions as compared with those in the AA phase. Statistical analysis between-groups, throughout all the phases, is displayed at the bottom of the graphs while withingroup, comparing window 5 vs. windows 12 (start of TPH) to 31, is depicted at the top. For saving the space of article, selected indices are included in Figure 1, while others are only discussed. Also, the indices in Figure 1 were ordered as they are used throughout the manuscript.

Clinical and Functional Measurements
According to Table 1, the IPF group was characterized by moderate reduction in vital capacity and FEV1, mild reduction in diffusion capacity, hypoxemia, mild hypercapnia, and high respiratory rate. In contrast, the CON group had no signs of pulmonary alterations; although the PaCO 2 and PaO 2 were not registered at the time of the study, it is plausible to estimate them at the altitude of Mexico City (Vazquez-Garcia and Perez-Padilla, 2000). Consequently, the estimated normal values of PaCO 2 and PaO 2 for Mexico City residents are around 32.7 and 65.9 mmHg, respectively; there were no statistically significant differences in age, anthropometric measures, hematocrit, and hemoglobin values. As can be seen in Figure 1A, peripheral oxygen saturation (SpO 2 ) increased importantly from TPH towards the steady SupplO 2 phase for both groups until a saturation of 96% was reached, but no statistically significant differences between-groups were found during SupplO 2 . In contrast, within-group statistical analysis revealed highly significant differences. Furthermore, Rico et al. observed that S P O 2 tends to decline during the aging process that is in line with the SpO 2 in our Control group (Rico et al., 2001). On the other hand, the IPF group showed a greater dispersion of SpO 2 during AA phase that may be explained by age, disease stage, and pulmonary condition. Also, it is relevant to point out that the population of this age was under medication for other comorbidities. A comprehensive discussion of this aspect for the groups under study can be found in (Santiago-Fuentes et al., 2021).

Linear and Nonlinear Univariate Analysis in AA
For mean BBI index, statistically significant differences between CON and IPF groups were found from windows 8 to 12, Figure 1B, while for the mean respiratory frequency highly significant differences were found for the whole phase, Figure 1C. Also, the mean SYS and DIA BP for CON tended to be higher than for IPF, Figure 1D. For linear indices as rmssd, there were no significant differences for BBI, SYS or DIA, Figures 1E,F, but for BBI the IPF group showed a tendency to lower values. Also, for BBI power in the LF and HF bands was significantly different (p < 0.05); particularly, the CON group showed higher BBI-HF power and consequently, higher cardiac vagal influence than IPF, Figure 1G. For SYS-nLF, IPF showed a tendency to higher values, i.e., increased sympathetic influence on the vasculature than in CON. In the case of SD nonlinear analysis, for IPF the indices BBI-phvar2 as well as SYS-plvar3 ( Figure 1H) showed a tendency to lower probability values than CON, i.e., low BBI variability and higher SYS variability, respectively. Nevertheless, in the DFA analysis, BBI-α 2 was statistically different (p < 0.05) between groups, where IPF patients showed higher sympathetic activity than CON, Figure 1I. Consequently, by univariate analysis, statistically significant differences between groups were found just for BBI, indicating a significant cardiac sympathetic modulation in IPF patients during AA scenario.

Linear and Nonlinear Univariate Analysis in TPH and Steady SupplO 2
At the firsts windows in TPH, BBI showed statistically significant differences between-groups for mean, BBI-HF, and BBI-α 2 indices, Figure 1. Specifically, the mean cardiac frequency decreased in both groups, slowly for CON and more drastically for IPF, whereas for HF and α 2 , the CON group showed higher and lower values than IPF, respectively. The former behavior is in line with the expected effect of SupplO 2 on the cardiac frequency. The mean respiratory frequency was reduced for both groups however, the between groups differences were kept along both phases, Figure 1C. Although betweengroups no statistically significant differences were found for mean SYS and DIA in TPH or steady SupplO 2 , the mean SYS showed a tendency to increase in CON while the mean DIA decreased in IPF.
To evaluate the influence of O 2 within each group, statistical comparisons were achieved between AA and TPH or steady SupplO 2 phases. The effect of O 2 was significantly different for CON and IPF. Specifically, significant differences were obtained for BBI, mean respiratory frequency, and DIA. For the IPF group, the BBI-rmssd index showed highly significant differences (p < 0.001) from windows 16 to 31, i.e., the cardiac vagal influence in IPF was increased to a greater extent by O 2 , Figure 1E. The former time course was supported by the nonlinear BBI-phvar2 index (p < 0.001). Furthermore, BBI-LF and BBI-α 1 indices provided significant within-group differences. Regarding BBI-LF index, only for the IPF group, the LF power increased from TPH and throughout steady SupplO 2 phase (p < 0.001) while the BBI-α 1 index decreased. It is worthy to note that for SYS, in none of the groups neither linear nor nonlinear indices provided significant differences. Although, for IPF, SYS-rmssd had a tendency to increase while SYS-plvar3 to decrease, indicating that SYS increased its variability with O 2 , Figures 1F,H. In the case of DIA, just for IPF, mean DIA decreased due to oxygen from TPH to the end of the steady SupplO 2 phase in a moderately significant way (p < 0.01), that may be related to the lower cardiac frequency with O 2 , i.e., a vagal modulation effect, Figure 1D. Also, the DIA-plvar4 index showed moderately and highly significant differences along the two phases, the percentage of words of low variability was reduced, i.e., in the IPF group the DIA variability was increased by SuppO 2 , Figure 1J.

Bivariate Analysis in AA, TPH and Steady SupplO 2
For DSM analysis, different pattern lengths and delays were tested, and the results indicated that statistical differences between groups occurred with lengths of 2-3 samples and shifts between 0 and 3 beats. The bivariate indices during AA, associated with arterial baroreflex sensitivity, did not show statistical differences between groups, only a tendency to lower values for IPF. The former behavior is in line with the literature indicating that hypoxia produces a resetting of arterial baroreflex, without changing the sensitivity, to higher heart rates and systolic blood pressures due to stimulation of peripheral chemoreceptors (Halliwill et al., 2003). In contrast, the IPF patients of the present study were characterized by significantly higher cardiac frequency and cardiac sympathetic modulation but similar systolic pressure. It is worthy to mention that the interaction between baroreflex and peripheral chemoreflexes remains controversial (Kronsbein et al., 2020). During steady SupplO 2 , from window 21, bslope and tslope indices provided significant differences (p < 0.05 and p < 0.01) between groups, Figures 1K, 1L. For both indices, IPF showed lower values than CON, i.e., IPF decreased their baroreflex in comparison to CON. In the case of withingroups analysis, in TPH and steady SupplO 2 phases, bslope tends to increase for CON while for the IPF group tends to decrease. Systolic pressure, being a variable to control, tends to oscillate with greater amplitude when feedback mechanisms are deficient, showing an inverse relationship with variations in heart rate, as a variable to regulate (Mancia et al., 1986;Lanfranchi and Somers, 2002). In fact, it should be noted that the significant reduction in cardiac output in IPF observed during SupplO 2 (Santiago-Fuentes et al., 2021) could be a consequence of alteration of the cardiovascular regulation. The former behavior may be explained, on the one hand, by statistically significant changes (attenuation) of the sensitivity of the baroreceptors reflected by bslope and tslope indices ( Figures  1K,L), increase in BBI ( Figure 1B), and in total peripheral resistance and, on the other hand, by statistically nonsignificant trends in systolic blood pressure and its variability. Furthermore, diverse research pointed out that clinicians should be aware of the prognostic implications of increased blood pressure variability, a marker of cardiovascular decompensation, which may lead for example to organ damage (Hoecht, 2013).

TFR i Magnitude Distribution in AA
An averaged time-frequency representation of RESP → SYS interaction by MVAR 1 , associated with the information flow from RESP to SYS, is showed in Figure 2A, for CON (above) and IPF group (below). The TFR i magnitude distribution points out differences between groups for the LF and HF bands, the significant differences are displayed in Figure 2B. The power in the LF band has been associated with the sympathetically mediated BP vasomotor modulation, the so-called Mayer waves around 0.1 Hz. The origin of Mayer waves has not been elucidated, but some authors agreed that these waves could be related to an oscillatory sympathetic activation that, in the specific case of humans, is independent of factors such as gender, age, or posture and can be induced in BP by hypoxia (Julien, 2006;Ghali and Ghali1, 2020). According to the RESP → SYS interaction in the LF band, the TFRi magnitude for IPF is higher than for CON that is in line with the SYS univariate results presented above; the IPF group showed a tendency to an augmented sympathetic activity in AA, as revealed by SYS-rmssd ( Figure 1F). Furthermore, for the CON group the magnitude of the TFR i in the HF band is concentrated around 0.28 Hz while for IPF group is spread out around 0.42 Hz, i.e., around the respective mean respiratory frequency, Figure 1C. It is plausible that different mechanical effect on central blood volumes is produced by higher respiratory frequency and lower FVC in IPF than in CON. Also, in the HF band Traube-Hering waves (THW) have been reported associated with respiratory-related fluctuations in sympathetic outflow that promotes changes in vascular tone (Menuet et al., 2020). In fact, the IPF group showed a significant increased TPR during AA, i.e., increased vascular tone (Santiago-Fuentes et al., 2021). In the case of SYS → BBI interaction there were no statistically significant differences and consequently, the corresponding TFR i is not shown. A possible explanation could be associated to the fact that subjects were in supine position and then, baroreflex feedback control may be blunted. The RESP → BBI interaction, Figure 2C, shows that the magnitude for CON is spread out over the whole HF band and is higher than for IPF that is localized around the mean respiratory frequency, the corresponding statistical differences are displayed in Figure 2D. Therefore, the magnitude of the RESP → BBI interaction pointed out that respiratory sinus arrhythmia (RSA), associated with the cardiac vagal influence, occurs at different operating point in CON and IPF. The former result is in line with other studies which suggest that RSA decreases with hypoxia (Yasuma et al., 2001) and in fact, the IPF group of the present study was characterized by hypoxemia (Table 1). Furthermore, RSA was more relevant for CON as the magnitude of the RESP → BBI interaction is higher than the corresponding of the RESP → SYS interaction, i.e., the influence of RESP on the heart is higher than on the vascular subsystem. Regarding the BBI → SYS interaction, Figure 2E, associated with the mechanical feedforward influence, a relevant magnitude is located towards the VLF and LF bands for both groups; however, there were no statistically significant differences. Results by MVAR 2 , using DIA instead of SYS are shown in Figure 3. The TFR i of the RESP → DIA interaction is shown in Figure 3A and resembles the TFR i of RESP → SYS, but for both groups the magnitude spread further in the LF and HF bands, statistically significant differences are shown in Figure 3B. For the DIA → BBI interaction, Figure 3C, the TFRi displayed a magnitude with a trend to be higher in CON, particularly in the upper part of the HF band, while for the IPF group a more uniform distributed magnitude is observed across the frequency bands. Furthermore, for both groups the DIA → BBI interaction presents a higher magnitude than SYS → BBI interaction. The former behavior may be interpreted in terms of an impaired left ventricular (LV) diastolic filling in IPF in contrast to a preserved LV systolic function, as was shown in a previous study (Papadopoulos et al., 2008). Also, the lower magnitude in the DIA → BBI interaction for the IPF group may reflect peripheral vascular changes due to increased arterial stiffness and in conjunction with a high heart rate may be indicative of the prevalence of sympathetic tone. For the BBI → DIA interaction, the TFR i magnitude distribution resembles the corresponding to BBI → SYS of MVAR 1 , Figure 3E, but with a tendency to display lower values for both groups. Furthermore, a consistent activity in the LF band was found, the Mayer wave. For the DIA → RESP interaction, the corresponding TFR i did not show relevant magnitude in any phase and, consequently it was not included in the paper. Therefore during the AA phase, based on the TFR i magnitude, its spread and relevance of statistically significant differences between groups, IPF during the AA phase mainly impacted the RESP → SYS and RESP → BBI interactions.

TFR i Magnitude Distribution in TPH and Steady SupplO 2
In the case of RESP → SYS interaction, during TPH the SupplO 2 increased the TFR i magnitude in both groups, Figure 2A. For the CON group the magnitude spread further in the LF and HF (the Traube-Hering waves) bands as compared with the AA phase. In contrast, for the IPF group the magnitude was increased mainly in the upper part of the HF band, the statistically significant differences are shown in Figure 2B. Furthermore, in the steady SupplO 2 phase, the TFR i magnitude decreased in both groups but remains higher in IPF, it seems that the interaction for both groups reached an operation level like the level in AA phase. Presumably, in IPF the heightened arterial peripheral chemoreceptors activity/sensitivity by the hypoxia keeps the drive of RESP on SYS (Stickland et al., 2016). In fact, chronic hypoxia produces among other things, hyperplasia of the carotid body and increases its activity. As can be seen in Figure 1C, throughout the SupplO 2 phase the mean respiratory frequency decreased in the IPF group, with a moderate p-value with respect to AA, but the influence on SYS does not change significantly.
During TPH, the RESP → BBI interaction reveals that RSA was highly significant different between CON and IPF groups, Figures 2C,D, i.e., RSA in CON was stronger. The former result is in line with the fact that vagal modulation facilitates RSA at rest while sympathetic activation attenuates its magnitude. Moreover, for the RSA physiological purpose, there are three major hypotheses, i.e., 1) improvement of gas exchange, 2) minimization of the energy consumption for the heart, and 3) reduction of BP fluctuations (Buchner, 2019). One of the premises of hypothesis 3) is that suppression of RSA generally increases BP fluctuations in its HF band. In fact, the former hypothesis may explain the decrease in RSA for the IPF group that may be related to the tendency to increase of SYS-rmssd and SYS-phvar1 indices, through SupplO 2 as discussed above. During the steady SupplO 2 phase the RESP → BBI interaction showed highly significant differences. For CON, the TFR i significant magnitude in the HF band, almost in the whole phase, reveals that SupplO 2 affects importantly respiratory drive on BBI. In contrast, for the IPF group the RESP → BBI interaction remains at similar operation level as previous phases, in consequence, the RSA in IPF was lower than in CON. It is noteworthy that in SupplO 2 phase, for the CON group the RESP → BBI (RSA) interaction increased while RESP → SYS magnitude (THW) decreased (Barnett et al., 2020). In the case of BBI → SYS, for IPF during TPH, the TFR i magnitude showed a sustained level of driving indicating a higher influence of the cardiac mechanical Frontiers in Network Physiology | www.frontiersin.org March 2022 | Volume 2 | Article 834056 modulation, Figure 2E. In steady SupplO 2 , in the LF band there were no significant differences, but there was a trend in IPF towards a greater magnitude, which is associated with sympathetic modulation on the vasculature. By MVAR 2 using DIA instead of SYS, the RESP → DIA interaction reveals that the TFRi magnitude for CON and IPF is concentrated around the mean respiratory frequency of each group as occurs for RESP → SYS, Figures 3A,B. It seems that the oxygen delivery did not influence significantly the information flow from RESP to DIA. Regarding the DIA → BBI interaction, Figure 3C, the TFR i magnitude in steady SupplO 2 phase showed some differences between groups in the HF band, and mainly toward the end of the phase, statistically significant differences are depicted in Figure 3D. A plausible explanation is that the diastolic BP is related to the pumping activity of the aorta and large arteries of the circulating system. So factors such as systemic vascular resistance or arterial stiffness modify diastolic BP more than systolic BP, and therefore, the detriment of baroreflex activity at rest, as reflected by the tslope index in Figure 1L, is more evident in DIA → BBI than in SYS → BBI. In the case of the BBI → DIA interaction, Figure 3E, the TFR i resembles the behavior of BBI → SYS interaction in the LF band and it seems to be slightly affected by oxygen but there were no statistically significant differences, Figure 3F. It is worthy to note that stationarity is assumed for temporal windows of 5 minutes of cardiovascular and respiratory time-series. However, during TPH nonstationary behavior may be elicited and consequently, a modified time-frequency approach could be applied.

CONCLUSION
According to the results of the present study, the IPF group in the AA phase, as compared with a matched healthy group, was characterized by high sympathetic modulation to the heart (supported by BBI-rmssd, BBI-α 2 indices, and increased mean respiratory frequency, among other indices) with a decreased influence of the parasympathetic activity. It is well-known that chemoreceptors sense the partial pressure of oxygen in blood vessels, and they are important modulators of sympathetic activation in response to hypoxemia. The activation is also known as chemoreflex-mediated sympathetic activation and one of the consequences is the hyperventilation that leads to the inhibition of the chemoreflex due to the stretch mechanoreceptors activity in the thoracic cage (Kara et al., 2003). Also, the IPF patients during SupplO 2 had blunted baroreflex (confirmed by DSM results and almost inexistent SYS → BBI interaction), as well as RSA activity (confirmed by BBI-HF and RESP → BBI interaction). The former autonomic behavior could be explained partially by the chronic hypoxemia and hypercapnia in the IPF group. In this study, the patient group suffered mild hypoxia and some evidence of hypercapnia, so the enhanced sympathetic activity and reduced parasympathetic response may be a consequence of alterations in baroreceptors and peripheral chemoreceptors, which reduce the response to supplemental oxygen. In fact, Van Gestel and Steier showed evidence related to the high sympathetic activity in COPD patients due to the impaired baroreflex sensitivity altered by hypoxia but not by hypercapnia (Van Gestel and Steier, 2010). Additionally, an impairment of cardiac control and sympathetic overactivity may be related to age, but these responses seem to be higher for the IPF group. In fact, Porta et al. at supine rest phase, found a possible age-related impairment of the cardiac control and altered response to stressors in conjunction with a gradual decrease in SYS complexity and thus, an increase of the sympathetic activity (Porta et al., 2014). The effect of oxygen has been studied mainly by its hemodynamic effects and the analysis of heart rate variability (Lund et al., 1999;Gole et al., 2011;Smit et al., 2018a;Smit et al., 2018b). To the best of our knowledge, this is the first study that tackle the effect of oxygen supplementation in IPF with a new perspective from the point of view of linear and non-linear indices as well as the dynamics of cardiovascular and respiratory systems interactions. For IPF patients the results showed that during AA phase: 1) the mean BBI value and power of BBI-HF band, as well as the mean respiratory frequency were significantly lower (p < 0.05) and higher (p < 0.001), respectively, indicating a strong sympathetic influence, and 2) the RESP → SYS interaction was characterized by Mayer waves and diminished RESP → BBI, i.e., decreased respiratory sinus arrhythmia. In contrast, for IPF during short-term SupplO 2 phase: 1) oxygen might produce a negative influence on the systolic blood pressure variability (SYS-rmssd was increased among other indices), 2) the arterial baroreflex reduced significantly (p < 0.01), and 3) reduction of RSA (RESP → BBI) with simultaneous increase of Traube-Hering waves (RESP → SYS), reflected increased sympathetic modulation to the vessels. Our study in patients with IPF, compared with control subjects residing at the same altitude level, suggests that the autonomic alterations induced by the pathology persist or worsen despite the acute administration of oxygen. Based on our previous effort, indicating a relevant increase of TPR (Santiago-Fuentes et al., 2021), current research by the group is directed to analyze interactions including TPR. Finally, the proposed TFRi analysis may be used to better understanding the underlying physiological phenomena of different respiratory diseases.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by National Institute of Respiratory Diseases (INER). The patients/participants provided their written informed consent to participate in this study.