The Study of EEG Dynamics During Anesthesia with Cross-Recurrence Rate

Introduction Intra-operative awareness has a reported incidence of approximately 0.1-0.8%. Anestheticinduced changes in the patient electrical brain activity (EEG) can be monitored in order to avoid intra-operative awareness cases. In this work we investigate the effect of anesthetic administration on the observed patterns of nonlinear dynamic coupling between different brain areas.


Introduction
General anesthesia is a chemically-induced reversible state of unconsciousness and depression of reflexes to afferent stimuli. Modern anesthesia involves the administration of different drugs to achieve the desired components of unconsciousness, amnesia, analgesia and immobility. Awareness during general anesthesia is considered a rare event; however, the psychological consequences can be severe for those who experience it. The incidence of awareness ranges from 0.1-0.8% [1] and is affected by a number of factors, such as patient characteristics and the type of surgery [2]. Devices that monitor the depth of anesthesia are now commercially available; such devices could provide a valuable means of identifying awareness during surgery, particularly since the patient himself cannot communicate this to the anesthetist due to the induced immobility. These devices function by monitoring the frontal electrical brain activity (EEG) from 2 electrodes placed on the patient's forehead. Particular characteristics are extracted from this recorded activity, which are then converted into a number corresponding to the level of hypnosis (0-100; no activity -fully awake, respectively). Available monitors are based largely on changes in the spectral content of the EEG, which are non-unique to anesthesia [3,4].
It is clear that an important issue in monitoring awareness during surgery is the particular set of methods used to extract features from the EEG activity. This is evident from (i) the large variations in classification accuracy or prediction probabilities used to assess the performance of different algorithms (see, e.g., table 5 in [5]), and (ii) the disagreement between different commercially available depth of anesthesia monitors in assessing patient state of hypnosis [6][7][8]. In addition to commercially available monitors, other algorithms utilized so far for monitoring awareness include spectrum-based methods [9][10][11][12], entropy-based methods [13][14][15], and methods from non-linear dynamics and complexity [13,14,16]. Recently, the use of recurrence methods has been introduced to study how the EEG activity is affected by the administration of anesthetic agents during general anesthesia and results indicate that such methods are highly promising for monitoring anesthetic depth. Recurrence methods are advanced nonlinear techniques that study the recurrence of the phase space trajectories of a dynamical system [17]. Two of the most popular recurrence methods are order recurrence plots, a tool for visualizing the dynamics of the phase space trajectories of the system [17], and order recurrence rate (ORR), an adaptation of recurrence plots that improves robustness against amplitude distortions of the time series under consideration [18]. Recurrence methods are advantageous over other non-linear methods used in anesthesia studies: they are suitable for short time-series [19,20], they are robust to the presence of artifacts [18,21], they are invariant to arbitrary transformations of the amplitude [18], and no assumptions on stationarity or linearity are necessary [22]. In a study by Li et al. the application of recurrence methods revealed an increase in the determinism of the frontal EEG with increasing concentrations of sevoflurane [23]. In other studies it was also found that recurrence quantification analysis displayed an increased ability, as measured with prediction probability, to separate consciousness from unconsciousness under different anesthetic regimes utilizing information from single-channel or twochannel frontal EEG [21,24].
In this work we investigate whether the application of cross-recurrence rate [21,25], a bivariate extension of recurrence rate, yields additional information to the study of nonlinear coupling relationships between different brain areas during anesthesia. This study was encouraged by previous investigations, where the CRR was applied to EEG data of 10 patients during recovery of consciousness at the end of surgery [26]. Here we utilize whole-head EEG data from the entire duration of surgery and without imposing restrictions on the anesthetic regime in order to characterize the evolution of the EEG complexity between different brain areas over the entire scalp and under actual surgical conditions. The main difference with other recurrence-based methods is that we are characterizing the evolution of the bivariate dynamical relationships between two brain areas, as opposed to characterizing the evolution of activity in a single brain area.

Dataset
The dataset used in this study was collected from 22 patients of age 45.7±21.3 (1 female patient) undergoing routine general surgery at Nicosia General Hospital, Cyprus. The study was approved by the National Bioethics Committee of Cyprus and patients gave written informed consent for their participation. A detailed description of the dataset is provided in previous work (see, e.g. [5]). In summary, anesthesia was induced by the on duty anesthetist with a propofol bolus and was maintained either (i) with an intravenous administration of propofol (17 patients) at concentrations ranging between 20-50 ml h -1 (200-500 mg h -1 ) depending on patient characteristics and surgery requirements, or (ii) with a continuous inhalational administration of sevoflurane (3 patients), desflurane (1 patient), or isoflurane (1 patient). In most patients remifentanil hydrochloride (Ultiva®; 2 mg, dissolved in 40 ml) was also administered intravenously throughout surgery at a rate ranging between 2-15 ml h -1 (0.1-0.75 mg h -1 ). Patient lungs were ventilated with an air-oxygen or air-oxygen-N 2 O mixture. Some patients were also given boluses of neuromuscular blocking agents (tracrium, vecuronium or sisatracurium). The particular drugs administered (propofol and potent inhaled anesthetics) follow the basic anesthetic-induced EEG pattern, i.e. changes in both frequency and amplitude (amplitude increase and frequency slowing), with the possibility of burst suppression in high doses.
EEG data were collected with the TruScan32 system (Deymed Diagnostic) and 19 electrodes were placed according to the International 10/20 system: Fp1, Fp2, F7, F3, Fz, F4, F8, T3, C3, Cz, C4, T4, T5, P3, Pz, P4, T6, O1 and O2, with an FCz electrode reference. The sampling rate was 256 Hz and the data were analog lowpass filtered with cutoff at 100 Hz; thus, the sampling frequency of 256 Hz was more than adequate to ensure fulfillment of the Nyquist frequency consideration. No additional filtering was performed during or after data collection. Data recording usually commenced while patients were still awake prior to anesthetic induction, continued throughout the entire surgery and until patients regained consciousness at the end of surgery. Manual markers were inserted in the EEG records indicating important events, such as anesthetic induction, recovery of consciousness and start/end of anesthetic administration. These markers are important in subsequent identification of EEG windows corresponding to wakefulness and anesthesia. The point at which the patient began responding to verbal commands or tactile stimuli by the anesthetist at the end of surgery was defined as recovery of consciousness. This was not influenced by neuromuscular blocking agent administration as an antagonist was administered to cancel the neuromuscular blocking effects, if necessary, once anesthetic administration was discontinued. Depending on patient characteristics loss of consciousness occurred approximately 10-30 seconds after anesthetic induction; this was assessed by the anesthetist on duty either as a loss of eyelash reflex, lack of response to light tapping of the patient's forehead just above the nose between the eyes, or a lack of patient response to a verbal command.

Estimation of cross-recurrence rate
The first step in estimating the CRR is to represent the dynamics of the system as a symbolic sequence and to encode them as order patterns. The main advantage of a conversion to a symbolic sequence is that the order patterns remain invariant under strict monotonic amplitude transformations [18]. Considering a onedimensional EEG time series at a given electrode, , a phase-space representation of the EEG time series can be constructed by time delay embedding, with an embedding dimension m and embedding delay τ ed [27]: The process of time delay embedding with embedding dimension m effectively provides a phase-space representation of the system, without changing the underlying system dynamics. Each m-dimensional segment is ordered in terms of rank. For example, for m=3, an embedded segment with values [3, 6, 2] would correspond to the rank vector [1,2,0]. This sequence of ranks is the order pattern, , of the embedded segment. In general, there are m! possible ways (order patterns) , in which a sequence of m values can be ordered. The recurrence matrix, , identifies a match between the order patterns of two embedded EEG signals, and , separated by a delay, \tau_{ord} (order recurrence delay): where The cross recurrence rate, , is a bivariate measure representing the frequency of order pattern recurrences. CRR is estimated as: The factor is a normalizing factor such that . The CRR allows us to study the recurrent dynamics between different systems by estimating how many times a particular state occurs simultaneously in and [28]. This is very useful for the analysis of delayed interacting systems or systems with feedback [20].
The choice of the parameters m and is important and methods for estimating them should be used [29]. We chose here (4 ms) and (see Discussion for more details).

Cross-recurrence rate flow
The non-symmetric nature of the CRR measure allows us to identify the direction of the temporal evolution of the non-linear bivariate relationships. In this context, we interpreted as EEG signal leading EEG signal (and vice versa). Thus, we can investigate whether a given electrode, , is in general a leader (eq. 4) or a follower (eq. 5), by comparing the CRR values obtained when the specific electrode leads or lags other electrodes respectively.
(4) (5) where N e is the number of available electrodes for a particular patient. These two quantities allow us to characterize whether: (i) the activity from electrode leads or lags activity generated at electrode ; and (ii) whether a particular electrode is in general a 'leader' or a 'follower' . This information can also be used to obtain an overall picture of which brain areas (defined by the corresponding electrodes) act as 'leaders' or 'followers' during wakefulness and anesthesia. This is obtained by counting, for each patient, the number of windows, N, that displayed either higher lead or lag activity during wakefulness and anesthesia at each electrode, . If , then electrode is identified as a general 'leader'; otherwise, if , then electrode is identified as a 'follower'.

Data analysis
The EEG data from each patient were analyzed using a moving and non-overlapping window. Windows corresponding to wakefulness and anesthesia were defined based on the markers inserted in the continuous EEG records during data recording: anesthesia is the activity contained within the markers for 'induction' and 'recovery of consciousness'. Given that recovery of consciousness is an abrupt state change

Statistical analysis
Two tests for statistical significance were applied. Firstly, the method of phase randomized surrogate data was used to assess the statistical significance of each estimated CRR value [30]. The surrogate data has the same spectral properties as the original data, but (nonlinear) phase relationships are destroyed. This is important as it is clear from the CRR formulation that the CRR can also be related to the phase of the time series. The number of surrogate signals needed for a given significance level is a function of the probability of false rejection, , and can be estimated as [31]. We now have datasets a total, including both the surrogates and the original dataset. The probability that the original data has the largest CRR of this dataset merely by coincidence is exactly . Hence, if the estimated CRR exceeds the maximum surrogate CRR, we can accept the hypothesis that the estimated CRR is significant at a significance level . In this work 19 surrogate signal pairs (i.e. ) were generated for each electrode signal.
Secondly, the effects of condition (wakefulness vs anesthesia), recording location and type of anesthesia (intravenous vs inhalational) on the mean CRR values were assessed. To test for these effects over all patients a 3-way analysis of variance (ANOVA F-test, ) was performed (command ' anovan' in Matlab). The mean CRR values used in the ANOVA tests were obtained from all CRR values, without excluding those that were found to be non-significant from the surrogate method. Figure 1 shows the topography of the patient-wise average lead CRR flow during wakefulness and anesthesia (the lag CRR flow is simply the reverse). The scalp maps were obtained using the function 'topoplot' (part of the EEGLAB toolbox [32]). To obtain a smooth topography 'topoplot' uses interpolation for areas in between the available electrodes. The topographies were obtained by considering only the significant CRR values, i.e. the CRR values that exceeded the maximum surrogate CRR value at each time window. To obtain smooth topographies, Figure 1 shows the overall picture of which brain areas (defined by the corresponding electrodes) act as 'leaders' or 'followers' during wakefulness and anesthesia (as described in the section 'CRR Flow'). Thus, Figure 1 provides a broad (patient-wise average) picture of the general role of each electrode as a 'leader' or 'follower'.  4) and (5)) for patients S6, S2 and S5 respectively. As described in the previous section, the depicted CRR values at each electrode represent the sum (hence, values can be > 1) of incoming or outgoing CRR between the particular electrode and all other electrodes. Sharp increases in both lead and lag CRR correlate with wakefulnessanesthesia transitions and not with discontinuation of anesthetic administration. For visualization purposes, the estimated CRR values were smoothed with a moving average filter . Tables 1, 2 show the median CRR values (estimated over all electrodes) for each patient individually, for outgoing and incoming CRR respectively. The median is presented instead of the mean, as the median does not take into account outlier values in its estimation. The increase of CRR values during anesthesia compared to wakefulness is observed for all patients studied, except for patient S1 (we have not yet been able to identify any particular characteristics of this patient that can explain this behaviour).

FIGURE 1: Patient-wise mean tendency of CRR 'outflow' during (a) wakefulness and (b) anesthesia ('inflow' activity is simply the reverse).
Higher CRR values indicate larger similarities. During wakefulness the 'outflow' of activity is more spatially focused over centro-posterior areas, suggesting that flow of information initiates in posterior areas. Administration of anesthesia causes a more spatially widespread 'outflow' of activity, with the exception of central areas, which become the main recipient of this widespread activity. Values are interpolated over the entire scalp to obtain a smooth topography.

FIGURE 2: Lead and lag CRR activity at electrode Cz for patient S6, estimated using non-overlapping EEG windows of 2 seconds duration.
Anesthesia was induced and maintained with propofol. Vertical lines indicate anesthetic induction (LOC), recovery of consciousness (ROC) and termination of anesthetic maintenance (AN OFF). The use of diathermy equipment causes large outlier values.     show the corresponding scalp topographies of the average lead and lag CRR values (top and bottom rows respectively), averaged over all EEG windows that correspond to wakefulness or anesthesia (left and right columns respectively), for patients S1, S2 and S4 respectively. In contrast to Figure 1, which represents the overall patient-wise tendency of particular brain areas (electrodes) to act as leaders or followers, Figures 5-7 show individual patient lead and lag CRR values (top and bottom rows respectively) at each EEG electrode, averaged over all EEG windows that correspond to either wakefulness or anesthesia (left and right columns respectively) for the particular patient. Due to the inter-subject variability of the CRR values for individual patients different scales are used in plotting each figure. However, the patterns observed are the same, regardless of the actual CRR values. A video showing the estimated CRR topographies for the entire EEG record of patient S1 can be provided per request. In addition, the CRR patterns are not affected by the window size (2 or 4 seconds). The 3-way ANOVA test (F-test, ) showed significant effects on the CRR values. It was found that condition, recording location and anesthetic type had a significant effect on the CRR values. The mean CRR during anesthesia was significantly different from the mean CRR during wakefulness (lag CRR: F=384.8, p=0; lead CRR: F=374.4, p=0). In addition, a significant effect of recording location on the mean CRR values during wakefulness (lag CRR: F=51.8, p=0; lead CRR: F=50.0, p=0) and a significant effect of anesthetic type on the mean CRR values during anesthesia (lag CRR: F=5.2, p=0; lead CRR: F=5.11, p=0) were identified. The particular effects of these factors can be summarized as follows: (a) maintenance with propofol resulted in lower mean CRR values during wakefulness compared to maintenance with inhalational agents (mean lag CRR of 0.22 vs 0.29 respectively; values were identical for mean lead CRR); and (b) more posterior locations displayed larger mean CRR values than those at fronto-central locations (Figure 8 -showing mean lag CRR values; these were identical for mean lead CRR). The general increase in CRR values during anesthesia is more prominent at posterior locations.

Discussion
The most significant observation is that administration of anesthesia induces a spatially widespread increase in the CRR. This is independent of individual patient characteristics and of the particular anesthetic protocol followed. Even though the anesthetic protocol does not affect the observed CRR patterns, it does have a significant effect on the CRR values during wakefulness. The lower mean CRR values during wakefulness after propofol maintenance compared to maintenance with isoflurane and sevoflurane suggest that inhalational agents may have a larger suppressive effect on the observed EEG activity. In addition, frontal and centro-frontal areas display significantly smaller CRR values compared to other brain areas; however, the amount of lead and lag activity from these areas is at a similar level. During wakefulness we can also identify larger lead activity from posterior areas compared to all other areas (Figure 1), which represents the tendency of these areas to act as initiators of the observed activity.
Administration of anesthesia induces more spatially widespread nonlinear coupling relationships by entraining more areas into a synchronized activity, resulting in a more spatially widespread activity flow and the disruption of the 'leader' relationship identified during wakefulness. This is with the exception of fronto-central areas (F3, Fz, F4), which remain as 'followers'. This occurs consistently across all patients studied, although a considerable inter-subject variability in the actual CRR values was found. Considering the patient-wise average CRR for wakefulness and anesthesia ( Figure 1) some important implications can be suggested. During wakefulness there is a bigger tendency for centro-posterior areas to display larger lead activity, thus these areas can be considered as 'leaders'.
These observations can be related to other work where it is suggested that EEG is more irregular during wakefulness and becomes more predictable during anesthetic-induced unconsciousness, thus leading to a decrease of the EEG complexity. This is reflected in the low CRR values during wakefulness implying that there is no significant coupling of the phase-space trajectories. Administration of anesthesia induces strong bidirectional coupling of the phase-space dynamics, which is identified as a significant increase in the CRR values. In addition, parietal and frontal areas display higher lead and lag CRR values compared to anterior central areas. This weaker relationship of the EEG dynamics could reflect some form of a decoupling between anterior central and other areas. As seen from Figures 2-7, as well as from the supplementary figures SF1-SF3 and the supplementary video SV1, the observed patterns are similar regardless of individual patient characteristics and of the particular anesthetic regime followed. Evidence from recent imaging studies with isoflurane and halothane [33], sevoflurane [34], midazolam [35] and propofol [36] indicate a functional disconnection of the cortex from outside sensory experience. The disruption of sensory information flow to the cortex implies that no external input can reach the cortex resulting in a hypersynchrony of neuronal activity as the unstimulated brain areas begin to oscillate together. Consequently, different brain regions of the cortex oscillate together resembling coupled oscillators. The effects of this are reflected as an increase in the CRR values as observed in the broader context of the non-invasive EEG.
It is very likely that the observed changes are not a mere reflection of the concentration of the anesthetic agent in the body. This is supported by considering that the administration of anesthetics was switched off some minutes prior to recovery of consciousness. Should the estimated CRR values have been a reflection of the metabolic decrease of the anesthetic agent, a gradual decrease in the CRR values would have been observed rather than the sharp decrease observed here. In addition, it is known that surgical noxious stimuli tends to lighten the level of hypnosis [37]. Even though recurrence methods are robust to artifacts, it is likely that some outlier CRR values are also present given that no artifact removal has been performed. However, despite these outlier values the observed CRR patterns remain stable from the onset of loss of consciousness to its recovery at the end of surgery and are neither affected by, nor are a direct result of, the surgery itself.
The robustness and stability of the CRR patterns is important for the clinical utility of the proposed methodology. Despite the inter-subject variability in the actual CRR values, the observed patterns are similar for all patients. In a clinical setting, a patient-specific baseline CRR during wakefulness can be obtained prior to administration of the anesthetic bolus. Bolus administration will cause the CRR values to increase above this baseline level, where they will remain stable provided the patient is anesthetized. This is observed for 21 of the 22 patients studied (Tables 1, 2). The anesthetist can then monitor the CRR values and ensure that these do not return to the patient baseline level during surgery. Therefore, the known consciousness/unconsciousness states of the patient prior and immediately after anesthetic bolus administration allow system calibration for each patient. Thus, the proposed method is generalized, as the CRR patterns are similar for all patients, while at the same time individual variability in the CRR values is taken into consideration. This is an important advantage over existing monitoring methods that do not take into account inter-subject variability and employ a universal scale from 0 to 100. In addition, the robustness of the CRR patterns over all recording locations studied allows the use of a smaller number of electrodes without affecting performance.
The choice of the particular values of and was based on previous findings. Jordan et al. state that the ability to separate the two states decreases for order recurrence delays [21]. We also found a similar effect from our previous investigations during recovery of consciousness using the CRR [26], where we concluded that: (i) a value of provided the best discriminating ability between wakefulness and anesthesia; and (ii) the choice of m affects the estimated CRR values, but not the observed CRR patterns. Thus, based on our previous investigations, and given that large values of m constitute the estimations very computationally complex, we chose here and (a value of is also recommended in [18]).
As a final remark, the CRR patterns were estimated using a reference electrode montage (reference location at FCz). Even though the particular choice of reference may have some effects on the estimated activity, an "ideal" reference should effectively be placed far away from the underlying source generator. However, no such distant point is generally available in EEG recordings; thus, there is essentially no difference between recording and reference electrodes. Other popular montages include (i) the linked-ears reference, which is commonly used despite the minimal theoretical basis supporting this approach, (ii) bipolar recordings, which is useful for improving the spatial resolution of the EEG, and (iii) average reference, which is most likely a poor choice for use with the standard 10/20 system (for an extensive discussion of the various montages and their advantages/disadvantages see [38]).

Conclusions
We investigated the nonlinear bidirectional interactions during wakefulness and anesthesia and identified a spatially widespread increase of interactions during anesthesia that is likely a result of an increase in synchronization. The study has important implications when considering the feasibility of cross-recurrence methods in identifying between anesthesia and wakefulness. In a previous study CRR displayed an increased ability (98%) to separate wakefulness and anesthesia [26]. Taking this into consideration and also considering the speed of estimation and the statistically significant changes observed in the CRR the future application of CRR in a device for monitoring anesthesia during surgery is encouraged.

Additional Information Disclosures
Human subjects: All authors have confirmed that this study did not involve human participants or tissue. Animal subjects: All authors have confirmed that this study did not involve animal subjects or tissue.

Conflicts of interest:
In compliance with the ICMJE uniform disclosure form, all authors declare the following: Payment/services info: All authors have declared that no financial support was received from any organization for the submitted work. Financial relationships: All authors have declared that they have no financial relationships at present or within the previous three years with any organizations that might have an interest in the submitted work. Other relationships: All authors have declared that there are no other relationships or activities that could appear to have influenced the submitted work.