Windowed Multitaper Correlation Analysis of Multimodal Brain Monitoring Parameters

Although multimodal monitoring sets the standard in daily practice of neurocritical care, problem-oriented analysis tools to interpret the huge amount of data are lacking. Recently a mathematical model was presented that simulates the cerebral perfusion and oxygen supply in case of a severe head trauma, predicting the appearance of distinct correlations between arterial blood pressure and intracranial pressure. In this study we present a set of mathematical tools that reliably detect the predicted correlations in data recorded at a neurocritical care unit. The time resolved correlations will be identified by a windowing technique combined with Fourier-based coherence calculations. The phasing of the data is detected by means of Hilbert phase difference within the above mentioned windows. A statistical testing method is introduced that allows tuning the parameters of the windowing method in such a way that a predefined accuracy is reached. With this method the data of fifteen patients were examined in which we found the predicted correlation in each patient. Additionally it could be shown that the occurrence of a distinct correlation parameter, called scp, represents a predictive value of high quality for the patients outcome.


Introduction
The outcome of patients suffering from catastrophic neurological events such as traumatic brain injury (TBI) or subarachnoid hemorrhage (SAH) is influenced by two types of pathophysiological mechanisms: the primary injury sustained at the moment of impact and the secondary injury caused by a cascade of ischemic, ionic, and neurochemical insults which evolve over time after the primary insult [1]. Since the primary injuries are merely irreversible, the main focus of neurocritical care is to detect and prevent secondary brain injury [2]. Consequently, an array of different monitoring techniques for the assessment of intracranial pressure (ICP), oxygenation status, and metabolism has been developed to unmask occult local and systemic alterations of the brain [3][4][5][6]. However, the resulting high volume of multimodal datasets frequently exceeds the ability of the neurointensivist to process and integrate these data into an adequate treatment algorithm [7]. There is an urgent need for the development of advanced biomathematical tools to analyze the complex and interwoven data sets in real time guiding a goal directed therapy of these patients [8]. Our group has recently developed a mathematical model which allows simulating brain swelling and loss of cerebral autoregulation [9]. In case of a severe brain swelling, this model predicts a long-term positive correlation between arterial blood pressure (ABP) and ICP if the cerebral autoregulation is experimentally switched off, in contrast to a negative correlation for an intact autoregulation [10]. The main goal of our study was to develop a set of mathematical tools that reliably identify such positive and negative correlations in bimodal monitoring data. The correlations are described by three classification numbers, called selected correlation (sc), selected correlation positive (scp), and selected correlation negative (scn). The sc index is deduced from windowed multitaper coherence calculations and multitaper power spectra [11]. The phasing of the data windows is detected by the mean value of the windowed Hilbert phase difference [12] leading to the indices scp and scn. For all indices a statistical testing method is presented to determine the significance of the particular results. In a second step we calculated the percentage allocations of positive and negative correlations and related the results to the outcome of the analyzed patients. We hypothesized that patients showing predominantly positive correlations will have a poor prognosis. In addition we attempt to validate the predicted episodes of diminished cerebral compliance by analysis of CT scans.

Patient Population.
The study was conducted in accordance with the ethical guidelines of the University Regensburg Institutional Review Board. Informed consent was obtained from the patient's relatives, all study results were stored and analyzed in an anonymized fashion. We prospectively investigated a cohort of 15 adult patients (5 female, 10 male) who were treated at the neurosurgical intensive care unit for traumatic brain injury (TBI) or subarachnoid hemorrhage in 10 and 5 cases, respectively. The mean age was 43.7 years (range: 20.2-72.4), the median Glasgow Coma Scale GCS at the time of admission was 7 (range: 3-15). The mean duration of multimodal monitoring was 142 hours (range: 20.3-349.4). The detailed patient characteristics are summarized in Table 1. All patients were sedated and mechanically ventilated during the observation period and received an intra-arterial catheter for the continuous measurement of arterial blood pressure as part of the standard treatment procedure in our institution. ICP monitoring was performed continuously using either an external ventricular drain equipped with an electronic pressure device (EVD) or a parenchymal ICP probe (both from Raumedic, Helmbrechts, Germany). The ABP and ICP data were acquired continuously using a data logger (Daq USB 6210, National Instruments, Munich, Germany) with a sample frequency of 1000 Hz. For the correlation analysis, the data were resampled to 0.2 Hz (one data point every five seconds) to reduce noise effects and to smooth out fast oscillations or spikes. Computed tomography imaging was conducted when clinically indicated. Followup was completed up to July 2013 by reviewing outpatient records and contacting the patients family member or the patient's primary physician. The mean followup time was 77.9 months, no patient was lost for followup. The neurological outcome was measured by the Glasgow Outcome Scale (GOS) [13], and revealed a median score of 3 (range: 1-5) at the final followup exam.

Correlation Index Calculations.
To identify the above mentioned positive and negative correlations between ABP and ICP in monitoring data from the ICU, we used a Fourier based method, as described in detail previously [14]. In brief, we adopt coherence and power spectra calculations with the multitaper method (MTM) to determine the coherence between segments of two time series that were synchronously recorded with a sampling rate of 0.2 Hz. Additionally, we calculate the mean Hilbert phase difference between these segments to determine the phasing of the correlation. The MTM [11,15], based on fast Fourier transforms (FFT), is a sophisticated way to determine power spectra and coherence of data using a set of so-called Slepian tapers as window functions to handle noisy and short data. Additionally the MTM software comes with built-in statistical tests proving the significance of each frequency. The principle approach of the above mentioned windowed MTM calculations can be described as follows: From the isochronous time series , select windows , of fixed size with = 2 , ∈ N and potentially different starting points and . For this pair of windows we calculate the mtm-spectra mtms( ) and the mtm-coherence between the windows (mtmc( )): ( With the above mentioned setting for the window size the MTM-spectra and MTM-coherence contain exactly /2 frequencies . Using the built-in significance test of the MTM each single frequency is tested for significance. Thus, we can define the pointwise selected correlation (PSC) assuming a significance level for the built-in significance test of 99%: PSC , , := (psc , 1 ( , ) , . . . , psc , /2 ( , )) with: The PSC represents a list of length /2 with elements of value 1 if the correspondent frequency is considered significant in the spectra and the coherence; otherwise the value of the element is 0. The significance of in both spectra guarantees that only frequencies are considered that essentially contribute to the original signals, whereas the significance in the coherence assures that a specific implies correlation between the input signals. Repeating the PSC calculations for pairs of isochronous windows leads to the mean pointwise selected correlation (MPSC): The elements of the MPSC list represent the percentage of being significant in both spectra as well as in the coherence calculation for each single frequency . Therefore, using MPSC we are able to determine frequency intervals reflecting relevant correlations within a whole data set. If distinctive correlations occur in the data set, the related frequencies lead to an increase in the index MPSC ( ). Calculating the MPSC for a complete data set reveals in which frequency range correlations occur on average. After having identified such a frequency interval = ( , . . . , ) we want to determine periods in the data set where strong correlation with respect to occurs. To do so, we first estimate the degree of correlation of a distinct pair of windows with respect to the interval by calculating the sum of all elements psc of PSC Table 1: Characteristics of patients included in the study who were treated for either SAH (subarachnoid hemorrhage) or TBI (traumatic brain injury). The column GCS describes the level of consciousness at admission according to the Glasgow Coma Scale. Measuring time reports the duration of multimodal monitoring interval per patient. Over all correlation describes the percentage of significant correlations within the observation period. The last column indicates the outcome determined by the Glasgow Outcome Scale at the final follow up.

ID
A pair of windows is called selected correlated if sc > lsc for any chosen threshold lsc. A statistical test, to be described later on, is used to determine the significance of the lsc value. The sc value therefore serves as a measure for the degree of correlation of a pair of data windows with respect to a specific frequency range where the specified correlation occurs. Consequently, to gather time resolved information about the correlation between ABP and ICP time series in a specific frequency range, we determine the index sc , , ( , ) for isochronous windows while shifting the starting point along the time axis (Figure 1).

Hilbert Phase Differences.
Having identified a pair of windows exhibiting a sufficient high correlation index sc, we calculate, in a next step, the Hilbert phase difference [12] of the corresponding time series segments to detect the phasing of the signals. For two time series ( ), ( ) the Hilbert phase difference is defined by the components of the corresponding analytical signals, analytic = ( ) + * ̃( ) = ( ) * * ( ) with: Each time the sc value is higher than limit lsc a profound correlation between ABP and ICP is detected. If additionally the mhpd value surpasses either lmhpd pos or lmhpd neg a distinctive negative correlation (depicted as dots) or a distinctive positive correlation (depicted as crosses) is detected.
where P.V. stands for the Cauchy principal value of the of the integral considered. Prior to this calculation both segments are normalized to zero mean and unit standard deviation, just as it is done for the MTM-calculations. From the resulting hpd( ) we calculate the mean value of the Hilbert phase difference as a measure for the phasing of the two input windows.
A correlation is called negative if the mean Hilbert phase difference (mhpd) is higher a predefined limit lmhpd neg and correspondingly is called positive if lower a predefined lmhpd pos . Finally, a statistical test is used to determine the significance of lmhpd neg and lmhpd pos .

Statistical Test.
The statistical test, presented in the following, is a specific type of perturbation test ( [16]) and based on the model prediction that correlations between ABP and ICP occur isochronous. Consequently, two segments should not be correlated if their starting points are quite apart from each other and the corresponding sc value should be zero or at least very low. Assuming that a sc value is meaningful if it is higher a predefined specific threshold lsc, we can count how often this separated and therefore uncorrelated windows produce sc values higher than lsc. From this obviously wrong hit we deduce the error rate and finally the significance of the sc value with respect to the predefined threshold. In a first step we will determine an appropriate separation length, or time shift, to guarantee the uncoupling of the time series with respect to the chosen segment length . If at least one of the signals exhibits a distinctive autocorrelation the time shift between the segments has to be chosen big enough to minimize this source of irritation. First we construct test data sets ABP all and ICP all by simply joining the appropriate data sets of all patients involved. Within such a test data set we randomly determine a starting point for the first window of length . The starting point of the second window, then, results from the first starting point plus a specified time shift . Subsequently we can calculate the sc value for this distinct pair of windows. To cover the characteristics of the whole test data set this procedure is repeated times. Afterwards, the time shift is increased, and the whole procedure is repeated until ending in the calculation of the mean windowed autocorrelation (mwa): If the time shift if large enough to exclude autocorrelation artifacts, the subsequent mwa values should be small and stable. With the ability to exclude interference from autocorrelation effects, the above mentioned calculation of the error rate for sc with respect to a predefined threshold lsc could be done as follows: Using the ABP all and ICP all test data sets for the generation of the input windows, we repeatedly and randomly determine pairs of starting points ( , ) with time shift much bigger than autocorrelation effects (| − | > ac ). Again this is done to cover the characteristics of the whole test data set. For the associated pairs of nonisochronous windows, we calculate the selected correlation sc , , ( , ). As correlations should only occur for isochronous windows these values of sc reflect a threshold that should be passed by correlated windows. In other words, if we call a sc value meaningful if higher than a predefined threshold lsc, all sc values of such nonisochronous windows satisfying the lsc criterion are wrong hits. Accordingly, we deduce the error index ei , , ( , ), indicating whether the selected correlation sc , , ( , ) is higher than a predefined limit lsc and the error rate asc, that is, the rate of obviously wrong hits with respect to lsc:  The resulting error rate asc describes the percentage of obviously wrong hits in the above mentioned sense with respect to a predefined limit lsc. Thus asc represents the probability of a distinct value sc > lsc to be generated by accident and therefore describes the significance of sc with respect to lsc. The above described method can easily be modified to calculate the error rates of the mean Hilbert phase difference (mhpd) by substituting the lsc criterion in ei , , ( , ) through appropriate criteria lmhpd pos and lmhpd neg . If mhpd > lmhpd neg the correlation between the appropriate windows will be called negative; if mhpd < lmhpd pos the correlation will be called positive. Accordingly, the error rates will be named amhpd pos and amhpd neg , respectively, and will be used to determine the significance of the lmhpd pos and lmhpd neg . The above described tools will be used to determine for an arbitrary pair of input windows whether these windows are significantly correlated and, if so, whether the correlation is significantly positive or negative. Shifting the input windows synchronously along the time axis of the measured data produces a time resolved information about the different possible correlations (see Figure 1). Additionally we will calculate the percentage of correlated windows with respect to lsc and the percentage of distinct positive and negative correlations, with respect to lmhpd pos and lmhpd neg . The percentage of distinct positive correlated windows will be termed as selected correlation positive (scp) and correspondingly the percentage of distinct negative correlated windows will be termed as selected correlation negative (scn).

Results
For all subsequent calculation a window size of 2048 data points was used. In The mwa value decreases very rapidly for both signals. For offsets greater than 500 data points the mwa value remains stable at mwa < 0.011 (a). Error rate asc with respect to sc: choosing the criterion for a pair of windows to be correlated as lsc = 0.03 leads to a asc < 0.05 (b). Error rate amhpd neg with respect to the limit lmhpd neg for correlations to be called negative: a lmhpd neg of 130 deg for the mean Hilbert phase difference (mhpd) generates an amhpd neg < 0.05 (c). Error rate amhpd pos with respect to the limit lmhpd pos for correlations to be called positive: a lmhpd pos of 50 deg for the mean Hilbert phase difference generates an amhpd neg < 0.05 (d).
band Figure 2 only shows the MPSC value for frequencies lower than 0.015 Hz. For frequencies higher than 0.015 Hz MPSC values are mostly zero except of a view small regions with MPSC values lower than 0.00012. Therefore we choose a frequency interval of < 0.005 Hz for further analysis. The mean windowed autocorrelations (mwa) of ABP signal and ICP signal are depicted in Figure 3(a), where 5000 different pairs of input windows were used for the calculation of a single offset . If no time shift is applied (step size 0), the mwa value is highest but not one. This is due to the fact that in this case all exploited frequencies are sufficiently correlated but some of them prove to be not significant in the spectra. With growing step size the mean windowed autocorrelation decreases rapidly for both signals and at least after a step size of half the window size is reached and mwa is low and stable (mwa < 0.011). From this we can conclude that if the windows are simply nonoverlapping, the time shift is big enough to exclude autocorrelation effects in both signals. Therefore we used an offset of three times the window size ( = 3 * ) for the statistical testing. For the calculation of the different error rates 2 × 10 6 different pairs of input windows were used. The resulting error rates asc, amhpd pos , and amhpd neg are depicted in Figures 3(b)-3(d). From lsc > 0.03 an asc < 0.05 follows. If we choose lsc > 0.1, asc reduces to asc < 0.0015. Choosing the upper and lower limits for the mean Hilbert phase difference according to earlier work from our group [10] as 130 deg and 50 deg results in amhpd neg , amhpd pos < 0.05. Using 150 deg a 30 deg for the upper and lower limit reduces the error rates, amhpd pos , amhpd neg < 0.004 (see Figures 3(c) and 3(d)). To achieve a significance level higher than 95% for the detection of correlations we therefore used lsc = 0.03 and lmhpd neg = 130 deg, lmhpd pos = 50 deg for the detection of negative and positive correlations. Using the above defined thresholds, significant correlations between the parameters ABP and ICP were recorded in all patients; however, the time span of correlating parameters per observation period differed considerably between the individual patients (see Table 1). All patients, from which CT scans were obtained during positive selected correlation phases, indicating a loss of autoregulation combined with reduced intracranial compliance, showed dramatic changes indicating severe edema, hemorrhagic transformation, or brain ischemia (Figure 4).
In addition, the value of selected correlation positive (scp) during the observational period correlated highly significantly with the patients outcome as measured by the Glasgow Outcome Scale ( < 0.0001, see Figure 5). Interestingly, all patients showing scp values higher than 10% died during the course of the disease. Multiple regression analysis including age, diagnosis, scp, and the GCS score at admission revealed the selected correlation positive as independent factor for patients outcome ( < 0.001).

Discussion
Despite promising in vivo results in several animal models, investigating neuroprotective agents such as glutamate antagonists, free radical scavengers, or growth factors, clinical trials have failed to show any benefit for patients with TBI or SAH [17,18]. Since causal treatment is lacking so far, the main goal of neurointensive care is to provide an optimal physiological and cellular environment to prevent secondary injury and to facilitate endogenous regeneration of the injured brain [19]. To achieve this mission, a number of technically sophisticated monitoring devices measuring ICP [20], brain oxygenation [21], metabolism [22], cerebral blood flow [23], and EEG activity [24] have been developed recently. However, the interpretation of these multimodal data sets is until today based on experience and adequate intuition of the treating neurointensivist, who may be facing more than 200 variables per patient [6]. As shown previously, it is almost impossible to grasp the interrelation between more than two data sets without the support of biomathematical analysis tools [25]. For example, an elevated ICP may require the increase of the systemic blood pressure in a patient with functional autoregulation to provide adequate cerebral perfusion pressure (CPP), whereas it can be harmful to a patient with a loss of autoregulation [2]. It is therefore necessary to avoid dogmatic threshold values as a foundation for treatment algorithms [26] but to individually tailor therapeutic intervention for each patient [5,7]. This aspect becomes even more vital as recent studies have demonstrated significant adverse effects of treatment strategies targeting CPP or brain oxygenation if executed inadequately [27,28]. The main goal of our study was therefore to employ advanced biomathematical tools to identify patients at risk and to improve the pattern of care and prognostication in patients treated in the neurointensive care unit. With only 15 patients, the sample size of our cohort is rather small, which may be a potential limitation. Further studies are underway to prospectively validate our findings and their potential impact on treatment outcome.
In conclusion, the used methods are able to reliably detect correlations with different phasing as predicted by a mathematical model. By means of a statistical testing method the accuracy of the results could be tuned as needed for a reliable judgment of the patients status. Interestingly, the frequency of positive selected correlations highly was identified as an independent factor for the prognosis of our study patients.
Our results indicate that a real time application of this analysis would potentially allow for an early intervention improving the outcome of patients by preventing secondary injury in a more timely fashion [29]. Abbreviations amhpd neg : Error rate for limit of negative mhpd (percentage of false positives of mhpd with respect to lmhpd neg ) amhpd pos : Error rate for limit of positive mhpd (percentage of false positives of mhpd with respect to lmhpd pos ) asc: Error rate of selected correlation (percentage of false positives of sc with respect to lsc) ei: Error Index (ei is defined as 1 if the sc value is higher than lsc and 0 otherwise; for data windows of abp and icp with an offset, ei indicates false positives) FFT: Fast Fourier transform (method to calculate the discrete Fourier transformation of a dataset) hpd: Hilbert phase difference (difference of the Hilbert phase of two signals) lmhpd neg : Limit mean Hilbert phase difference negative (arbitrary limit for mhpd to indicate a negative phasing) lmhpd pos : Limit mean Hilbert phase difference positive (arbitrary limit for mhpd to indicate a positive phasing) lsc: Limit selected correlation (arbitrary limit for sc to indicate correlation) mhpd: Mean Hilbert phase difference (mean value of the Hilbert phase difference for a fixed data window) MPSC: Mean pointwise selected correlation (element by element mean of PSC vectors of a selection of data windows) MTM: Multitaper method (sophisticated fft based method to calculate power spectra and coherence from a data window) mtmc: Multitaper coherence (coherence calculated with multitaper method) mtms: Multitaper spectrum (power spectrum calculated with multitaper method) mwa: Mean windowed autocorrelation (mean value of sc for data windows of identical data (abp or icp) with a predefined offset but varying starting points) PSC: Pointwise selected correlation (a vector containing a 1 if the corresponding frequency is significant in the spectra and the coherence and 0 otherwise; the PSC is calculated for a fixed pair of corresponding data windows) sc: Selected correlation (sum over the elements of PSC normed by the length of PSC) scn: Selected correlation negative (percentage of sc with additionally a negative mhpd) scp: Selected correlation positive (percentage of sc with additionally a positive mhpd).