Association between autonomic control indexes and mortality in subjects admitted to intensive care unit

This study checks whether autonomic markers derived from spontaneous fluctuations of heart period (HP) and systolic arterial pressure (SAP) and from their interactions with spontaneous or mechanical respiration (R) are associated with mortality in patients admitted to intensive care unit (ICU). Three-hundred consecutive HP, SAP and R values were recorded during the first day in ICU in 123 patients. Population was divided into survivors (SURVs, n = 83) and non-survivors (NonSURVs, n = 40) according to the outcome. SURVs and NonSURVs were aged- and gender-matched. All subjects underwent modified head-up tilt (MHUT) by tilting the bed back rest segment to 60°. Autonomic control indexes were computed using time-domain, spectral, cross-spectral, complexity, symbolic and causality techniques via univariate, bivariate and conditional approaches. SAP indexes derived from time-domain, model-free complexity and symbolic approaches were associated with the endpoint, while none of HP variability markers was. The association was more powerful during MHUT. Linear cross-spectral and causality indexes were useless to separate SURVs from NonSURVs, while nonlinear bivariate symbolic markers were successful. When indexes were combined with clinical scores, only SAP variance provided complementary information. Cardiovascular control variability indexes, especially when derived after an autonomic challenge such as MHUT, can improve mortality risk stratification in ICU.

Various types of signal processing tools have been devised and applied to spontaneous fluctuations of cardiovascular variables recorded from short periods of time (i.e. about 5 minutes) in the attempt to characterize short-term cardiovascular control mechanisms 1,2 . Applications of these methods have been particularly successful in describing specific aspects of the cardiovascular control that are hardly quantifiable, such as vagal modulation 3 , or can be quantified with difficulty, such as cardiac baroreflex sensitivity 4 . The clinical utility of computing autonomic nervous system indexes in predicting the outcome of both cardiovascular and non-cardiovascular pathologies [5][6][7][8][9][10][11][12][13][14] and in stratifying the risk for adverse perioperative and postoperative cardiac events [15][16][17][18][19][20][21] has been proven by several studies. The clinical value of autonomic and cardiovascular control markers as risk predictors was confirmed in the emergency department 22,23 and critical care unit (ICU) [24][25][26][27][28] . However, it is unclear whether the application of an autonomic challenge in critical clinical settings might improve further their predictive and discriminative power.
The aim of this study is to test the association of autonomic control parameters derived from spontaneous fluctuations of cardiovascular variables recorded within the first day after admission in ICU with mortality in critically ill patients and their incremental value to traditional clinical scores in connection with an autonomic challenge. A large sample of techniques, commonly applied to the analysis of short beat-to-beat variability series complex on the ECG and locating with minimum jitters its peak using parabolic interpolation 39 , the temporal distance between two consecutive QRS complex apexes was computed and utilized as an approximation of the nth HP (HP n ). The maximum of arterial pressure within HP n was taken as the nth SAP (SAP n ). The amplitude of the first R-wave delimiting the HP n and measured from the isoelectric line provided the nth measure of R (R n ), as derived from the respiratory-related amplitude modulations of the ECG 40 . Fiducial points were carefully checked to avoid erroneous detections or missed beats. If isolated ectopic beats affected HP and SAP values, these measures were linearly interpolated using the closest values unaffected by ectopic beats. Sequences of about 300 consecutive HP, SAP and R values were randomly selected inside REST and MHUT sessions. From HP and SAP series we computed in time domain the HP and SAP means, μ HP and μ SAP . μ HP and μ SAP were expressed in ms and mmHg respectively. Then, HP, SAP and R series were linearly detrended before computing any additional cardiovascular variability index including HP and SAP variances termed σ 2 HP and σ 2 SAP . σ 2 HP and σ 2 SAP and expressed in ms 2 and mmHg 2 respectively. No additional preprocessing techniques were applied except for model-based complexity, model-based and model-free causality analyses requiring the division of the zero-mean HP, SAP and R series by their standard deviation to avoid any dependence of complexity and causality indexes on the amplitude of the spontaneous fluctuations.
Model-based spectral analysis. Power spectral analysis was carried out via a parametric method grounded on the identification of the coefficients of the autoregressive (AR) model and on a procedure decomposing the AR process into its basic components (see Supplement). The model order was optimized according to the Akaike figure of merit 41 in the range from 8 to 14. The power was computed in the low frequency (LF) band (i.e. 0.04-0.15 Hz) and in the high frequency (HF) band (i.e. 0.15-0.4 Hz) 1 . We computed the HF power of HP series (HF HP ) and the LF power of SAP series (LF SAP ) both expressed in absolute units (i.e. ms 2 and mmHg 2 respectively) as indexes of vagal and sympathetic modulation directed to the heart and vessels respectively 3,29 . We computed the HF power of the R series expressed in percent units (HF% R ) defined as the power of the R series in the HF band divided by the overall variance and the result multiplied by 100. Over the R series we estimated the respiratory rate (i.e. f R ) as the frequency of the dominant component in the HF band. f R was expressed in Hz.
Model-based and model-free complexity analyses. Complexity analysis was carried out through the computation of conditional entropy (CE) via a linear model-based and a nonlinear model-free approach (see Supplement). The model-based CE (MBCE) was again grounded on fitting HP and SAP series with an AR model 32 . The identification procedure and model order selection were carried out as for spectral analysis. As a quantity linked to the CE but computed via a nonlinear model-free method, we utilized the sample entropy (SampEn) 30 . At difference with MBCE, SampEn was able to account for the eventual presence of nonlinear dynamical features present in the series. SampEn was computed according to standard settings (i.e. embedding dimension d = 3, tolerance r = 0.2 × standard deviation of the series, lag between consecutive samples forming the pattern τ = 1, and Euclidean norm to calculate distances among patterns). Univariate symbolic analysis. Univariate symbolic analysis was carried out according to a method 31 using a uniform quantization approach over 6 bins, forming symbolic patterns by concatenating 3 consecutive symbols and grouping patterns into 4 families based on the number and sign of variations between adjacent symbols (see Supplement). The 4 pattern families are: i) no variation (0 V); ii) one variation (1 V); iii) two like variations (2LV); iv) two unlike variations (2UV). The percentage of 0 V, 1 V, 2LV and 2UV patterns (i.e. 0 V%, 1 V%, 2LV% and 2UV%) was obtained by dividing the number of patterns belonging to a class by the total number of patterns and, then, by multiplying the result by 100. Indexes were computed over HP and SAP series and labeled as 0 V% HP , 1 V% HP , 2LV% HP , 2UV% HP and 0 V% SAP , 1 V% SAP , 2LV% SAP , 2UV% SAP respectively. Model-based cross-spectral analysis. Cross-spectral analyses between HP and SAP and between HP and R were performed to estimate linear HP-SAP and HP-R dynamical relations. Cross-spectral analysis was carried out via a parametric method grounded on bivariate AR model (see Supplement). The model order was fixed 33 to 10. HP-SAP cross-spectral analysis allowed us to compute the transfer function from SAP to HP indicated as H HP-SAP (f) and the squared coherence function between HP and SAP indicated as K 2 HP-SAP (f). The averaged modulus of H HP-SAP (f) in the LF band, indicated as |H HP-SAP (LF)|, was taken as an estimate of baroreflex sensitivity in the LF band 34,42 . |H HP-SAP (LF)| was expressed in ms·mmHg −1 . The averaged K 2 HP-SAP (f) in the LF band, labeled as K 2 HP-SAP (LF) was taken as an estimate of the strength of the linear association between HP and SAP 43 . We computed also the averaged phase of H HP-SAP (f) in the LF band, denoted as PhH HP-SAP (LF). PhH HP-SAP (LF) was expressed in radians (rad). Negative phase values suggested that HP changes lagged behind SAP variations. The computation of |H HP-SAP (LF)| was performed without checking the prerequisites 34 of significant K 2 HP-SAP (LF) and negative PhH HP-SAP (LF). An analogous analysis was carried out to characterize the HP-R dynamical relation in the HF band. We computed |H HP-R (HF)|, K 2 HP-R (HF) and PhH HP-R (HF) taken respectively as markers of the gain, strength and phase of the cardiorespiratory relation 44 . Negative PhH HP-R (HF) values suggested that HP changes lagged behind R fluctuations.

Bivariate symbolic analysis.
To account for the eventual presence of nonlinear interactions between two series, the bivariate joint symbolic analysis 35 was applied over HP and SAP and over HP and R series. The bivariate joint symbolic analysis method exploited the univariate symbolic approach described in the Section "Univariate symbolic analysis" to build symbolic patterns over the two series (see Supplement). Then, joint symbolic schemes were formed by associating a pattern of one series with the pattern of the other starting τ cardiac beats ahead. Given the fastness of the baroreflex and cardiopulmonary interactions τ was set to 0 when both the dynamical interactions between HP and SAP and between HP and R were studied 45 . After defining as coordinated scheme as the one associating two patterns belonging to the same pattern family (i.e. 0V-0V, 1V-1V, 2LV-2LV and 2UV-2UV), their number was counted and their percentage over the class of coordinated patterns was tracked and indicated as 0V-0V%, 1V-1V%, 2LV-2LV% and 2UV-2UV%.

Model-based Granger causality (MBGC) analysis.
We exploited a traditional linear MBGC approach in the time domain (see Supplement) to assess the strength of the interactions from an input signal to an output one, while disambiguating the effect of a third one 46 . An AR model with two exogenous inputs was set in Ω = {HP,SAP,R} taking the HP as a target signal and SAP and R as exogenous ones. The variance of the HP prediction error in Ω was compared to the variance of the prediction error in Ω after excluding SAP through the log-causality ratio (logCR) (see Supplement). This index was labeled as logCR SAP → HP . The larger the logCR-SAP → HP , the greater the involvement of cardiac baroreflex 37 . An analogous index compared the variance of the HP prediction error in Ω with the variance of the prediction error in Ω after excluding R and was denoted as logCR R → HP . The larger the logCR R → HP , the greater the strength of the cardiopulmonary coupling 37 . The delays from SAP to HP and from R to HP were set to 0 beats to account for the fast vagal arm of baroreflex and cardiopulmonary pathway 45 . The model order was optimized in Ω in the range from 4 to 14 according to the Akaike figure of merit for multivariate processes 41 and, then, maintained in the restricted Ω obtained after excluding the presumed cause, while the coefficients of the model were identified again.

Model-free Granger causality (MFGC) analysis. MFGC approaches 47 extended the traditional MBGC
technique 46 by avoiding the a priori decision about the structure of the model and accounting for the eventual presence of nonlinear relations among signals in Ω. In Ω = {HP,SAP,R} we exploited the k-nearest neighbors approach 48 with a zero order local approximation 36 in a nonuniform embedding space 49,50 built incrementally according to a procedure maximizing the correlation between the original series and the predicted one 36 . The variance of the prediction error in the restricted Ω was computed after excluding from the optimal multivariate embedding space built in full Ω the components relevant to the presumed cause. The same indexes as those computed via MBGC method describing in the Sect. "Model-based Granger causality (MBGC) analysis" were computed (i.e. logCR SAP → HP and logCR R → HP ). The number of nearest neighbors k was fixed to 60. The image of the reference vector was excluded from the set of images of the k nearest neighbors utilized to predict the current values to limit overfitting. At difference from the MBGC approach the delays from SAP to HP and from R to HP were optimized according to the procedure exploited to construct the multivariate embedding space 36 . The maximum multivariate embedding dimension was fixed to 15 and the maximum number of lagged variables considered over the same signal was set to 10. Table 1 were tested according to χ2 test in the case of categorical variables and Mann-Whitney rank sum test in the case of continuous variables. Two-way repeated measures analysis of variance (one factor repetition, Holm-Sidak test for multiple comparisons) was used to check whether indexes exhibited between-group differences within the same experimental condition (i.e. REST or MHUT) and between-condition differences within the same group (i.e. SURV or NonSURV). Univariate Cox regression analysis was carried out by considering mortality as an outcome and the time elapsed from the admission to ICU to death for NonSURVs or to hospital discharge for SURVs, expressed in days, as independent time variable. Univariate Cox regression analysis was carried out over all variables, including the clinical ones that were found able to distinguish SURVs from NonSURVs. Those parameters found to be associated with mortality entered in a multivariate Cox regression model to assess their complementary value. Statistical analysis was carried out using a commercial statistical program (SPSS 22, IBM, Chicago, IL, USA). A p < 0.05 was always considered as significant.

Discussion
The main findings of the study can be summarized as follows: i) autonomic markers assessed during the first day after admission in ICU were associated with mortality; ii) even though the association was visible at REST, the exploitation of an orthostatic challenge, such as MHUT, amplified the possibility offered by autonomic markers; iii) SAP variability indexes were more likely to be associated with mortality, while HP variability indexes were useless; iv) nonlinear indexes computed via model-free approaches (i.e. symbolic analysis and SampEn) were more powerful than linear ones in stratifying the risk of patients in ICU; v) while markers derived from the HP-SAP variability interactions appeared to be helpful, those assessing HP-R variability relation were useless; vi) causality markers along cardiac baroreflex and cardiopulmonary pathway were not linked to the mortality and this conclusion held regardless of the method (i.e. MBGC or MFGC); vii) when all indexes showing a significant association with mortality were tested together in a multivariate Cox regression model, only SAPS II and SAP variance during MHUT carried complementary information.

Study hypotheses and rationale underlying autonomic challenge in ICU.
We hypothesized that cardiovascular and cardiorespiratory control indexes closely linked to the autonomic nervous system state could be valuable candidates to be included in a model predicting the risk of mortality in patients admitted to ICU and that an orthostatic maneuver operated at patient's bedside in the ICU could be fruitfully exploited to unveil the association of autonomic markers with mortality. If SURV and NonSURV groups could be separated according to one of the proposed autonomic control markers, that parameter should be a valuable candidate in a multivariate clinical model predicting the risk of mortality for people admitted to ICU and its incremental value to clinical scores traditionally utilized in ICU, such as the SOFA score and SAPS II, should be tested. The separation between SURVs and NonSURVs cannot be trivially assumed given that individuals admitted to ICU might exhibit an impaired autonomic control in relation to their pathological state and might have undergone a pharmacological treatment, such as the administration of catecholamines, and/or an intervention, such as sedation, interfering with the ability of autonomic function to govern physiological variables. Given the presumed low level of residual ability of the autonomic nervous system to regulate physiological variables in critical patients, we supposed that the application of an autonomic challenge such as the MHUT 38 could be helpful to amplify the difference between those subjects who preserved autonomic control and those who did not. However, if the separation could be obtained just at REST, the candidate variable should be considered more powerful than the one leading to a separation solely during MHUT given that no stressor is needed, thus simplifying the test procedure.
Distinguishing SURVs from NonSURVs during their first day in ICU. We found that autonomic nervous system markers, although depressed in both SURVs and NonSURVs, can be fruitfully exploited to distinguish the two groups and MHUT can favor this separation. The list of the parameters allowing this distinction comprises σ 2 SAP , SampEn SAP , 0 V% SAP , 1 V% SAP , 0V-0V% HP-SAP solely during MHUT, 2LV% SAP and 2LV-2LV% HP-SAP both at REST and during MHUT. Since more variables distinguished SURVs from NonSURVs during MHUT than at REST, we conclude that an orthostatic maneuver carried out at the patient's bedside is helpful toward the stratification of the mortality risk in people admitted to ICU. We also observed that SAP variability is more helpful than the HP one. Indeed, none of the variability markers involving HP series could distinguish NonSURVs from SURVs. This result stresses the relevance of analyzing SAP series in addition to the HP one 51 . The variance of SAP during MHUT was significantly larger in NonSURVs than in SURVs, thus indicating a greater difficulty of NonSURVs in keeping under control SAP fluctuations and a reduced capacity of NonSURVs in coping with external and/or internal disturbances. Also symbolic analysis of SAP variability evidenced that SAP values are less stable in NonSURVs than in SURVs: as a matter of fact, the likelihood of more stable SAP patterns decreased and the rate of more variable SAP features increased in NonSURVs. Complexity indexes describing the HP-SAP variability interactions via joint symbolic analysis showed a remarkable power in distinguishing SURVs from NonSURVs. Indeed, during MHUT the HP-SAP coupling at slow time scales, as monitored via 0V-0V% HP-SAP , was significantly smaller in NonSURVs than in SURVs, while the opposite situation was observed at faster time scales as indicated by the trend of 2LV-2LV% HP-SAP . This finding once again emphasized the pathological response of NonSURVs to the orthostatic stimulus that usually leads in healthy population to an increase of 0V-0V% HP-SAP and a decrease of 2LV-2LV% HP-SAP as an indication of the activation of the cardiac baroreflex 52 .
Remarkably, nonlinear model-free markers appear to be more powerful than linear model-based ones in distinguishing SURVs from NonSURVs. For example, complexity markers computed via SampEn were more powerful than those calculated via MBCE approach and univariate symbolic indexes were more effective than spectral ones. A model-free approach based on joint symbolic analysis was able to separate the two groups, while linear model-based methods assessing the gain and phase of the transfer function and Granger causality indexes could   Table 2. Results of univariate Cox regression analysis. SOFA = sepsi-related organ failure assessment; SAPS = simplified acute physiology score; REST = at rest in supine position; MHUT = modified head-up tilt; HP = heart period; SAP = systolic arterial pressure; σ 2 SAP = variance of SAP; SampEn SAP = sample entropy of SAP; 0 V% SAP , 1 V% SAP and 2LV% SAP = percentage of 0 V, 1 V, and 2LV patterns respectively as detected by univariate symbolic analysis of SAP; 0V-0V% HP-SAP and 2LV-2LV% HP-SAP = percentage of 0V-0V and 2LV-2LV patterns respectively as detected by bivariate symbolic analysis of HP and SAP. not. Moreover, the list of parameters helpful in distinguishing SURVs from NonSURVs stresses that the more elaborated the method, the weaker its ability in separating SURVs from NonSURVs. Indeed, markers of cardiac baroreflex, assessing the HP-SAP gain, HP-SAP phase and strength of global association between HP and SAP, and indexes of causality along cardiac baroreflex, estimating the strength of the linear association in the temporal direction from SAP to HP, could not distinguish SURVs from NonSURVs and, at the best, could detect the effect of MHUT. The same consideration holds for the analysis of HP-R variability interactions. This observation is in line with a recent contribution 14 suggesting that in clinical applications and in pathological groups the robustness of the index could be more important than its methodological design. At this regard joint symbolic analysis might provide a good tradeoff between the complexity of the approach allowing the description of dynamical interactions among physiological variables and the robustness of the method.
Association of cardiovascular control indexes with mortality and their complementary value to traditional risk scores. All cardiovascular control markers able to differentiate SURVs from NonSURVs were found significantly associated with mortality via univariate Cox regression analysis, thus suggesting that they can be potentially helpful in stratifying the mortality risk of critical patients admitted to ICU. However, only σ 2 SAP during MHUT carried complementary information to traditional risk indicators that were able to distinguish the two groups as well (i.e. the presence of septic shock, SOFA score and SAPS II). This result stressed the relevance of monitoring short-term SAP variability and of applying an autonomic challenge in ICU. The hypovolemic condition evoked by MHUT 38 is likely to have an impact much stronger in NonSURVs who could not appropriately deal with it, as suggested by the larger amplitude of SAP oscillations that were not efficiently buffered due to a more impaired cardiovascular control. We recommend the application of this autonomic challenge to gain additional information and the inclusion of σ 2 SAP measured during the challenge in mortality risk scores.

Effects of MHUT in patients during their first day in ICU.
MHUT induces an increase of SAP fluctuations in the LF band and a decrease of cardiac baroreflex sensitivity in healthy subjects with age similar to those of the present study, being these findings compatible with a sympathetic activation directed to the vessel and with a cardiac baroreflex unloading associated with the challenge 38 . In healthy subjects the central hypovolemia evoked by an orthostatic challenge led to an increased strength of the HP-SAP variability relation along cardiac baroreflex 37,52,53 , likely related to the augmented involvement of this reflex in preserving adequate arterial pressure levels 54,55 , a decreased strength of the cardiorespiratory coupling 45,53 and a diminished complexity of cardiac control 32,56 , likely related to the vagal withdrawal associated with the challenge 54,57-59 . Conversely, the complexity of vascular control in healthy individuals was not affected by the orthostatic challenge 56,60 likely because arterial pressure regulatory mechanisms are not under vagal control. Remarkably, in SURVs a subset of these modifications in response to MHUT could be observed, thus suggesting that autonomic control preserved a certain ability to govern the changes of physiological variables in this subgroup of critically ill individuals. For example, cardiac baroreflex sensitivity and complexity of HP series decreased significantly in SURVs during MHUT. Conversely, the decrease of the amount of linear association between HP and SAP in a frequency band typical of the baroreflex functioning (i.e. LF band) during MHUT, observed in both SURVs and NonSURVs, suggests an impairment of the cardiovascular control in both groups and a greater isolation of physiological systems typical of an organism that lost its ability to provide an integrated response to a stressor. An additional sign of the autonomic control impairment in both SURVs and NonSURVs is the inability of observing the expected modifications of the strength of causal relation along cardiac baroreflex and cardiopulmonary pathway during MHUT 37,45,53 . MHUT highlighted some peculiar differences between SURVs and NonSURVs. For example, the larger increase of the magnitude of the SAP variability during MHUT in NonSURVs should be interpreted as a greater inability of the regulatory systems in this group in limiting arterial pressure instabilities compared to SURVs more than a sign of the sympathetic activation evoked by the stimulus. In addition, the decrease of complexity of the vascular control observed only in NonSURVs indicated an excessive simplification of the mechanisms of arterial pressure regulation, thus stressed again the limited regulatory resources of this group.

Conclusions
The stratification of mortality risk at the admission in ICU is commonly carried out by several scores that do not account for autonomic control markers. The rationale underlying this choice is that autonomic control mechanisms might be impaired and/or depressed in the category of individuals admitted to ICU as a consequence of their pathological condition and/or pharmacological intervention to such a level that the information carried by autonomic control indexes might be negligible or meaningless. Conversely, this study suggests that markers of autonomic regulation computed over SAP variability are associated with mortality in patients admitted to ICU and this association is even stronger whether an orthostatic challenge is applied at the bedside to stimulate an autonomic reaction. More specifically, we recommend the application of an autonomic challenge at the bedside in ICU and the inclusion of the amplitude of SAP variability into more traditional scores of mortality risk given its complementary value.
Data availability. Data are fully available through the corresponding author.