The meaning of “coherent” and its quantification in coherent hemodynamics spectroscopy

We have recently introduced a new technique, coherent hemodynamics spectroscopy (CHS), which aims at characterizing a specific kind of tissue hemodynamics that feature a high level of covariation with a given physiological quantity. In this study, we carry out a detailed analysis of the significance of coherence and phase synchronization between oscillations of arterial blood pressure (ABP) and total hemoglobin concentration ([Hbt]), measured with near-infrared spectroscopy (NIRS) during a typical protocol for CHS, based on a cyclic thigh cuff occlusion and release. Even though CHS is based on a linear time invariant model between ABP (input) and NIRS measurands (outputs), for practical reasons in a typical CHS protocol, we induce finite “groups” of ABP oscillations, in which each group is characterized by a different frequency. For this reason, ABP (input) and NIRS measurands (output) are not stationary processes, and we have used wavelet coherence and phase synchronization index (PSI), as a metric of coherence and phase synchronization, respectively. PSI was calculated by using both the wavelet cross spectrum and the Hilbert transform. We have also used linear coherence (which requires stationary process) for comparison with wavelet coherence. The method of surrogate data is used to find critical values for the significance of covariation between ABP and [Hbt]. Because we have found similar critical values for wavelet coherence and PSI by using five of the most used methods of surrogate data, we propose to use the data-independent Gaussian random numbers (GRNs), for CHS. By using wavelet coherence and wavelet cross spectrum, and GRNs as surrogate data, we have found the same results for the significance of coherence and phase synchronization between ABP and [Hbt]: on a total set of 20 periods of cuff oscillations, we have found 17 coherent oscillations and 17 phase synchronous oscillations. Phase synchronization assessed with Hilbert transform yielded similar results with 14 phase synchronous oscillations. Linear coherence and wavelet coherence overall yielded similar number of significant values. We discuss possible reasons for this result. Despite the similarity of linear and wavelet coherence, we argue that wavelet coherence is preferable, especially if one wants to use baseline spontaneous oscillations, in which phase locking and coherence between signals might be only temporary.


Introduction
Recently, our group proposed a novel method, named coherent hemodynamics spectroscopy (CHS), for studying microvascular and microcirculation integrity by inducing controlled hemodynamic changes in living organisms. 1 These hemodynamic changes may be associated with perturbations in the arterial blood pressure (ABP) that, in human studies, can be induced by a forcing mechanism such as paced breathing, 2,3 cyclic occlusions and releases of pneumatic thigh cuffs, 4 head-up tilt protocols, 5 or squat-stand maneuvers. 6 For brain studies, the ensuing changes in cerebral blood volume and blood flow can be measured by a number of techniques, including functional magnetic resonance imaging, near-infrared spectroscopy (NIRS), transcranial Doppler ultrasound, and diffuse correlation spectroscopy. 7 The basic hypothesis of CHS is that tissue with healthy microvasculature and perfusion differs from diseased tissue in terms of the dynamic relationship between blood flow and blood volume (or their covariates such as the tissue concentrations of oxy-[HbO 2 ], deoxy- [Hb], and total hemoglobin [Hbt]) in response to hemodynamic perturbations such as those induced by systemic ABP changes. In particular, when perturbations in ABP are the driving force of the observed hemodynamic changes, one can study the phenomenon of cerebral autoregulation (CA), which is critically important for brain health. 7 The idea of studying the covariation of [Hb] and [HbO 2 ] with NIRS to obtain diagnostic or functional information is not new. For example, spontaneous cerebral co-variations of [Hb] and [HbO 2 ] were characterized in infants 8,9 and adults, [10][11][12] and the phase relationship of paced-breathing-induced oscillations in cerebral [Hb] and [HbO 2 ] was associated with the level of CA. 13 Furthermore, covarying representations of cerebral [Hb] and [HbO 2 ] were proposed for functional brain studies. [14][15][16] The novelty of the CHS technique is that it targets a specific kind of cerebral hemodynamics, namely, those that feature a high level of coherence and phase synchronization with a specific driving physiological process. CHS is not only based on the characterization of multi-frequency oscillations or dynamic transients of cerebral [Hb] and [HbO 2 ], but it also translates such frequency-or time-resolved dynamic information into physiological quantities on the basis of a CHS mathematical model that treats the cerebral microvasculature as a linear, time-invariant system. 1 The CHS model shows that oxy-and deoxy-hemoglobin concentrations oscillations, in particular their phase and amplitude relationships, depend on ABP oscillations (and associated blood volume and flow oscillations) through six physiological parameters, including the capillary and venous blood transit times, relative blood volumes in the arterial, capillary, and venous compartments, a measure of CA, etc. These parameters are fitted for in an inversion procedure using our CHS model as a forward solver. 17 The strengths of our CHS model and its current limitations with respect to others proposed in the literature are described concisely in our previous publications. 7,18 The CHS treatment of the microvasculature as a linear time-invariant system requires a high level of coherence or covariation between the measured hemoglobin concentration dynamics and the physiological changes of blood volume, blood flow, and oxygen consumption in the investigated tissue. In fact, when a process is stationary, low coherence levels between measured hemoglobin concentrations and underlying physiological processes may reflect a nonlinear relationship, a dependence on additional variables not considered in the CHS model, or contributions from confounds or noise. 19 This work considers the coherence or covariation between cerebral total hemoglobin concentration [HbT] and systemic ABP, which are directly relevant for CHS, and presents methods to define proper metrics and their critical values for significance.
We note that the study of the significance of metrics for coherence or covariation between signals is a prerequisite for studying CA. For example, in transcranial Doppler, it is required high coherence between ABP and blood flow velocity (BFV) in the middle cerebral artery. 20 Some NIRS published studies have used ABP and the difference [HbO 2 ]- [Hb], which has been considered a surrogate of cerebral blood flow. 21 The importance of coherence between ABP and cerebral concentrations of oxy-and deoxyhemoglobin was also reported in a transfer function analysis of NIRS signals. 10 In all of these studies, a coherence greater than a threshold value of about 0.5 was considered to be significant. More generally, a threshold coherence value for significance should result from a statistical test of a null hypothesis of zero coherence, which depends on the degrees of freedom associated with the specific case of data collection and analysis. Following this statistical procedure, lower threshold values in the range 0.1-0.3 have been reported. 19,22,23 The thresholds of significance for coherence reported in these works were demonstrated by Koopmans 24 for the case of stationary Gaussian processes and nonoverlapping segments used for the Welch's periodogram method. A step forward in coherence estimation was achieved by Gallet and Juliene, 25 in which the authors derived a formula valid also for overlapping segments according to Welch's method. We also point out the work of Faes et al., 26 in which the authors used surrogate data for defining a threshold of significance of coherence in cardiovascular variability analysis and showed that such a threshold could be frequency-dependent.
In this work, we find critical values of significance for coherence and phase synchronization of ABP and [HbT] under typical conditions of CHS measurements, and we discuss their relevance for CHS. We have used [Hbt] instead of [HbO] and [Hb], because usually it is the most robust NIRS measurand. However, the results of this study can be applied also for the covariation between [HbO] and ABP, and [Hb] and ABP. These two NIRS measurands would be used in a CA study. A key assumption of the CHS technique and model is that the dynamics of the concentrations of hemoglobin species be mainly driven by a single physiological or functional process, say changes in ABP, in which case they are expected to feature a high level of coherence and covariation with such process. In reality, in addition to ABP, there are a number of factors that can affect the concentrations of oxy-and deoxyhemoglobin in brain tissue: for example, fluctuations in blood CO 2 levels, changes in oxygen consumption during brain activation, local vasomotion, spontaneous changes in respiratory and heart rates (which do lead to associated changes in blood pressure), Mayer waves, etc. Therefore, in order to minimize the effects of other confounds, we have proposed a protocol for CHS that uses cyclic thigh cuff occlusion and release maneuvers at multiple frequencies, in which the subjects are resting and not engaged in any mental or physical challenge. This protocol of cyclic inflation and deflation of pneumatic thigh cuffs is known to induce cyclic oscillations in ABP. 4 However, spontaneous oscillations in ABP may also contribute to the observed dynamics of hemoglobin concentration. In transcranial Doppler studies, in which the autoregulation in the macrocirculation is studied, no distinction is made between fluctuations in BFV that are caused by spontaneous or induced oscillations in ABP (as long as the amplitude of blood pressure oscillations is above a certain threshold; see method of coherent averaging 27 ). Similarly, in our study, we have not investigated whether blood pressure oscillations are induced by thigh cuff oscillations (this could be done by using the methods discussed later) but we only have focussed on a metric of coherence or covariation between oscillations in ABP and in cerebral hemoglobin concentration.
As a metric of covariation, one could use the magnitude square coherence (linear coherence) and test whether its value at the frequency of cuff pressure oscillations is above a certain threshold of significance, as discussed in the literature quoted earlier. However, we must assume that the whole process under observation is stationary. In our study, we consider a protocol that involves an initial baseline acquisition for 2 min followed by five 2 min time windows in each of which we induce cuff pressure oscillations at different frequencies in the range 0.045-0.083 Hz. Twenty seconds of baseline separates two consecutive time windows. Because of the nature of our protocol, the time series of blood pressure and total hemoglobin is nonstationary (meaning that they have time-varying frequency components), and therefore we have considered more appropriate the wavelet coherence and the method of analytic signal based on the Hilbert transform to define metrics of covariation between ABP and total hemoglobin oscillations. Wavelet coherence was used to define a metric for the timefrequency-resolved value of coherence between ABP and total hemoglobin concentration. Since the coherence between two signals is affected by the covariation of both amplitude and phase of the signals, [28][29][30] we also define a metric only based on phase covariation (or phase synchronization) of the two signals. In fact, as a result of a strong amplitude covariation, the coherence of two signals may be significant even when the phase covariation is not. 29 In summary, in this work, we have designed an ad hoc method, tailored to our CHS protocol, in order to investigate the level of coherence and phase synchronization between ABP and cerebral total hemoglobin concentration. Since our protocol induces nonstationary oscillations of ABP and NIRS measurands, we used two metrics suitable for nonstationary data to determine whether total hemoglobin oscillations were directly associated with ABP oscillations: (a) wavelet coherence and (b)phase synchronization index (PSI) which are suitable for nonstationary data. For PSI, the phase between two signals was defined by using the wavelet cross spectrum and also the method of analytic signal, based on the Hilbert transform. We have also used linear coherence for comparison with wavelet coherence. The thresholds or critical values of significance for these two metrics were defined by the method of surrogate data, which involve the derivation of the distributions of the two metrics under a null hypothesis of zero coherence or zero phase synchronization, respectively. Among surrogate data, we have also tested Gaussian random numbers (GRNs), because of the practicality to derive a threshold for our metrics that is only dependent on the protocol and data processing method, but that is not specific to the particular experimental data set.
Random numbers have already been used in the literature as surrogate data for statistical inference. 31 We observe that published NIRS studies on functional connectivity 32-36 and autoregulation [37][38][39][40][41] have used the concept of coherence and/or phase synchronization between signals, by using continuous wavelet transform, wavelet coherence, and the analytic signal methods. In a few studies, [32][33][34]36,37,41 a method of surrogate data was used in order to assess the significance of the signals covariation. Other NIRS studies in which the method of surrogate data was used for statistical inference are: Bernjak et al., 42 in which blood flow and oxygen saturation were studied on the skin of human subjects and Bu et al., 43 in which the authors studied the effect of sleep deprivation. To the best of our knowledge, this is the first NIRS study where the question of covariation between signals has been addressed more critically by using different metrics of covariation and different methods for the estimation of threshold values for significance of those metrics.

Experimental protocol and data acquisition in vivo
The NIRS measurements on human subjects were performed with a commercial frequencydomain NIRS instrument (OxiplexTS, ISS Inc., Champaign, IL, USA). Optical probes connected to the spectrometer delivered light at two wavelengths, 690 nm and 830 nm, at a source-detector distance of 35 mm. The probe was placed against the left side of the subject's forehead, to access brain tissue in the prefrontal cortex, and secured with a flexible headband. Continuous ABP was recorded with a finger plethysmography monitoring system (NIBP100D, BIOPAC Systems, Inc., Goleta, CA, USA). Pneumatic thigh cuffs were wrapped around both of the subject's thighs and connected to an automated cuff inflation system (E-20 Rapid Cuff Inflation System, D.E. Hokanson, Inc., Bellevue, WA, USA). The air pressure in the thigh cuffs was continuously monitored with a digital manometer (Series 626 Pressure Transmitter, Dwyer Instruments, Inc., Michigan City, IN, USA). Analog outputs of the ABP monitor and the thigh cuff pressure monitor were fed to auxiliary inputs of the NIRS instrument for concurrent recordings with the NIRS data. All signals were sampled synchronously at an acquisition rate of 12.5 Hz. The protocol consisted of 2 min of baseline acquisition followed by five groups of thigh cuff occlusion and release at the frequencies: 0.0454, 0.0555, 0.0625, 0.0714, and 0.0833 Hz, corresponding to oscillation periods of 22, 18, 16, 14, and 12 s, respectively. Each group consisted of six consecutive cuff inflation/deflation cycles followed by 20 s of rest (where the cuffs stayed deflated). At the end of the last cycle, 1 min of recovery was collected. The order of the cuff frequencies was scrambled in order to have a minimum overlap among the induced frequencies between consecutive groups of oscillations (see Sec. 2.4.2). Specifically, the sequence of frequencies of cuff oscillations was: 0.0625, 0.0833, 0.0454, 0.0714, and 0.0555 Hz. The maximum cuff pressure during the legs occlusion was 200 mmHg. The Tufts University Institutional Review Board approved the experimental protocol, and the subjects provided their informed consent. In this study, we show the results obtained on four subjects.

Wavelet coherence and analytic signal methods
In this section, we describe the wavelet coherence and the analytic signal approach, which we have used in the analysis of our data to obtain measures of coherence (in addition to linear coherence by Fourier transformation) and phase synchronization, respectively. Given a real function x(t), the continuous wavelet transform C x (cwt) is defined as: where φ(t) is the mother wavelet, and a and b are the scale and location parameters, respectively. In Eq. (1), the * and ⊗ symbols denote the complex conjugate and the convolution operator, respectively. 44 We can describe intuitively the wavelet coefficients C x (a, b) as measuring the local temporal similarity of the original function x(t) and a timescaled version of the original mother wavelet φ(t). A typical mother wavelet used in cwt is the complex Morlet wavelet, φ(t):  Note that it can be misleading to associate one frequency to a scaled wavelet. In fact, due to its time localization, a wavelet is characterized by a spectrum of frequencies and f represents only the value where the maximum of the absolute value of its Fourier transform (FT) is found. Equation (3) is the formula that allows one to shift between scale and characteristic frequency. Time scaling of the mother wavelet is therefore achieved by shrinking (a < 1) or dilating (a > 1) the wavelet in time, which corresponds to probing a signal locally in a frequency band centered at higher or lower frequencies, respectively.
Given two functions x(t) and y(t) the wavelet cross spectrum C x,y and wavelet coherence Coh x,y are defined by: where S is a smoothing operator in time and scale. Here we have reported the formulas provided in MATLAB documentation (https://www.math-works.com/help/wavelet/ref/ wcoherence.html). The time-scale resolved phase difference between the two signals is defined as: A different method to define the phase of a time-dependent function is the analytic signal method. 45 Given a function x(t), its Hilbert transform (HT) is defined as H x (t): where p.v. denotes the Cauchy principal value of the integral. One important property of the HT is that it shifts the phase of the Fourier components by ±π/2: where the tilde (~) denotes the FT and sgn(ω) is the signum function. The analytic signal of a function x(t) is defined as: and it has the property that its FT contains only positive frequencies. However, the definition of analytic signal by itself does not guarantee a unique definition of instantaneous phase and frequency of a signal, the reason being that the analytic signal may contain a broad spectrum of positive frequencies. The method of analytic signal leads to a unique (noncontradictory) definition of instantaneous phase and frequency of a signal only if the Bedrosian's product theorem holds true, 45 which requires that a signal is defined in a "narrow" frequency band. If this requirement is not met, paradoxical results for the instantaneous phase and frequency of a signal are found 46 (pp. 913-914). Therefore, for a correct definition of instantaneous phase by the analytic signal method, first a narrow band pass filter centered at the frequency of interest is applied to the signal x(t) and a filtered signal x F (t) is obtained. Afterwards, the analytic signal is defined by: where H x F is the Hilbert transform applied to the filtered signal. The instantaneous phase is defined as: . (11) Both wavelet transform and the analytic signal method allow one to define an instantaneous phase difference between two signals which are particularly useful under nonstationary conditions. Usually, given a time interval in which the phase difference of two signals is studied, both methods yield a distribution of phase values (one phase value for each time point), and the question arises if the distribution is consistent with the hypothesis of a constant, well-defined phase difference between the signals. If the distribution is peaked around a certain phase value, we could argue in favor of a phase covariation or synchronization of the two signals. On the contrary, if the phase distribution is uniform or rather broad, we conclude that there is no phase covariation or synchronization of the signals. In this study, we have used the concept of PSI which is based on the definition of entropy of a random variable. 47 Given a discrete random variable with N possible values, the definition of Shannon entropy is: where P i is the probability of the ith value of the random variable. The PSI is defined as: where E max = ln(N). Strictly speaking, the Shannon entropy of a random variable is defined by using log 2 , since the maximum entropy coincides with the number of bits necessary to code the maximum information, i.e., when all the values are equiprobable. In general, we can easily demonstrate that, regardless of the base for the logarithm d (d > 1), the maximum entropy occurs when all the outcomes of the random variable are equiprobable and it coincides with log d (N) On the contrary, the entropy of a random variable is zero only when one outcome is certain (i.e., P i = 1; P j = 0, j ≠ i). From the definition of entropy, it follows that the PSI can take values in the range 0-1 and, specifically, PSI = 1 when E = 0, and PSI = 0 when E = E max . Other metrics for phase synchronization are also found in the literature.
We mention a metric inversely related to the standard deviation of an angular variable, 44,48 an index based on conditional probability, 47 and an index based on mutual information. 49 We note that different metrics may yield different synchronization values because they are sensitive to different features of the phase distribution of the two signals. For example, the PSI of the relative phase of two signals can have a high value also for a multimodal distribution, i.e., having distinct narrow peaks. However, in this case, the distribution might have a large standard deviation and therefore the related synchronization index might not be significant. More complex measures of synchronization have been considered in the field of electroencephalography (EEG). 50

Method of surrogate data
We used surrogate data to estimate the statistical distribution that our experimental data would follow under conditions of no coherence or no phase synchronization. The method of surrogate data was initially proposed for assessing conditions of nonlinear dynamics in both univariate and multivariate time series. 49,[51][52][53][54][55][56][57] It has also been used for testing significance for linear coherence 26,58 and phase synchronization of signals. 49,[59][60][61][62] Surrogate data are simulated data that are generated under a null hypothesis of zero coherence or zero synchronization between two signals. There are different ways to generate pairs of surrogate data having zero coherence or zero phase synchronization. For example, we could assume that the surrogate data are derived from two independent white noise processes, or instead that they are derived from two independent linear Gaussian processes. Therefore, the null hypothesis of zero coherence or zero phase synchronization will be more specifically linked to the underlying stochastic processes the data are generated from. The method to generate surrogate data from different processes is based on retaining different features of the experimental signals, but in all the cases the temporal structure and/or synchronization present in these signals is destroyed by following a process of randomization. The process of randomization is repeated, and the metrics (coherence and PSI in this study) are recalculated for each pair of surrogate data to generate a distribution of the metrics under the null hypothesis. Afterwards, the values of coherence and PSI obtained for the actual experimental signals are compared with a threshold value derived from the distribution obtained under the null hypothesis (usually the 95th percentile, which corresponds to a typical error of type I, α of 0.05).
Some of the most common surrogate data methods are the following, where we used the nomenclature of Paluš. 49 (1) IID1 surrogates are realization of independent identically distributed (IID) white noise stochastic processes which are obtained by scrambling the order of the original signals with two independent permutations for the samples of each pair of surrogate data. The null hypothesis presented by the IID1 surrogates is that the original signals result from independent white noise processes, so that any temporal structure (and synchronization) present in the original signal is destroyed, and only the distribution of the values is preserved.
(2) IID2 surrogates are obtained by using the same permutation to scramble the order of the original signals. In this case, the null hypothesis is that the original signals are mutually dependent white noise processes, where the distribution of the values and the Pearson correlation coefficient between the original signals are preserved during the process of randomization.
(3) FT1 surrogates are obtained by randomizing the phase of the FT of the signals. The two samples in a pair of surrogate data are obtained by using two independent random sequences for the phase of their FTs. In this case, the spectra (or the autocorrelation functions) of the signals are preserved, and the surrogate data realize the null hypothesis of two independent linear Gaussian stochastic processes that asynchronously oscillate at the same frequencies of the original signals. 49 (4) FT2 surrogates are obtained in a similar manner as FT1 surrogates but the phases are randomized by using the same sequence in the samples of a pair of surrogate data. In this case, each randomization preserves not only the auto spectra but also the cross spectra. The null hypothesis in this case is that of two linear stochastic processes showing a linear synchronization, which is usually weaker than that from nonlinear stochastic processes. 49 (5) AAFT1 (amplitude-adjusted FT) surrogate data preserve both the auto spectra and the distribution of values of the signals (which is not preserved by FT1). The null hypothesis in this case is that the signals are a monotonic (static) nonlinear transformation of linear Gaussian processes that oscillate asynchronously at the same frequencies of the original signals. 51 (6) Similar to FT2, one can define AAFT2, which preserves auto and cross spectra, and distribution of the signals.
(7) GRN (Gaussian random numbers) surrogate data are a sequence of independent Gaussian random numbers. In this case, no feature of the original signals is preserved and these surrogate data present the null hypothesis that the signals result from independent Gaussian white noise processes. Several of these surrogate data were discussed in the work of Paluš, 49 others like AAFT1 are discussed in the work of Theiler et al., 51 whereas GRN was adopted in the work of Xu et al. 31 The method AAFT2, to the best of our knowledge, was not discussed previously. Because FT2 (and AAFT2) preserves the cross spectrum of the original data, it was argued that these surrogate data preserve some degrees of linear synchronization and that the rejection of FT2 implies the detection of nonlinearity in the abp T phase synchronization. 49 In this work, we have not investigated this point (i.e., if the type of synchronization between abp and T is linear or nonlinear); therefore, we have applied the IID1, IID2, FT1, AAFT1, and GRN surrogate data methods.

Data processing
2.4.1. Optical data and ABP-Oxyand deoxyhemoglobin concentration changes were calculated from the measured optical intensity changes at two wavelengths (690 nm and 830 nm) by the modified Beer-Lambert law 63,64 with differential pathlength factors (DPFs) of 7.8 at 690 nm and 7.1 at 830 nm. These values are within a range of values that can be found in the literature for brain tissue. Even though the real DPFs may be different from the ones we assumed, the results of this study (about the significance of interdependence between two signals) are not expected to be significantly affected. The changes of intensity were calculated with respect to a baseline value which was defined as the average intensity value during baseline (first 2 min of data acquisition). Total hemoglobin concentration changes (T) were calculated as the sum of oxy-and deoxyhemoglobin concentration changes. ABP was de-meaned and normalized by its mean value (average value during baseline, ABP 0 ). The normalized ABP which will be used for studying the covariation with T is identified with lower case letters (abp) and defined as

Wavelet coherence and analytic signal-
We used the MATLAB built-in function "wcoherence" for the calculation of wavelet coherence and cross spectrum between ABP and T. We used standard settings of the function, including the complex Morlet wavelet as mother wavelet (https://www.mathworks.com/help/wavelet/ref/wcoherence.html). The output of the function was the time-frequency resolved wavelet coherence and cross spectrum between abp and T. More precisely, the coherence was a matrix of real numbers (Coh T, abp WT (a, b); the superscript WT stands for wavelet transform) in the range [0, 1] with 145 rows and a number of columns equal to the number of time points of the original abp and T signals. We remind that the parameter a is the scale and b is time. Each row corresponded to a different scale, therefore to a different characteristic frequency in the range 0.001457-5.968 Hz. The frequencies are obtained by defining the scales using 12 octaves (this number depends on the length of our data) and 12 voices per octave (by default). More specifically, the scales are defined by a = 2 , where in our case i o = 0, 1, …, 12, i v = 0, 1, …, 11, and N v = 12. The characteristic frequencies are derived from the scales by the formula: f = f c /a where f c = 6/(2π). F s /2, where F s is the sampling frequency. The cross spectrum (C T, abp WT (a, b)) was a matrix of complex numbers, having the same size as the previous one.
From the cross spectrum matrix, it was possible to define the time-frequency resolved phase difference between signals (Eq. (6)). In order to assign a value and an error to the coherence and phase difference between abp and T in each group of six thigh cuff oscillations, we used the following procedure: (a) first we identified a range of frequencies associated with the group of six cuff oscillations. In the assumption of sinusoidal oscillations, one group has the absolute value of the FT with a main lobe defined by the frequencies (c) we defined two temporal arrays obtained by averaging the phase differences ( Δ φ T, abp WT (b) ) and coherence ( Coh T, abp WT (b) ) across the scales in the main lobe. We used circular statistics for averaging the phase differences 65 ; (d) we identified the time range associated with the beginning and the end of a group of six cuff oscillations, namely (t b , t e ), and calculated the average value and standard deviation of (13)), we arranged in a circular histogram a subset of the values Δ φ T, abp WT (b) , those defined in the interval (t b , t e ), by dividing the 360° angle into 72 angular intervals of 5° and calculating the frequency of occurrence as estimate of P i (see Eq. (12)) in each interval.
PSI T, abp WT was calculated according to Eq. (13). Note that the calculation of PSI T, abp WT was not associated with an error, since the data are used only one time for the estimation of one value of PSI T, abp WT .
The analytic signals of abp and T were defined by first applying a narrow band pass filter to the original signals. More precisely, we generated the filter's coefficient by using a Parks McClellan algorithm corresponding to a pass band of 0.01 Hz, a stop band of 0.02 Hz and centered at f cuff . Note that the stop band had similar width of the main lobe of the cuff frequencies, which was in the range 0.015-0.028 Hz. The ripple in the pass band and the attenuation between the stop band and the pass band were less than 0.05. The coefficients of the filter were used as input of a zero-phase filtering procedure ("filtfilt" function in MATLAB). The phase between T and abp was defined as: Arg( ) where Z T F (t) and Z abp F (t) are the analytic signals defined after filtering T and abp (Eq. (10)). An average phase difference and standard deviation were defined by using circular statistics in the time range where the cuff oscillations occurred (t b , t e ), therefore obtaining Δ φ T, apb where the superscript HT stands for Hilbert transform. We also calculated PSI T, abp HT in a similar manner as for PSI T, abp WT .

Example of raw data and phase difference plots
An example of raw data is shown in Fig. 1  As it is visible from the raw abp, T, and CP signals, the cuff oscillations drive the amplitude and entrainment of abp and T oscillations at its own frequency. These oscillations in abp and T were not present (or present to a much lesser extent) during baseline. For these reasons, the signals obtained in our protocol are strictly nonstationary, so that nonstationary methods (like wavelet and Hilbert transforms or short time FT (STFT)) are in principle more appropriate. For example, from Fig. 1

Author Manuscript
Author Manuscript

Author Manuscript
Author Manuscript value Coh T, abp FT = 0.90 ± 0.05. The significance or non-significance of coherence or PSI was true regardless of the null hypothesis for surrogate data. We have found that Fourier, Hilbert, and wavelet-based methods always agreed within errors for the calculation of the phase differences between T and abp. One example is shown in Fig. 2 for the same subject of Fig.   1. On the left panel are plotted the average phase differences and standard deviations calculated at the cuff frequencies with the three methods. For the cross power spectral density (FT), these quantities are obtained by considering the frequencies across the main lobe of the cuff frequencies (Sec. 2.4.3), whereas for the analytic signal (HT) and wavelet cross spectrum (WT) methods, after averaging across the frequencies, the final averages and standard deviations were calculated in the time intervals where cuff oscillations occurred (Sec. 2.4.2). On the right panel, the same time-frequency calculation of averages and standard deviations was carried out for the analytic signal method (HT) and wavelet cross spectrum (WT) during baseline. As one can see from these two nonstationary methods, the standard deviations are much smaller during cuff oscillations than during baseline, which can be interpreted as enhanced phase synchronization or entrainment between T and abp driven by the cuff. We remind that one of the synchronization indices defined in the literature (not used in this work) is inversely related to the standard deviation of the phase difference distribution. 44,48 Therefore, the content of Fig. 2 is in agreement with the different values of coherence and PSI found during cuff oscillations and baseline (Fig. 1). Also, from Fig. 2 (right panel), we can see that the phase differences calculated with the analytic signal and wavelet methods do not always coincide (at least when the synchronization of the signals is poor). We will comment on this fact in the discussions. These results confirm that Fourier methods tell us only if abp and T are overall coherent or entrained, whereas analytic signal and wavelet methods tell us also when this happens.

Significance of phase synchronization and coherence
A comparison between real data and threshold values for coherence and phase synchronization is presented in Fig. 3 )). In each panel, the x axis is defined as in Fig. 2. The threshold values for these parameters were calculated by using five different methods to generate surrogate data: (1) Gaussian random numbers (GRN, diamond symbol); (2) independent identically distributed white noise stochastic processes (IID1, square symbol); (3) mutually dependent white noise processes (IID2, cross symbol); (4) linear stochastic processes that asynchronously oscillate at the same frequencies of the original signals (FT1, circle symbol); and (5) monotonic (static) nonlinear transformation of linear Gaussian processes that oscillate asynchronously at the same frequencies of the original signals (AAFT1, asterisk symbol).
The values of these parameters calculated for the experimental data are also shown (data, triangle symbol). We remind that the PSI values for the real data are calculated only one time, and therefore no error is associated to them. In this work, we have not applied statistical tests between experimental data and threshold values; rather, we define an experimental data value as significant if it is not in the error range of the corresponding threshold value. If an experimental data value has an associated error range (as is the case for Coh T, abp FT and Coh T, abp WT ), we define it as significant if its range does not overlap the range of the threshold value. As we can see from Fig. 3 (subject 1), almost all of the four covariation parameters are significant regardless of the null hypothesis of surrogate data.
The only exception is PSI T, abp WT (Fig. 3(d)), which is not significant at f cuff = 0.056 Hz.
However, this value becomes significant if we choose IID1 and GRN as surrogate (we note that for this case the symbol size is larger than the error bars).
In Fig. 4 (subject 2), we see some discrepancies about the significance of the parameters of covariation. For example, if IID2 is chosen for threshold, Coh T, abp FT ( Fig. 4(a)) has no significant values, whereas Coh T, abp WT ( Fig. 4(b)) has two significant values (at f cuff = 0.071 Hz and 0.083 Hz). PSI T, abp HT (Fig. 4(c)) has one significant value (at f cuff = 0.083 Hz),

Surrogate data thresholds across subjects
In this section, we show a comparison of the threshold values obtained by different surrogate data across subjects. If possible, one would want to disentangle the threshold values from the particular data set collected on a subject, so that a single threshold value may be used. In this sense, the threshold values obtained with GRN are not data-dependent, but depend only on the details of the protocol and on the details of the calculation method of the metric of covariation. Once those are fixed, the thresholds obtained with GRN are also fixed and could be used in a lookup In Fig. 7, we note that the variance across subjects with GRN and IID1 seems to be negligible, but it is not the case for IID2 and FT1. For IID2, the values scale according to the Pearson correlation coefficient between abp and T: subject 2 (square symbols in Fig. 7) had the lowest value of 0.33, whereas subject 4 (circle symbols in Fig. 7) had the highest value of 0.59. This discrepancy of the threshold values between subjects obtained with surrogate data IID2 and FT1 is also confirmed by repeating the calculations with higher statistics. We used 1,000 pairs of surrogate data and we extracted Coh T, abp WT for each set; this was repeated three times and we computed Coh T, abp WT ± σ Coh WT . For this calculation, we also chose a more refined partition of the interval (0,1) which was divided into 50 sections.
The standard deviation σ Coh WT is defined as the maximum between the one obtained from the three runs and 0.02 (width of each section). The results are shown in Fig. 8, in which we can see that subject-dependent thresholds can happen at all frequencies (IID2; left panel) or only at specific frequencies (FT1; right panel). In Fig. 9, results are shown for subject 3 (diamond symbols) and subject 4 (square symbols), and they were calculated by using better statistics (1,000 pairs of surrogate data) as in Fig. 8.
As one can see in Fig. 9, also for other parameters of covariation, the threshold values calculated from surrogate data may depend on the subject; therefore, they should be recalculated for better accuracy when a new data set on a different subject is collected.
Finally, for all the four parameters of covariation the widest variance across subjects is found by using IID2 surrogate data (results not shown).

How many frequencies are induced?
Last, we discuss the effectiveness of our protocol to induce coherent oscillations of T and abp. One direct way (not used in this work) is to assess the coherence/synchronization between abp and cuff signal and T and cuff signal separately (at a given f cuff ), by using one of the surrogate data described before. One indirect way is to compare the coherence/ synchronization between abp and T during cuff oscillations and during baseline. If the coherence/phase synchronization is higher during cuff oscillations than baseline (at a given f cuff ), we say that a particular frequency f cuff was induced by the cuff operation. This point can only be studied with nonstationary methods; therefore, we used three of the four parameters of covariations, namely, Coh T, abp WT , PSI T, abp HT , and PSI T, abp WT . These parameters were calculated during the periods of cuff oscillations and also during baseline, and their significance was decided by assuming the threshold values calculated with FT1 surrogate data (Figs. 3-6). The results are summarized in Table 2

Discussion
In this section, we address some questions arising from the choice of surrogate data and about the significance of the four parameters of covariation considered by us.
In our data analysis, different surrogate data yield comparable threshold values for significant covariation. However, on average, GRN and IID1 yield slightly lower threshold values, followed by FT1 and AAFT1, and last by IID2. Except GRN, this hierarchy of threshold values according to different null hypothesis reflects the information of the original data which is preserved after different processes of randomization. For IID1, only the distribution of the original data is preserved in the surrogate data, but oscillations and any kind of temporal feature are absent. For FT1, the auto spectra of the original data are preserved, and therefore the surrogate data are independent realization of the original data with de-synchronized oscillations (if they were present in the original data). AAFT1 preserves both auto spectra and the distribution of the original data. It is expected that the threshold values calculated with FT1 and AAFT1 are greater than those found with GRN and IID1 and that the difference is related to the broadness of the spectrum of the original signals. In fact, if two signals are narrow band (ideally only one or few Fourier components), the randomization of the Fourier components will not destroy (or destroy only partly) the coherence/synchronization of the signals calculated at the center frequency of the band. On the contrary, GRN is characterized by a constant power spectrum (independent on frequency) where the Fourier components have unrelated phases, therefore yielding lower values of PSI. Therefore, for the case of a narrow band spectrum, the threshold values obtained with FT1 would be progressively larger (as the band width of the spectrum is decreased) than those obtained with GRN. We verified this point by using finite sinusoidal oscillations as test signals (results not shown). Similar results were found by Faes et al. 26 by comparing IID1 and FT1 for studying the threshold of linear coherence in cardiovascular variability analysis (in our study we have found the same thresholds with IID1 and GRN). The authors concluded that the advantage of using FT1 (for controlling precisely error type I) was clear only for narrow band signals. While for broad band signals, the two surrogate data methods could be used interchangeably.
The method IID2 yields threshold values greater than those obtained with FT1 and AAFT1. This result was not found in the work of Paluš 49 ; on the contrary, the author found similar threshold values for IID1 and IID2 when studying the synchronization of the number of solar spots and the temperature in Prague, as function of the calendar year. Since in our study ABP and T are measured independently, any correlation of the signals cannot be the result of correlated noise, but only due to a true covariation of the signals. Therefore, we argue that the process of randomization that preserves the Pearson correlation coefficient might be too strict and should not be chosen. In fact, in the limit case in which the Pearson coefficient is "1" (one signal is the scaled version of the other), IID2 yields a threshold value of "1", causing an error type II of 100%.
The practical importance of using GRN as surrogate data is that the threshold values depend only on the details of the protocol and data acquisition (sampling frequency, cuff frequencies, time span of cuff oscillations) and data processing (filter bandwidth for the HT, mother wavelet choice for WT) but not on the particular data collected on a subject. with FT1. The reason for comparing GRN and FT1 is because they yield (on average) slightly lower and slightly higher threshold values (Fig. 7), respectively. These results suggest that GRN surrogate data can be used to obtain reliable threshold values of coherence and phase synchronization in a typical protocol used for CHS. We have also carried out PSI calculation by using the MATLAB function "cwt" (continuous wavelet transform) by using both Morse wavelet and complex Morlet wavelet (having f b = 2, f c = 6/(2π)), obtaining comparable threshold values by using different surrogate methods as those found with the function "wcoherence". significance viewpoint, even though the two metrics did not always agree, they shared a similar number of significantly induced oscillations for the four subjects investigated in this study (Table 1). However, we remind that the meaning of the significance is different for the two metrics: a significant Coh T, abp FT at a frequency f cuff means that during the entire experiment the two signals were overall coherent at that frequency; in contrast, a significant Coh T, abp WT at f cuff means (at least for this study) that the two signals were coherent during the time period when the cuff was oscillating at the frequency f cuff (note that Coh T, abp WT can be studied, in principle, in any time interval). One possible explanation of the similarity of these two metrics is that our particular protocol, due to the wide overlapping of frequency bands induced by the cuff oscillations (at least for some cuff oscillation frequencies), can be considered weakly nonstationary. In fact, by using  (Table 1). Given the relevance of time-resolved measures of coherence, for example, to identify coherent hemodynamics oscillations at rest (when no ABP perturbations are induced), we favor measures of coherence based on the wavelet transform rather than the FT for CHS. However, the nonstationary approach of the STFT is also suitable for CHS.
Another interesting question we want to address is about the relationship between coherence and phase covariation (or synchronization). There is a large amount of literature (especially in the field of EEG and magnetoencephalography (MEG)) pointing out that coherence is affected by both phase synchronization and covariation of amplitudes. [28][29][30]66 For example, Srinath and Ray 29 found that even when the phase relationship between two EEG channels is random, a strong amplitude covariation may introduce a significant coherence. Also, in an MEG study by Tass et al., 47 the authors found a significant coherence between channels without a detection of significant phase synchronization. These results, observed in two different research fields, point to the conclusion that a significant phase synchronization might be more difficult to achieve than a significant coherence. In our study, we can compare the coherence and phase synchronization metrics by releasing the assumption of  (1)), we can see that the WT is a convolution and that at each scale (i.e., each characteristic frequency, f = f c /a), the scaled mother wavelet acts as a filter (at the frequencies centered around f) to the entire signal. We can prove that for any scaled wavelet of the Morlet mother wavelet Δf/f (where Δf is the bandwidth of the scaled wavelet) is a constant. Therefore, for the wavelet "filter" Δf is directly proportional to f ≅ f cuff and Δf · Δt is constant for different f cuff . About the significance, if one compares the two PSI metrics across the subjects (Table 1)  for FT1, respectively. In the work of Bruns, 30 the author presented theoretical arguments to infer the similarity among STFT, wavelet transform, and Hilbert transform. Bruns argued that the differences between the three non-stationary methods could be narrowed down to the choice of a window function in the three approaches, and if one chose the same window function there would be no difference among these methods. The author, after choosing similar window functions for STFT, WT, and HT, showed a qualitative comparison of the spectral amplitudes of some EEG data with the three methods. The subject of the similarity of phase synchrony with WT and HT was addressed, at least partly in the work of Le Van Quyen et al. 28 In the work of Le Van Quyen et al., the authors presented a figure (Fig. 7(d)) in which the number of significant synchronous events found with HT and WT was comparable; however, it is not clear whether they applied the same metric or two different metrics to calculate synchronicity of events. In their work, they described the metric of phase locking value (PLV) and single-trial PLV (SPLV) when using wavelet transform. These metrics are related to angular standard deviation of the phase difference between two signals. Instead, they described PSI and mutual information when using Hilbert transform. We have already discussed in the introduction that the metric based on angular standard deviation and that based on PSI may not always agree. In our work, we have chosen the same metric (PSI) for the two ways (based on HT and WT) of calculating the phase difference between two signals. It is possible that the discrepancies found between the significance of the two PSI metrics in our study are due to the different window functions, which in the wavelet approach is the shape of either the real or imaginary part of the wavelet and for the Hilbert transform is the shape of the impulse response function of the band pass filter.
Similar to Bruns' work, 30 we want express our ideas about a comparison among STFT, WT, and HT for studying coherence and synchronization in CHS. There are basically two reasons to prefer STFT or WT instead of HT: (1) STFT and WT yield two time-frequency resolved parameters (coherence and phase difference between two signals), whereas HT yields only one (the phase difference); (2) STFT and WT imply less computational effort than HT, at least if we want to determine the instantaneous phase at a number of distinct frequencies, since HT requires previous application of narrow band filtering at each frequency which can be a computationally heavy process. We note though that these drawbacks of HT might be reduced if one applies the empirical mode decomposition 67 (EMD), previous to the application of the HT to the data in order to find the so-called Hilbert spectrum (where both the amplitudes and the instantaneous frequencies of the intrinsic mode functions (IMF) are expressed as a function of time). However, in principle, one would want to decompose the original signals in IMFs which are associated with clearly identifiable physical processes which still seem to be an issue for EMD. 68 About the comparison between WT and STFT, research in our group applied to CHS has found that these two methods are similar and could both be used equally well for data analysis. In STFT, given a certain window size and type, it is more straightforward to assign a time and frequency resolution. On the contrary, as we have repeatedly discussed in this work, WT deals with scales and the "translation" of scales into frequencies should be done always with caution. In particular, one should not confuse the characteristic frequencies as the only frequencies of the scaled wavelets and their spacing (Δf) as the frequency resolution associated with WT. For example, in our study, for a typical characteristic frequencies list output by the MATLAB function "wcoherence", the spacing Δf becomes smaller than the inverse of the total experiment duration, which obviously cannot be the true resolution.
In this study, we have focussed on the covariation between two independent data set like abp and T. However, in NIRS, we are often interested to study the covariation between oxy-and deoxyhemoglobin concentrations. These hemoglobin species are often derived by using mBLL, therefore, by using a linear combination of the same normalized changes of intensities at two wavelengths. The linear combination introduces a spurious correlation between the two hemoglobin species. For example, when we use λ 1 = 690 nm and λ 2 = 830 nm, we can show that even when two sequences of random numbers are chosen as normalized intensities at the two wavelengths, the phase difference between oxy-and deoxyhemoglobin is weakly peaked around 180°. Given a fixed protocol, the study of coherence and PSI in this case shows that higher thresholds are found for oxy-and deoxyhemoglobin than for abp and T.

Conclusions
In this work, we have investigated the significance of covariation between ABP oscillations (abp) and cerebral concentrations of total hemoglobin (T) in human subjects during a protocol involving cyclic inflation and release of pneumatic thigh cuffs. We focussed on the question of covariation, in the sense of both coherence and phase synchronization, between abp and T, which is directly relevant in the new technique of CHS. Even though the data processing method was tailored to our protocol and CHS, the scope of our work is broader because it aims at defining suitable metrics of significant covariation for any pair of signals.
The main goal of this work was to determine criteria to identify coherent/synchronous hemodynamic oscillations that are suitable for CHS. On the basis of our results, we propose to use independent Gaussian random numbers (GRN surrogate data) to generate a nullhypothesis statistical distribution of coherence, whose 95th percentile is taken as a critical (or threshold) value for significance. Here, we have focussed on coherence between abp and T, but the methods presented are applicable to the characterization of coherence between any physiological and cerebral hemodynamics measures. However, we have pointed out that if one were to consider the coherence between the concentrations of oxy-hemoglobin (O) and deoxyhemoglobin (D) as measured with NIRS, care must be taken to properly take into account the intrinsic correlation introduced by the linear relationship between optical intensities and O and D. In general, the existence of an intrinsic correlation between two signals results in a greater critical value of significance for their coherence.
We also propose to use a wavelet-transform-based (or STFT-based) measure of coherence, which allows for a time-resolved assessment of significant coherence during the course of a CHS measurement. This latter result is particularly important because it enables the identification of coherent hemodynamics during periods of induced physiological perturbations as well as during rest periods, in which only spontaneous hemodynamic oscillations are present. The ability to apply CHS to resting subjects, without a need to induce controlled perturbations in ABP, may result in a broader applicability of CHS for monitoring microvascular reactivity, cerebral perfusion, and vascular brain health.  Phase difference between T and abp at the cuff frequencies. On the left panel the phase differences were calculated as averages across the main lobe of the cuff frequencies with the analytic signal (HT), wavelet cross spectrum (WT), and cross power spectral density (FT) methods. For the wavelet (WT) and analytic signal (HT) methods, the averages were also carried out in the time intervals when the cuff was oscillating at each specific frequency. On the right panel the time-frequency averages of the phase differences were calculated with analytic signal and wavelet methods during baseline. The vertical lines define the different cuff frequencies. The error bars are the standard deviations across the frequencies of the main lobe of the cuff frequencies (FT) or across the temporal ranges of cuff oscillations (WT and HT).  The threshold values of coherence and PSI are shown for five different methods to generate surrogate data together with the values calculated for an experimental data set in vivo (subject 1). The x axis is defined as in Fig. 2. Specifically, in the four panels are plotted: linear coherence (a), wavelet coherence (b), PSI calculated with analytic signal (c), and PSI calculated with wavelet cross spectrum (d). The methods for surrogate data are: GRN (diamonds), IID1 (squares), IID2 (crosses), FT1 (circles), and AAFT1 (asterisks). The experimental data are represented by triangles.

Author Manuscript
Author Manuscript