Evaluation of the impact of surgical aortic valve replacement on short-term cardiovascular and cerebrovascular controls through spontaneous variability analysis

We assessed the effect of surgical aortic valve replacement (SAVR) on cardiovascular and cerebrovascular controls via spontaneous variability analyses of heart period, approximated as the temporal distance between two consecutive R-wave peaks on the electrocardiogram (RR), systolic, diastolic and mean arterial pressure (SAP, DAP and MAP) and mean cerebral blood flow (MCBF). Powers in specific frequency bands, complexity, presence of nonlinear dynamics and markers of cardiac baroreflex and cerebral autoregulation were calculated. Variability series were acquired before (PRE) and after (POST) SAVR in 11 patients (age: 76±5 yrs, 7 males) at supine resting and during active standing. Parametric spectral analysis was performed based on the autoregressive model. Complexity was assessed via a local nonlinear prediction approach exploiting the k-nearest-neighbor strategy. The presence of nonlinear dynamics was checked by comparing the complexity marker computed over the original series with the distribution of the same index assessed over a set of surrogates preserving distribution and power spectral density of the original series. Cardiac baroreflex and cerebral autoregulation were estimated by assessing the transfer function from SAP to RR and from MAP to MCBF and squared coherence function via the bivariate autoregressive approach. We found that: i) orthostatic challenge had no effect on cardiovascular and cerebrovascular control markers in PRE; ii) RR variance was significantly reduced in POST; iii) complexity of SAP, DAP and MAP variabilities increased in POST with a greater likelihood of observing nonlinear dynamics over SAP compared to PRE at supine resting; iv) the amplitude of MCBF variations and MCBF complexity in POST remained similar to PRE; v) cardiac baroreflex sensitivity decreased in POST, while cerebrovascular autoregulation was preserved. SAVR induces important changes of cardiac and vascular autonomic controls and baroreflex regulation in patients exhibiting poor reactivity of cardiovascular regulatory mechanisms, while cerebrovascular autoregulation seems to be less affected.


Introduction
Since sympathetic activation and vagal withdrawal lead to the decrease of the variability of heart period, approximated as the temporal distance between two consecutive R-wave peaks on the electrocardiogram (RR), and to the increase of the fluctuations of systolic arterial pressure (SAP), the RR and SAP variances have been exploited to infer vagal and sympathetic neural controls [1][2][3][4][5][6]. Spectral analysis allowed a more precise link of RR and SAP variabilities with the state of the autonomic nervous system given that oscillations in the high frequency (HF, from 0.15 to 0.4 Hz) band of the RR series are more specifically associated to vagal control directed to the sinus node [1][2][3][4][5][6], while those in the low frequency (LF, from 0.04 to 0.15 Hz) band of the SAP series are more specifically linked to sympathetic control directed to the vessels [4][5][6]. Complexity and nonlinear content of the RR and SAP variability series are frequently explored to gain insight into peculiar features of autonomic regulatory mechanisms that cannot be fully described by power spectral density [7][8][9][10][11][12][13]. Complexity of the cardiac control is reduced during sympathetic activation and vagal withdrawal [7,8]. The development of nonlinear dynamics appears to be favored by an enhancement of the vagal control [9][10][11][12] and by a weakened action of the cardiac baroreflex [13]. Also the different response of regulatory mechanisms to positive or negative SAP variations and circuits imposing a certain degree of cardiorespiratory phase coupling might play an important role in producing nonlinear RR variability patterns [14][15][16].
The integrative regulation of cerebral blood flow (CBF) comprising chemoreflex, neuronal metabolism, cerebral autoregulation and autonomic control [17] is routinely assessed via noninterventional techniques based on the analysis of the spontaneous variability of physiological variables related to cerebral circulatory system such as mean arterial pressure (MAP), acquired via volume-clamp photoplethysmography from the middle finger [4], and mean CBF (MCBF) [18], estimated via transcranial Doppler device from middle cerebral arteries [19,20]. The assessment of the MAP and MCBF powers and of the MCBF-MAP dynamical relation in specific frequency bands is the basis for the evaluation of the ability of cerebral vasculature to buffer MCBF via suitable modifications of cerebral resistance [18]. Given the important contribution of the autonomic nervous system to integrative regulation of CBF [17,21], the concurrent evaluation of cerebrovascular control with autonomic markers derived from RR and SAP variability series might provide more insight into cerebrovascular control mechanisms and complex interactions among cardiovascular and cerebrovascular regulations.
A deep analysis of the impact of cardiac surgery on cardiovascular and cerebrovascular controls is missing. Information is mainly limited to cardiovascular control. Indeed, spectral analysis of RR variability suggested that cardiac surgery depresses vagal control [22][23][24][25] as well as baroreflex function [26], thus exposing patients to a higher risk of postoperative adverse events such as atrial fibrillation and acute kidney dysfunction [26,27]. Complexity analysis of RR and SAP variability series was less frequently performed after cardiac surgery. Complexity analysis of RR and SAP variability series suggested a more variable response to cardiac surgery across individuals and groups [26,28]. Information about the impact of cardiac surgery on cerebrovascular control is even more limited [29]. Given the postoperative impairment of the autonomic control [22][23][24][25][26] and the relevant impact of the autonomic function on integrative regulation of CBF [17,21], we hypothesize a dramatic influence of cardiac surgery on the magnitude of MAP and MCBF variations and/or MCBF-MAP dynamical relation. The assessment of cerebrovascular control after cardiac surgery might provide more information about the patient's post-operative state and, if this characterization was carried out in association with that of the cardiovascular one, the description might be more complete and insightful.
The aim of this study is to characterize cardiovascular and cerebrovascular controls via spectral and complexity analyses and to test the presence of nonlinear patterns from variability series recorded before and after surgical aortic valve replacement (SAVR). Variability of RR, SAP, and diastolic arterial pressure (DAP) were analyzed to infer the state of cardiac and vascular controls, while the variability of MCBF and MAP was assessed to infer that of cerebrovascular regulation. Nonlinear dynamics were detected through a local nonlinear prediction marker [10,30] in association with a surrogate data approach [10,31]. The analysis was completed with the description of the baroreflex control and cerebral autoregulation via the transfer function method applied to SAP and RR variability series [32] and MAP and MCBF variability series [18] respectively. Preliminary results were presented at the 11 th meeting of the European Study Group of the Cardiovascular Oscillations [33] and at the 42 nd Annual International Conference of the Engineering in Medicine and Biology Society [34].

Ethics statement
The study was in keeping with the Declaration of Helsinki. The study was approved by the ethical review board of the San Raffaele Hospital, Milan, Italy (approval number: 68/int/2018; approval date: 05/04/2018) and authorized by the Policlinico San Donato, San Donato Milanese, Milan, Italy (authorization date: 13/04/2018). Written signed informed consent was obtained from all subjects.

Population and experimental protocol
We enrolled 11 patients (age: 76±5 yrs, 7 males) undergoing SAVR at the IRCCS Policlinico San Donato, San Donato Milanese, Milan, Italy. Demographic and clinical data of the SAVR group were reported in Table 1. They did not feature either atrial fibrillation, overt autonomic nervous system pathologies or cerebrovascular diseases. We acquired electrocardiogram (ECG) from lead II (BioAmp FE132, ADInstruments, Australia), non-invasive finger arterial pressure (AP) by volume-clamp photoplethysmography (CNAP Monitor 500, CNSystems, Austria), respiration (RESP) via a thoracic piezoelectric belt (ADinstruments, Australia) and CBF velocity via a transcranial Doppler device (Multi-Dop X, DWL, San Juan Capistrano, CA, USA) from the left or right middle cerebral artery. Signals were sampled at 400 Hz through a commercial acquisition system (Power Lab, ADInstruments, Australia). Signals were recorded 1 day before SAVR (PRE) and within 7 days after SAVR (POST). Acquisition sessions comprised recordings at rest in supine position (REST) and during active standing (STAND). REST and STAND lasted 10 minutes and REST always preceded STAND. In PRE REST and STAND sessions were carried out in all subjects. In POST REST was performed in 8 individuals and STAND in 6 subjects due to the physical and psychological debilitation of some patients. ECG and AP were recorded in all the subjects present in a given experimental session. Conversely, CBF was recorded in PRE in 10 subjects out of 11 at REST and in 8 individuals out of 11 during STAND and in POST in 6 subjects out of 8 at REST and in 4 individuals out of 6 during STAND. These figures were due to the difficulty in locating either left or right middle cerebral artery.

Extraction of beat-to-beat variability series and time domain indexes
After detecting the R-wave peaks from the ECG with a classical method based on a threshold applied to the first derivative, the ith RR, where i is the progressive cardiac beat number, was derived as the time distance between the (i-1)th and the ith cardiac beat. The ith SAP and DAP were measured, respectively, as the AP maximum within the ith RR and the AP minimum following the ith SAP. The ith MAP was obtained as the integral of AP between the (i-1)th and ith DAP fiducial points and by dividing the result by the interdiastolic time interval. The ith RESP was obtained by sampling RESP signal at the apex of the first R-wave peak delimiting the ith RR. The ith MCBF was obtained as the integral of CBF between the (i-1)th and ith minima detected over CBF and closest in time to (i-1)th and ith DAP fiducial points and by dividing the result by the time distance between the two minima [35]. The RR, SAP, DAP, MAP and MCBF series were manually checked and corrected in case of missing beats or misdetections. Effects of ectopic beats or isolated arrhythmic events were mitigated via linear interpolation. Spectral, cross-spectral, complexity analyses as well as detection of nonlinear dynamics were carried out directly over synchronous sequences lasting 256 consecutive beats randomly selected within the whole recordings. In time domain we computed the means, indicated as μ RR , μ SAP , μ DAP , μ MAP and μ MCBF , and the variances, denoted with σ 2 RR , σ 2 SAP , σ 2 DAP , σ 2 MAP and σ 2 MCBF . Means and variances were expressed, respectively, in ms, mmHg, mmHg, mmHg and cm�s -1 and ms 2 , mmHg 2 , mmHg 2 , mmHg 2 and cm 2 �s -2 . After computing the mean, the series were linearly detrended before computing variance.

Univariate model-based frequency domain analysis
Univariate model-based frequency domain analysis was carried out via a traditional parametric power spectral method (see S1 Appendix). Briefly, variability series were described as a realization of an autoregressive (AR) process modeling the variation of the most recent value of the series about its mean as a linear combination of p past changes weighted by constant coefficients plus a sample drawn from a realization of a zero mean white noise, where p is the order of AR model [36,37]. The coefficients of the AR model and the variance of the white noise were identified directly from the series by solving the least squares problem via Levinson-Durbin recursion [36]. The number of coefficients p was chosen according to the Akaike's figure of merit in the range from 8 to 16 [38]. Power spectral density was computed from the AR coefficients and from the variance of the white noise according to the maximum entropy spectral estimation approach [36]. The power spectral density was factorized into a sum of terms, referred to as spectral components, the sum of which provides the entire power spectral density [37]. Power spectral decomposition provided the central frequency of the components expressed in normalized frequency units, namely cycles per beat. Central frequency ranged from 0 to 0.5 cycles per beat and was converted into Hz by dividing the value by the average sampling period T of the variability series, i.e. T = μ RR , expressed in s [37]. A spectral component was attributed to a specific frequency band if its central frequency lay in that band. If multiple spectral components were found to belong to the same frequency band, their powers were summed up and the weighted central frequency was computed.

Bivariate model-based frequency domain analysis
Bivariate model-based frequency domain analysis was carried out via a traditional parametric cross-spectral method (see S1 Appendix). This method allows the assessment of the input-output relation in the frequency domain (i.e. the transfer function) between two series when one is assumed to be the cause and the other is taken as the effect [32]. The dynamics of input and output series about their mean values were jointly described as a bivariate AR process [39] modeling the variation of the most recent value of one series about its mean as a linear combination of p past changes of the same series and p past variations of the other series weighted by constant coefficients plus a sample drawn from a realization of a zero mean white noise, where p is the order of the bivariate AR model. Although the bivariate AR model described the closed loop relation via the representation of the feedforward and feedback arms [37], we focused our attention on a specific direction of interactions by assigning a priori one series as the input and the other as the output [32]. The cross-spectrum from the input to the output and the autospectra of both input and output were computed from the coefficients of the bivariate AR model and from the variance of the white noises [39]. The model order p was fixed at 10 and the coefficients of the bivariate AR model were identified via least squares approach solved using Cholesky decomposition method [37,39]. The transfer function (TF) was estimated as the ratio of the cross-spectrum computed from the input to the output to the power spectrum of input series [40]. The TF modulus (TFM) and phase were calculated as a function of the frequency. The TFM could not be negative and was expressed in units being the ratio of the unit of the output to that of the input. The phase was expressed in radians (rad) and ranged between -π and +π, with negative value indicating that the output lagged behind the input. The ratio of the squared cross-spectrum modulus to the product of the power spectra of the input and output series provided the estimation of the squared coherence function (K 2 ) as a function of the frequency [40]. The K 2 was dimensionless and ranged between 0 and 1, where 0 indicated full uncoupling and 1 perfect association between the input and the output. The TFM, phase and K 2 functions were sampled in correspondence of the frequency where the K 2 peaked the maximum value within the considered frequency band [41].

Cardiovascular variability and cardiac baroreflex markers
Cardiovascular regulation was assessed in the low frequency (LF, from 0.04 to 0.15 Hz) and high frequency (HF, from 0.15 to 0.4 Hz) bands [2]. The power of the RR series in the HF band expressed in absolute units (HFa RR ) was taken as a marker of the vagal modulation directed to the heart [1,4,6] and the power of the SAP and DAP series in the LF band expressed in absolute units (LFa SAP and LFa DAP ) was taken as an index of sympathetic modulation directed to the vessels [5,6]. The central frequency of the dominant component of the RESP series in the HF band was taken as an estimate of the respiratory rate (f RESP ). HFa RR , LFa SAP , LFa DAP and f RESP were expressed in ms 2 , mmHg 2 , mmHg 2 and Hz respectively. The baroreflex function was assessed through the computation of the TFM from SAP to RR in the LF and HF bands [TFM RR-SAP (LF) and TFM RR-SAP (HF)], the TF phase from SAP to RR in the LF and HF bands [Ph RR-SAP (LF) and Ph RR-SAP (HF)] and K 2 between SAP and RR in the LF and HF bands [K 2 RR-SAP (LF) and K 2 RR-SAP (HF)]. TFM RR-SAP (LF) and TFM RR-SAP (HF) were taken as an estimate of the baroreflex sensitivity [32,40,42] and expressed in ms�mmHg -1 . While Ph RR-SAP (LF) and Ph RR-SAP (HF) were considered markers of the delay/advance of SAP on RR [43], K 2 RR-SAP (LF) and K 2 RR-SAP (HF) were measures of the RR-SAP coupling strength [44].

Complexity analysis based on local nonlinear prediction
From the series y = {y i , i = 1,. . .,N}, where N is the series length, we built the pattern set y = {y i = (y i , y i-1 , ���, y i-L+1 ), i = L,. . .,N}, where L is the pattern length and y i is an ordered vector of L consecutive samples taken from y. Given the reference pattern y i , its k nearest neighbors {y k } were searched for and the images of the reference pattern and its k nearest neighbors, namely y i+1 and {y k+1 }, were retained. We defined the best prediction of image y i+1 of the reference pattern as the weighted average of the images {y k+1 } of its k nearest neighbors, where the weights were the distances of each k nearest neighbor from the reference pattern [30]. The ability to predict was quantified as the complement to 1 of the squared correlation coefficient between the original series and the predicted one [10]. According to [10,48] k was set to 30, the distance was computed using the Euclidian norm and solely the reference vector was excluded from the set of its nearest neighbors. This cost function exhibited a minimum resulting from two contrasting tendencies, namely the ability to predict future values increasing with L and the growing dispersion of the patterns in the L-dimensional space resulting from the wider and wider volume occupied by the points with L [10]. The minimum of this cost function over L was taken as normalized complexity index (NCI) [10]. The greater the NCI and the closer to 1, the more unpredictable and complex the series. The smaller the NCI and the closer to 0, the more regular and predictable the series. NCI was computed over all the series and the indexes were labelled NCI RR , NCI SAP , NCI DAP , NCI RESP , NCI MAP and NCI MCBF .

Testing the null hypothesis of linear dynamic
We tested the null hypothesis of linear dynamic using a surrogate data approach in connection with a nonlinear index such as NCI [10,48]. One hundred surrogate series were generated from the original series according to the amplitude-adjusted iteratively-refined Fourier transform method [31]. This method iterates the following steps [31]: i) the Fourier transform is applied to the series; ii) the Fourier phases are substituted with numbers drawn from a uniform distribution bounded between 0 and 2π; iii) the inverse Fourier transform is applied to return in the time domain, thus generating the surrogate; iv) the distribution of the surrogate is constrained to be exactly equal to that of the original series by substituting any sample of the surrogate with that of the original series that occupies the same position according to a rank ordering procedure applied to both surrogate and original series; v) all the previous steps are repeated until the power spectrum of the surrogate and that of the original sequence match or differences between them are below of a given threshold. It was demonstrated in [31] that the procedure converges to produce a surrogate with the same power spectrum and distribution as the original series and a very good approximation of the original power spectrum with perfect preservation of the original distribution could be achieved after 100 iterations. The NCI computed over the original series was compared with the distribution of the NCI values calculated from the surrogate series. If the NCI computed over the original series was below the 5 th percentile of the NCI distribution calculated over the surrogates, the null hypothesis of linear dynamic was rejected and the alternative hypothesis of nonlinear dynamic was accepted. The rationale is that the deterministic nonlinear features, destroyed in the surrogates by the phase randomization procedure, were predicted better over the original series by the local nonlinear prediction method than over the surrogates. The percentage of nonlinear dynamics (NL%) was monitored over all the series and the indexes were labelled NL% RR , NL% SAP , NL% DAP , NL % RESP , NL% MAP , and NL% MCBF .

Statistical analysis
Two-way analysis of variance (Holm-Sidak test for multiple comparisons) was applied to variability indexes to detect the effect of SAVR within the same experimental condition (i.e. REST or STAND) and the response to postural challenge within the same period of analysis (i.e. PRE or POST). A p<0.05 was always considered statistically significant. χ2 test (Fisher exact test) was applied over the categorical variables (i.e. presence/absence of nonlinear dynamics) to assess the effect of SAVR regardless of the posture and the influence of the postural challenge regardless of the period of analysis. If appropriate, the level of significance (i.e. 0.05) of the χ2 test was reduced according to the number of comparisons (i.e. 4) to account for the multiple comparison issue. Statistical analysis was carried out using a commercial statistical program (Sigmaplot, v.14.0, Systat Software, Inc., Chicago, IL, USA). Table 2 reports the results of time and frequency domain analyses carried out over RR, SAP, DAP and RESP series at REST and during STAND in PRE and POST sessions. Regardless of the experimental condition (i.e. REST or STAND), μ RR significantly decreased in POST compared to PRE. σ 2 RR decreased and μ DAP increased during POST with respect to PRE but the rise was significant solely at REST and during STAND respectively. No additional time domain markers were able to detect POST-PRE changes. The sole time domain index that was able to highlight a significant effect of the orthostatic challenge was μ DAP that increased during STAND compared to REST in the POST session. None of the frequency domain markers, including f RESP was able to separate either experimental conditions (i.e. REST and STAND) or periods of analysis (i.e. PRE and POST). Table 3 reports the results of time and frequency domain analyses carried out over MAP and MCBF series at REST and during STAND in PRE and POST sessions. None of time domain indexes was able to detect either the effect of the postural challenge or the impact of the cardiac surgery. The sole frequency domain markers modified by the cardiac surgery were HFa MAP and HF% MAP : HFa MAP increased in POST compared to PRE solely during STAND, while HF% MAP raised in POST compared to PRE both at REST and during STAND.

Results
The grouped error bar graphs of Fig 1 show the NCI computed over RR (Fig 1A), SAP ( Fig  1C), DAP ( Fig 1E) and RESP (Fig 1G) series, while the grouped bar graphs of Fig 1 show the NL % over RR (Fig 1B), SAP (Fig 1D), DAP ( Fig 1F) and RESP (Fig 1H) series during PRE and POST sessions. Markers were computed at REST (solid black bars) and during STAND (solid white bars). NCI RR and NCI RESP did not vary in response to cardiac surgery. NCI SAP rose significantly after SAVR both at REST and during STAND and the POST-PRE increase of NCI DAP was significant only during STAND. Within the same period of analysis (i.e. PRE or POST) posture modification did not affect any NCI marker. NL% SAP was significantly higher at REST in the POST session compared to PRE and lower in POST session during STAND compared to REST. NL% RR , NL% DAP and NL% RESP did not vary with postural challenge and SAVR surgery. Fig 2 has the same structure of Fig 1 but it shows the NCI computed over MAP (Fig 2A) and MCBF (Fig 2C) series and NL% over MAP ( Fig 2B) and MCBF (Fig 2D) series during PRE and POST sessions. NCI MCBF did not vary in response to cardiac surgery. NCI MAP rose significantly after SAVR both at REST and during STAND. Within the same period of analysis (i.e. PRE or POST) posture did not affect any NCI marker. NL% MAP and NL% MCBF did not vary with postural challenge and SAVR surgery.
When all series were pooled together regardless of the period of analysis, NL% was significantly higher at REST than during STAND. Conversely, when all series were pooled together regardless of the experimental condition, no effect of the surgery over NL% was detectable.  The grouped error bar graphs of Fig 3 show K 2 RR-SAP (LF) (Fig 3A), K 2 RR-SAP (HF) (Fig 3B), Ph RR-SAP (LF) (Fig 3C), Ph RR-SAP (HF) (Fig 3D), TFM RR-SAP (LF) (Fig 3E) and TFM RR-SAP (HF) (Fig 3F) during PRE and POST sessions. Markers were computed at REST (solid black bars) and during STAND (solid white bars). Orthostatic stimulus and cardiac surgery did not influence either K 2 or phase markers. Regardless of the frequency band (i.e. LF or HF) SAVR surgery decreased TFM both at REST and during STAND. Assigned the period of analysis (i.e. PRE or POST) the impact of postural maneuver on the baroreflex sensitivity as measured via TFM was negligible in both LF and HF bands. Fig 4 has the same structure as Fig 3 but it shows K 2 MCBF-MAP (VLF) (Fig 4A), K 2 MCBF-MAP (LF) (Fig 4B), K 2 MCBF-MAP (HF) (Fig 4C), Ph MCBF-MAP (VLF) (Fig 4D), Ph MCBF-MAP (LF) (Fig 4E), Ph MCBF-MAP (HF) (Fig 4F), TFM MCBF-MAP (VLF) (Fig 4G), NL% RESP (h) in PRE and POST. The solid black and white bars are relevant to markers computed at REST and during STAND respectively. Data in the grouped error bar graphs are given as mean plus standard deviation. The symbol � indicates a significant modification between PRE and POST within the same experimental condition with p<0.05, while the symbol § indicates the significant REST-STAND variation within the same period of analysis with p<0.05. https://doi.org/10.1371/journal.pone.0243869.g001

PLOS ONE
TFM MCBF-MAP (LF) (Fig 4H) and TFM MCBF-MAP (HF) (Fig 4I) during PRE and POST sessions. Regardless of the frequency band (i.e. VLF, LF or HF) K 2 , phase and TFM of the MCBF-MAP dynamical relation did not vary with posture and SAVR surgery.

Discussion
The major findings can be summarized as follows: i) orthostatic maneuver did not modify cardiovascular regulatory indexes in PRE, thus suggesting that SAVR patients featured a depressed vagal autonomic and baroreflex controls just before the cardiac surgery; ii) vagal autonomic and cardiac baroreflex impairment were observed after SAVR; iii) cerebrovascular variability and cerebral autoregulation appeared to be less affected by SAVR; iv) at the vascular control level the autonomic control impairment took the form of a post-surgery increased complexity of SAP, DAP and MAP variabilities and an augmented presence of nonlinear dynamics of SAP after SAVR; v) STAND reduced the likelihood of finding nonlinear dynamics; vi) respiratory rate was not affected by either SAVR surgery or orthostatic challenge.

Effect of STAND on a cohort of patients undergoing SAVR surgery
Orthostatic challenges did not provoke the expected decrease of vagal modulation directed to the heart [1,3,4,6,49], the expected increase of sympathetic modulation directed to the vessels [4,6] and the expected diminution of baroreflex sensitivity [4,43,50,51]. Indeed, STAND induced negligible effects over cardiovascular markers in PRE. This result can be taken as a hallmark of cardiovascular control dysfunction in our cohort of patients just before SAVR surgery. The effect of STAND on cerebrovascular variability and cerebral autoregulation is analogously limited. This result might suggest a possible impairment of the cerebrovascular control mechanisms given that a significant increase of the MCBF-MAP coupling strength was found when this marker of the dynamical MCBF-MAP relation was assessed in a healthy group [35].

Effect of SAVR on the amplitude of cardiovascular and cerebrovascular variability markers
It is well-known that SAVR reduced cardiac vagal control [22][23][24][25][26]. This observation is confirmed in this study by the dramatic postoperative reduction of σ 2 RR observed at REST. This observation is important because our subjects exhibited an impaired cardiac vagal control already in PRE as indicated by the negligible STAND-REST variation of the HFa RR power [1,3,4,6,49]. The post-surgery depression of vagal regulation might expose our group to possible postoperative arrhythmic events [26,27]. The unaltered value of σ 2 SAP , σ 2 DAP and σ 2 MAP during POST regardless of the experimental condition points toward a more preserved sympathetic control directed to the vessels. These findings are compatible with the postoperative pharmacological therapy as well as the physical debilitation and the high perceived level of stress associated to the POST condition. The amplitude of the MCBF fluctuations seems to be preserved as indicated by the invariable amount of σ 2 MCBF and the unchanged MCBF powers in the various frequency bands in POST compared to PRE. This finding could be related to a better preservation of the sympathetic control after SAVR. Indeed, MCBF regulation is under sympathetic control [21], even though different mechanisms might contribute to preserve cerebral autoregulation including the postoperative improvement of the valve hemodynamics and aortic blood flow profile [17].

Effect of SAVR on the complexity and nonlinearities of cardiovascular and cerebrovascular variabilities
In spite of the maintenance of the magnitude of the SAP, DAP and MAP oscillations, important modifications of the SAP, DAP and MAP dynamics were observed. Indeed, complexity of SAP, DAP and MAP variability increased after SAVR and this increase was accompanied by a significant rise of SAP nonlinear dynamics, especially visible at REST. Even a tendency toward an increased likelihood of observing nonlinearities in MCBF series was visible in POST. We conclude that even sympathetic control seems to be affected by SAVR. An increased SAP variability complexity could be attributed to the inability of the sympathetic control to produce synchronous oscillations of peripheral resistances by producing rhythmical alternations between suitable vasoconstriction and vasodilatation episodes [52] as suggested in some pathological situations [7]. It is worth noting that, in spite of the limited response to STAND, we observed the typical decrease of the likelihood of detecting nonlinear dynamics in cardiovascular variability reported in healthy individuals during orthostatic challenge [9,10], thus indicating a residual cardiovascular control in both PRE and POST.

Effect of SAVR on the baroreflex control and cerebral autoregulation markers
Markers of cardiac baroreflex and cerebral autoregulation were evaluated via bivariate frequency domain metrics assessing phase and TFM from an input to an output and their degree of association via K 2 [32]. In the case of the assessment of the cardiac baroreflex the input is SAP variability and the output is the RR variability [40,[42][43][44], while in the case of the evaluation of cerebral autoregulation the input is MAP variability and the output is the MCBF variability [18,45,46]. It is well-known that SAVR reduced baroreflex control [26]. This observation is confirmed in this study by the dramatic postoperative reduction of TFM RR-SAP in the LF and HF bands both at REST and during STAND. Signs of an impaired cardiac baroreflex were observable just before SAVR given that the expected decrease of TFM RR-SAP in the LF and HF bands during STAND compared to REST in PRE was not visible [4,43,50,51]. The post-surgery depression of the baroreflex regulation in association with the weakness of this reflex just in PRE might expose SAVR patients to possible postoperative arrhythmic events [26,27]. Cerebral autoregulation mechanisms seem to be preserved after SAVR as indicated by the negligible post-surgery changes of TFM, phase and K 2 between MCBF and MAP in the various frequency bands and this result is tightly linked to the preservation of the MAP and MCBF powers. Again, this finding could be related to the maintenance of the sympathetic control after SAVR, even though the contribution of different factors including the beneficial effect of surgery over the arterial blood flow profile cannot be dismissed. This finding suggests that SAVR patient might be more at risk of developing postoperative arrhythmic episodes than postoperative cerebrovascular adverse events.

Conclusions
The major depression of the vagal autonomic control and cardiac baroreflex after SAVR might expose SAVR subjects to a greater risk of postoperative arrhythmic adverse events. Cerebrovascular regulation seems to be less affected by post-surgery autonomic regulation derangement and this finding suggests that SAVR patients feature internal resources that might maintain the risk of postoperative cerebrovascular adverse events at the preoperative level. We remark that the direct association between cardiovascular adverse events after SAVR and the depression of vagal control and cardiac baroreflex regulation deserves specific future studies involving an adequate number of subjects and classification of type and rate of events. Given the important dispersion of the cerebrovascular variability and cerebral autoregulation indexes, individual analysis is necessary to check whether people with markers compatible with a cerebrovascular impairment might have a greater probability of developing postoperative cerebrovascular adverse events. We remark the need of increasing the power of the study by enlarging the size of the group to verify whether the inability to distinguish differences could the mere consequence of the small sample size. Future studies should assess whether indexes could indicate any crosstalk between cardiovascular and cerebrovascular regulations and could be helpful to identify risky subjects.
Supporting information S1 Appendix. Spectral and cross-spectral methods. It contains supporting information relevant to spectral and cross-spectral methods. (PDF) S1 Dataset. Database of the extracted variability markers. It contains spectral, cross-spectral, complexity indexes as well as the result of the test checking for the presence of nonlinear dynamics computed at REST and during STAND before and after SAVR. (ZIP)