Time-variant coherence between heart rate variability and EEG activity in epileptic patients: an advanced coupling analysis between physiological networks

Time-variant coherence analysis between the heart rate variability (HRV) and the channel-related envelopes of adaptively selected EEG components was used as an indicator for the occurrence of (correlative) couplings between the central autonomic network (CAN) and the epileptic network before, during and after epileptic seizures. Two groups of patients were investigated, a group with left and a group with right hemispheric temporal lobe epilepsy. The individual EEG components were extracted by a signal-adaptive approach, the multivariate empirical mode decomposition, and the envelopes of each resulting intrinsic mode function (IMF) were computed by using Hilbert transform. Two IMFs, whose envelopes were strongly correlated with the HRV’s low-frequency oscillation (HRV-LF; ≈0.1 Hz) before and after the seizure were identified. The frequency ranges of these IMFs correspond to the EEG delta-band. The time-variant coherence was statistically quantified and tensor decomposition of the time-frequency coherence maps was applied to explore the topography-time-frequency characteristics of the coherence analysis. Results allow the hypothesis that couplings between the CAN, which controls the cardiovascular-cardiorespiratory system, and the ‘epileptic neural network’ exist. Additionally, our results confirm the hypothesis of a right hemispheric lateralization of sympathetic cardiac control of the HRV-LF.


Introduction
Epileptic seizure activity can cause alterations of the autonomic nervous system (ANS) in different ways and heart rate variability (HRV) is most frequently used to investigate both longterm and short-term alterations of the ANS in response to the type of epilepsy and to the evolution of the epileptic seizure [1,2]. HRV is an indicator for neuronal influences on the cardiac pacemaker and is one of the important functions of the ANS. Therefore, seizure-related HRV reactions may provide more information on the organization of the ANS [3] and the mechanisms supporting ANS changes.
The main divisions of the ANS are the parasympathetic and the sympathetic nervous systems, which typically function in opposition to each other. The coordination of 'autonomic functions and the ongoing behavioral needs of the organism through the activities of the somatomotor, endocrine and autonomic systems' [4] can be attributed to a central autonomic network (CAN) [5], which involves several interconnected areas distributed throughout brainstem and forebrain. According to the ANS's structure some of these CAN areas regulate parasympathetic output whereas others regulate sympathetic output, e.g. in order to control visceral function and homeostasis as well as to adapt to internal or external challenges [6].
From the point of view of network physiology [7], the temporal lobe epilepsy (TLE)related 'epileptic network' and the CAN can be considered as coupled networks. Taking this into account, couplings between EEG and HRV signals before, during and after an epileptic seizure can be considered as indicators of couplings between the 'epileptic network' and the CAN, i.e. this opens the possibility to study the time evolution of couplings between physiologic networks during an intense change of the physiological situation. By the term coupling a variety of interaction effects can be subsumed. Linear and nonlinear correlative approaches (linear correlation, spectral coherence, mutual information etc, e.g. [8]), methods for the analysis of directed interaction (Granger causality, transfer entropy, partial directed coherence etc, e.g. [9]), phase analysis techniques (time delay stability [10], phase-locking, n:m synchronization etc, e.g. [11][12][13][14]), and cross-frequency approaches (envelope-envelope, envelope-frequency, frequency-frequency coupling etc, e.g. [15][16][17]) can be used to quantify the specific type of physiological coupling in an appropriate way.
Our approach comprises the investigation of couplings between EEG and HRV signals before, during and after an epileptic seizure in children with TLE. Regarding the lateralization of the epileptic seizure during TLE, it has been shown that the heart rate (HR) is significantly increased in the preictal period (preictal tachycardia), when the focus of the TLE is in the right hemisphere [18]. No statistically significant changes could be observed in a left-focus group. Therefore, a right hemispheric lateralization of the sympathetic cardiac control can be assumed. Mayer et al [19] confirmed these results for children with TLE and stated a 'significant differences in HR evolution depending on location and side of seizure onset'. They found that an early and high HR increase was primarily associated with right hemispheric TLE. Consequently, our strategy is to divide the group of investigated children into a left hemispheric TLE ('left-focus group') and a right hemispheric TLE ('right-focus group') group.
In a previous study [20] we demonstrated that statistically significant coherence values between the HRV's low-frequency (HRV-LF) oscillation at about 0.1 Hz and the envelope of EEG delta-activity before and after the TLE seizure occur. Significant coherence values testify significant correlative coupling (synchronization) between two oscillatory signals of the same frequency. The findings that preictal and postictal EEG delta-activity is rhythmically modulated and that the modulation rhythm is correlated with HRV-LF were new at this time. However, time-variant coherence (tvCOH) analysis between the HRV and rhythmical EEG activity comprises a massive amount of results generated by the number of EEG channels, time points and frequency bins. Therefore, for feature extraction we defined regions of interest (ROIs) in the time-frequency range. The ROIs offer an opportunity to simplify the detection of relevant patterns, but reduce the benefit of our time-variant, frequency-selective analysis approach.
Our present study is based on two main methodological advancements. Firstly, we considered that the neurophysiological processes which contribute to an epileptic seizure are very individual, i.e. they depend on the type of seizure and for each patient the seizure process varies from seizure to seizure. Thus, empirical mode decomposition (EMD) was applied for the signal-adaptive selection of individualized frequency bands instead of fixed, group-related frequency bands as used in the previous study. As all the electrodes are considered to be a network, the multivariate empirical mode decomposition (MEMD) is used, i.e. the EEG is considered as a multivariate process. Multivariate approaches (multivariate autoregressive models) are usually used for connectivity analyses of the brain and of the cardiovascularrespiratory system [21,22].
Secondly, we extended the analysis strategy by using tensor decomposition of the results of tvCOH analysis which can be represented by a tensor. Tensor decomposition allows a reorganization and a dimension reduction of the analysis results. From the time-frequency analysis perspective, tensor decomposition is very attractive, because topographic, time, and frequency information (can be extended to trials, subjects, and groups) are simultaneously taken into account [23,24]. The results of tensor decomposition are used for the identification of topography-time-frequency profiles (patterns) of coherence analysis and in order to complement the topographic analysis on the basis of ROIs (mean coherence within subsequent time-intervals and over a defined frequency range). Additionally, the results of tensor decomposition can be used as a signal-adaptive feedback for the decision making of the ROIs in the time-frequency range.
The integration of both MEMD and tensor decomposition into a framework of a timefrequency analysis leads to an advanced processing and analysis concept. As expected we can show that the coupling effects are more pronounced and that additional information concerning the coupling profiles can be obtained.

Data material
The investigated data set consists of 18 children with TLE obtained during pre-surgical evaluation performed at the Vienna pediatric epilepsy center following a standard protocol. The protocol was approved by the local ethical committee of the University Hospital Vienna. For each patient, EEG (23 channels; 20 channels were used for our investigations) were recorded referentially from gold disc electrodes placed according to the extended 10-20 system with additional temporal electrodes. One-channel ECG was recorded from an electrode placed under the left clavicle. EEG and ECG data were recorded referentially against electrode position CPZ, filtered (1-70 Hz), converted from analog to digital (sampling frequency 256 Hz, 12 bit), and stored digitally for further data analysis. Further information concerning the classification of seizure type, onset and termination can be found in [25]. One recording per patient was analyzed containing 5 min before (preictal period) and 5 min after the seizure onset (ictal and postictal periods). The group of 18 patients was divided according to the side of the seizure's focus, i.e. patients with left or right hemispheric TLE, with 9 patients in each sub-group, referred to as 'left-focus' and 'right-focus group'.
After down sampling the EEG data from 256 to 64 Hz, an artefact removal procedure was applied, using the independent component analysis (ICA) provided by the Field Trip toolbox [26]. For each patient, the EEG signals were decomposed with ICA and the resulting components were visually inspected. The independent components (ICs) most affected by ocular or movement artefacts were removed and the signals were reconstructed using the remaining ICs. A detailed description of results of ICA decomposition of epileptic EEG can be found in [27]. That was followed by a re-referencing of the EEG to the average reference montage. Finally, 20 channels of EEG were selected for further analysis. QRS detection was performed after digital band pass filtering (10-50 Hz) of the ECG and an interpolation by cubic splines (interpolated sampling frequency 1024 Hz) to detect the time of the maximum amplitude of each R-wave. The resulting series of events was used for the HR computation, i.e., this series of events was low-pass filtered by means of a FFT-filter (cutoff frequency ⩽ half of the mean HR). The procedure is known as the French-Holden algorithm [28], which leads to the lowpass filtered event series (LPFES), a standard HRV representation in the time domain. The final HRV representation was obtained via multiplication of the LPFES with the sampling rate and with 60 beats per minute and down sampled to 8 Hz. An artefact rejection was performed manually to minimize the influence of false QRS triggering.

Methods
The analysis concept presented in figure 1 is based on and considerably advances our previous work [20]. In that study tvCOH between EEG envelopes of fixed frequency band components and the HRV was computed. In our new study, the fixed frequency band components are replaced with 'natural' frequency EEG components, obtained with the MEMD. Furthermore, tensor decomposition is adapted for individual or group results of achieved tvCOH.
2.2.1. Multivariate empirical mode decomposition. MEMD was introduced by Rehman and Mandic [29] as a multivariate extension of the standard EMD. EMD is a signal-adaptive algorithm for multiscale decomposition and time-frequency analysis of multicomponent signals [30]. EMD decomposes the multicomponent signal into a finite number of amplitude and frequency modulated zero mean oscillatory signals [31], called intrinsic mode functions (IMFs), and a monotonic residual. Details of the MEMD approach [29] can be found in the appendix section.
MEMD calculates the IMFs for all channels at the same time. Therefore, MEMD has the advantage of resulting in the same number of IMFs for all investigated channels, whereas EMD computation of the IMFs for the single channels results in a different number of IMFs for each EEG channel. This would lead to the correspondence problem, when a certain frequency band component has to be selected for further analysis. Using MEMD, for each patient the same number of IMFs is obtained for all channels. However, the number of IMFs might differ between patients, and for group analysis, IMFs were selected which are of interest for our analysis followed by a determination of the corresponding IMFs of different patients with similar frequency characteristics.
Time-invariant spectra of all IMFs were computed with the Fourier transform over 10 min epochs per patient and per channel. Consideration of power spectra is a standard method to investigate the properties of specific IMFs [32,33]. Distribution of the frequency characteristic was used to quantify the frequency range of IMFs by estimating the median as well as the 25th and 75th percentile of frequencies of the IMFs of each channel and each patient (see figure 3). We selected two IMFs that have the frequency characteristics most related to the delta band for further analysis.

Computation of IMF envelopes and tvCOH.
After selection of the two IMFs, their envelopes were computed for all channels in each patient using the Hilbert transform. The Hilbert transform gives the imaginary part of the analytical signal that is used to calculate the envelope. Having the IMFs, the tvCOH was computed pair-wise between the selected channelrelated IMF envelopes and the HRV, for each patient. For the computation of the coherence, the IMF envelopes were down sampled to 8 Hz to match the sampling frequency of the HRV; the time-variant spectra of the envelopes and HRV were obtained with the continuous Morlet wavelet transform (CMWT) (frequency range 0-0.5 Hz) and a time smoothing of the timevariant cross-spectrum and the two time-variant spectra was performed using a 512 points window. Details of the Hilbert transform and the computation of tvCOH can be found in [20] and in the appendix section.
For each patient, 20 (number of channels) tvCOH maps with the corresponding number of time points (10 min, sampling frequency 8 Hz) and frequency bins (0-0.5 Hz, 201 frequency bins) result. Having the tvCOH results for all patients, two directions of analysis were chosen. The first was a topographic representation of the mean value for the tvCOH in the 0.08-0.12 Hz frequency band and in subsequent 20 s intervals for all channels. The second was a statistical analysis of the mean value of the tvCOH in the 0.08-0.12 Hz frequency range (mean over frequencies) for specific electrode sites (T3 and T4). This step includes the computation of 1000 surrogate data by means of phase randomization [34] and the 95th (90th) percentile for the average tvCOH from patients of the same group (left-seizure/right-seizure) for each channel. A 5(10) % statistical threshold was computed as the mean over time of the 95th (90th) percentile.

Tensor decomposition.
Tensor decomposition techniques are used for the reorganization and reduction of the dimension of our tvCOH results. Using parallel factor analysis (PARAFAC) [35,36], the tensors were decomposed using a number of 2-4 factors. The PARAFAC algorithm decomposes the tensor into a sum of multi-linear terms (factors) without imposing orthogonality constraints, giving a unique solution of the PARAFAC model under very mild conditions [37]. The contribution of particular factors at each level of a particular mode of the dataset is given by the weights of the factors [36]. Tensor decomposition is uniquely identifiable up to scaling and permutation. Therefore, ordinate of each factor is in [a.u.]. Details of tensor decomposition algorithm can be found in the appendix section.
For each patient and both IMFs of interest the tvCOH results were arranged as a 3rd order tensor having the structure 'topography' × 'time' × 'frequency'. For group analysis, a 4th order tensor with the structure 'topography' × 'time' × 'frequency' × 'patient' was created for the two investigated groups and both IMFs.

MEMD decomposition and selection of relevant IMFs
For each patient, the MEMD was applied to the 20 channels, resulting in the same number of IMFs for each channel. An example of the first five IMFs for one patient and one channel is presented in figure 2. The purpose of the decomposition is to obtain the 'natural frequency' components of the EEG. To identify the IMFs that correspond to the EEG frequency bands (e.g. beta, alpha, theta and delta) and/or the individual seizure-related frequency characteristic of each patient, the associated Fourier spectrum was computed and analyzed. Figure 2 illustrates that the frequency bands are overlapping and that the IMFs do not correspond exactly to the standard EEG frequency bands, e.g., the spectrum of the 4th IMF (figure 2(D)) contain the frequency band 2-7 Hz with maximum around 3.5 Hz and represents a combination of delta (1.5-4 Hz) and theta (4-7 Hz) activity. The spectrum of the 5th IMF (figure 2(E)) have the frequency band 0.5-3.5 Hz with maximum around 2 Hz and is mainly associated with delta-activity.
Since the focus of this study is on the interaction between EEG delta-activity and HRV, the 4th and 5th IMF were selected for further analysis.
As already mentioned, using the MEMD approach, the correspondence between IMFs of different channels is not an issue. However, correspondence between patients has to be confirmed. In order to select the same 'natural frequency components' for all patients, the spectra of the two selected IMFs were compared. In figure 3, the distribution of the frequencies for two selected channels and the two IMFs (IMF4 and IMF5) for all patients is characterized by the median value and the 25th and 75th percentiles ( figure 3(A)). It can be seen that the selected IMFs have a consistent frequency band distribution over all patients. Figure 3(C) presents the median value, the 25th and 75th percentiles of the frequency distribution for the two IMFs and all channels for one patient (#02). A detailed representation of the spectra of the two channels and IMFs for the same patient is depicted in figure 3(B). The values for all channels are comparable for both IMFs. This allows the selection of the two IMFs for the computation of the corresponding envelope in the next processing step. Here, the time-invariant investigation of spectra of the IMFs serves as a simple tool to discriminate between different frequency ranges. For more detailed analysis it can be replaced by time-frequency representations achieved by CMWT.

Results of the tvCOH
The envelopes of the two selected IMFs were computed for all channels and all patients by means of the Hilbert transform. Examples of signals that were used in the study are depicted in figure 4. The EEG signal of one channel (A), the resulting IMF4 (B) and IMF5 (D) and their envelopes (C, E) along with the HRV (F) are shown. The two envelopes and the HRV rise immediately after the seizure onset (300 s). The HRV starts to decrease approximately 50 s after toward the preictal mean HRV value. The same decrease after 50 s can be observed for the envelope of the IMF4. However, the postictal values for both envelopes remain higher than the pretictal levels.
In our previous study, we found significant coherence patterns between the delta-envelope and the HRV in the 0.08-0.12 Hz frequency band. A similar approach was followed for this study. First, the tvCOH between the IMF-envelope corresponding to one electrode position and the HRV was computed in the 0-0.5 Hz frequency band. One example of tvCOH for a rightfocus patient is presented in figure 5(A). From the tvCOH, only the 0.08-0.12 Hz band was selected and the time was divided into subsequent 20 s intervals. By computing the mean in this restricted time and frequency interval one value was obtained, representing the mean coherence between the IMF-envelope of one channel and the HRV. Computing the same for all channels a topographic representation of the mean coherence in the selected frequency band and time interval ( figure 5(B)) was obtained, resulting in 26 topographic representations (interval from 40 to 560 s was used) for one patient. This representation offers the possibility to roughly analyze the topographic changes in time of the coherence patterns.
A view of a restricted interval (180-400 s) is presented in figure 6 as 11 subsequent topographic representations for the mean coherence values between the HRV and IMF4envelope (A) and IMF5-envelope (B) for one patient (#02) of the right-focus group (mean over time and frequency).
The coherence patterns are comparable to the ones obtained for the coherence between HRV and the delta-envelope [20]. For this patient, high coherence values indicate a strong correlative coupling (synchronization) between HRV and the two envelopes before the seizure onset. Starting with the seizure (300 s) the coherence values are very low for almost all channels, and a recuperation of the coupling starts after 20 s for some electrodes. We previously identified two electrodes (C4, T4) that have a consistent coupling pattern over a longer period of time. In the current study, it seems that the two electrodes have the same consistent pattern (designated by white frame). However, C4 has higher coherence values for IMF4 and T4 for IMF5. This could indicate that the two electrodes are correlated with the HRV at different frequencies (different activities).  The results for the group-related tvCOH analysis at two electrode positions (T4 and T3) and the statistical threshold are presented in figure 7. The 5 (10)% threshold was computed as the mean over time of the 95th (90th) percentile obtained with surrogate data.
The results for the right-focus group at electrode site T4 (on the side of the focus) of the mean tvCOH between the HRV and the IMF5-envelope (black) and the IMF4-envelope (gray) show significantly high values around 150-250 s (seizure onset at 300 s). For the left-focus group, the tvCOH at the same electrode position (opposite to the focus side) presents high coherence values that exceed the 5% threshold in the postictal period 480-580 s. These results are consistent with our previous ones. The mean tvCOH exceeds the 5% threshold for the other electrode position (T3), however, the respective time intervals are different for the two envelopes.
The individual characteristics of the epileptic seizure (severity, focus localization, activity spreading) impairs a correct assessment of a group analysis following this approach, since smearing effects can appear from averaging the coherence values over the entire group of patients (left-focus or right-focus group). Also, it may be that the mean value in the restricted frequency and time interval is also not truly representative for such a dynamic process as correlation. In order to analyze all results from the tvCOH, for each patient and to perform a group analysis without averaging the results, tensor decomposition was considered. Group-related results of tvCOH analysis for the electrode sites T3 (left) and T4 (right) achieved by averaging over the corresponding 9 children. Seizure onset is designated by the dashed vertical line. An overlay of tvCOH between HRV and IMF4envelope (gray) and IMF5-envelope (black) is represented. The 5% threshold is designated by the red line and the 10% threshold by green. The mean tvCOH for the right-focus group for electrodes T3 (A) and T4 (B), and for the left-focus group for T3 (C) and T4 (D) is shown. Abscissa are in (s).

Tensor decomposition of tvCOH
The advantage of tensor decomposition is that tvCOH can be analyzed at the same time for all channels without restricting the frequency band and/or the time interval. Thus, a first step was to create for each patient a 3rd order tensor with modes 'topography' (channel), 'time' and 'frequency' using the results for all channels of the tvCOH between the HRV and IMF4envelope and IMF5-envelope respectively, in the frequency range 0-0.5 Hz. An example of tensor decomposition for one right-focus patient is presented in figure 8.
The decomposition of the tvCOH between HRV and IMF4-envelope (A) and IMF5envelope (B) in a number of three factors was sufficient in order to obtain the topography-timefrequency profiles of the tvCOH that are related to the 0.1 Hz HRV-LF oscillations which are described by the second factor (b) for both decompositions. The clear peak around 0.1 Hz, the time profile with higher weights in the preictal period and the topography mode which includes electrodes on the side of the seizure suggest that the second factor describes best the patterns of the tvCOH between the EEG activity and the HRV-LF components before the seizure onset.
However, the frequency mode of the second factor also includes very low frequency components with weight as high as the 0.1 Hz component. The question arises whether the low frequency components influence the profiles of time and topography modes. Therefore, we restricted the tvCOH in the frequency band around the 0.1 Hz and repeated the tensor decomposition.
In figure 9, the tensor decomposition of the tvCOH between the HRV and the IMF4envelope (A) and IMF5-envelope (B) in the restricted 0.05-0.15 Hz frequency band is presented. The tensors were decomposed in a number of 2 factors, in order to extract the HRV-LF (0.1 Hz). The first factor of both decompositions has the time and frequency profiles of the tvCOH in the frequency range of the HRV-LF. The time weights are higher in the preictal period and are comparable to the group results for the right-focus group. The topography mode indicates the electrode positions that are the most related to the time and frequency profiles and it can be seen that the selected electrodes for the group analysis T3 and T4 are included. The topography-time-frequency profiles of the second factor are also consistent with the profiles of the second factor in the previous decomposition ( figure 8). However, the more refined analysis shows better results.
This proves that the tensor decomposition can be beneficially used to extract the topography-time-frequency profiles of the tvCOH.
Including all patients that belong to the right-focus group results in a 4th order tensor with modes 'topography', 'time', 'frequency' and 'patient'. This allows a group analysis for all channels using tensor decomposition without averaging the tvCOH results. The results of the tensor decomposition of the tvCOH between the HRV and the IMF4-envelope and IMF5envelope for the right-focus group are presented in figure 10. Two factors were used for both tensor decompositions.
The second factor resulting from the decomposition of the tvCOH between HRV and IMF4-envelope (A) has a frequency mode corresponding to the 0.1 Hz of the HRV-LF that is associated with the Mayer waves of the arterial blood pressure (ABP). The time mode has a high peak approximately 50 s before the seizure onset (300 s) that starts to decrease with the beginning of the seizure. The topography mode indicates the electrode positions most related to the time-frequency patterns and the patient mode suggest that the patterns are not associated to all patients.
The decomposition of the tvCOH between HRV and IMF5-envelope (B) results in a circumscribed 0.1 Hz peak (frequency mode) represented by the second factor. The patterns of the associated time and topography modes are comparable to the group results in figure 7. However, the patient-related weights of the patient mode indicate that the patterns are not consistent for all investigated patients.

Discussion
We have shown that correlative couplings between the (modulated) delta-activity of the EEG and the HRV-LF can be observed in the preictal and postictal periods which speak for a coupling of the seizure-associated cortical processes ('epileptic network') and the CAN. The CAN controls the Figure 10. Group analysis of the tvCOH results for the right-focus group by means of tensor decomposition. (A) Tensor decomposition of the group results for the tvCOH between HRV and IMF4-envelope using a number of two factors. The frequency mode of the second factor indicates that the associated time and topography profiles are most related to the 0.1 Hz (HRV-LF). The time mode has a very high peak approximately 50 s before the seizure onset (designated by the red line at 300 s) that starts to decrease with the seizure onset. The patient mode indicates that this specific topography-timefrequency profile is not associated with all patients, and is more specific to only two patients. (B) Second factor of the tensor decomposition for the tvCOH group results between HRV and IMF5-envelope. (B) represents the best demonstration of the tvCOH patterns in the 0.1 Hz frequency range and confirms the choice of the frequency band for the ROI (0.08-0.12 Hz) in the time-frequency range. The time and topography profiles are comparable to the group results obtained in figure 7, however are not found in all patients. Ordinates are in (a.u.). cardiovascular-cardiorespiratory system(s), i.e. these systems together with the 'epileptic network' form a higher-level, integrative network. Bashan et al [38] argued that 'the human organism is an integrated network where complex physiological systems, each with its own regulatory mechanisms, continuously interact, and where failure of one system can trigger a breakdown of the entire network'. The sudden unexpected death in epilepsy (SUDEP) can be seen as such an 'entire breakdown'. Possible causes for SUDEP are related to cardiac factors (e.g. arrhythmias), respiratory factors (e.g. suppression of respiratory centers in the brain stem), an ANS dysregulation, etc. Therefore, our results contribute to both to the development of the new field 'network physiology' [7] and to the elucidation of the mechanisms which cause SUDEP.
Further discussion is subdivided into three main parts: (1) a comparison of our results to results of other HRV-EEG coupling studies, (2) a discussion of the physiological mechanisms which are possibly involved in the coupling chain between HRV and EEG, and (3) a discussion of the methodological aspects of our analysis strategy.

HRV-EEG coupling studies
HRV-EEG coupling studies are mostly related to sleep research (polysomnographic data). There are studies which show a time-variant correlation between frequency band activity and HRV, e.g. during paradoxical sleep (between HRV and EEG delta-theta bands [39]) and across sleep stages (EEG spectral edge frequencies and HRV-LF (Pedemonte, Rodriguez-Alvez et al 2005)). Additionally, Jurysta et al [40] demonstrated that a closed connection between cardiac autonomic activity and spectral EEG bands exists. The delta band shows the highest variations in response to high-frequency HRV components and ANS activity precedes changes in the EEG during sleep in healthy young men. In one of our previous studies [41] a correlation between quadratic phase couplings courses of the EEG and HRV before and during EEG burst activity during quiet sleep (NonREM sleep in preterm and fullterm newborns) in preterm newborns was shown. It is known that preterm neonates have a deficit in ANS activity and a sympatheticparasympathetic imbalance characterized by sympathetic predominance exists [42]. The synchronous changes of EEG and HR we have discussed as an indication for a coupling between cortical, thalamocortical and central autonomic brain areas.
The polysomnographic analysis concept of P Ch Ivanov widens the scope further. His analysis strategy aims at the identification of 'a network of interactions between diverse physiological organ systems, to quantify the hierarchical structure and dynamics of this network' [7]. This higher-level network consists of ten nodes representing six physiological systems: brain activity (five EEG frequency bands), cardiac system (HR), respiration (respiratory movements), chin muscle tone, leg and eye movements.
As only EEG and ECG recordings were used in our studies, a broader investigation to test Ivanov's network approach was not possible. However, our results will aid in further refinement of the analysis methodologies for 'physiological networks'.

Physiological mechanisms
Between HRV-LF and EEG activity modulation a chain of subsystems with own regulatory mechanisms and mutual interactions to other systems has to be considered. A discussion should start with the phenomenon ('enigma' [43]) of Mayer waves in systemic ABP and their connection with HRV-LF. Mayer waves are strongly correlated with the oscillations of efferent sympathetic nervous activity. In fact these waves actually result from oscillations of the sympathetic arterial smooth muscle tone [44] and 'are almost invariably enhanced during states of sympathetic activation' [43]. Mayer-wave associated HRV-LF oscillations include most probably both sympathetic and parasympathetic (vagal) influences [45].
The baroreflex plays a major role in the generation of Mayer waves. In TLE patients the baroreflex function is chronically impaired, e.g. the transfer function gain between ABP and HRV-LF is reduced which determines the baroreflex function [46,47]. Other studies have shown that TLE patients are characterized by a dysfunction of the cardiovascular autonomic regulation (autonomic instability [3]), manifesting as impaired HR responses to certain stimuli [48]. These impaired regulatory mechanisms and the resulting instabilities might be responsible for the augmentation of the HRV-LF before and after the seizure (acute influence) which we have demonstrated in another study [25]. This augmentation (enhanced phase-locking) takes place during the time periods in which we have found the statistically significant coherence values. And vice versa, when the HRV-LF is suppressed during the seizure [25] the coherence values are non-significant ( figure 7).
The CAN can be directly influenced by a TLE seizure. It is known that seizure activity which arises from or propagates to areas in the CAN (e.g. cortical limbic areas which are connected with subcortical regions) can mimic stimulation of autonomic afferents [49] by which ANS functions can be altered. Activation of the amygdala, which is often the source or associated with mesial TLE seizures, increases the sympathetic activation. It is also well-known that the HR can increase before the TLE seizure (e.g. [50]). In our data such a HR increase can also be detected, in particular in the right-focus group [19,25]. However, in the preictal period premonitory brain processes must exist, that support such coupling effects between the HRV-LF and delta-envelope. It was recently shown that focal hemodynamic changes, i.e. cerebral blood flow (CBF) increases and hemoglobin oxygenation decreases, precede seizure onset by approximately 20 s [51,52]. These changes could be elicited by subtle neuronal or glial events, astrocyte-or pericytemedicated signaling or local potassium, and local neurotransmitter/neuropeptide release. Severe postictal disturbances (dysregulations) of the ANS over a time range of 5-6 h are described by Toth et al [53] (HRV analysis). Therefore, it is not surprising that we found stronger postictal couplings in both groups in comparison to those of the preictal period.
The neurovascular coupling, i.e. the relationships among neuronal activity, metabolism, tissue oxygenation, and CBF, is another mechanism which can contribute to such phenomena [51]. Additionally, cerebral auto-regulation, the process by which CBF is kept at a constant level, as well as sympathetic regulation of CBF may also participate in these effects. It is known that Mayer waves in systemic ABP create variations in cerebral blood flow velocity (CBFV) in the intracranial arteries of the same frequency [54]. This establishes a link between Mayer waves associated HRV-LF and CBF. Chen et al [55] investigated dynamics of CA regulation from the relationship between the peripheral blood pressure in the finger and the CBFV for a group of healthy and post-stroke subjects during different physiologic conditions. This partial synopsis of functional chains and regulatory mechanisms demonstrates the high complexity of interactions between systems (networks) and subsystems. Lastly, a physiological aspect which remains to be addressed is the influence of the TLE focus side. Jansen and Lagae [3] noted that 'due to the hemispheric specific organization of the central ANS, autonomic symptoms in epileptic seizures can provide lateralizing and localizing information'. It has been shown that the HR already significantly increases in the preictal period (preictal tachycardia), if the focus of the TLE was in the right hemisphere [18]. In contrast, no statistically significant changes could be observed in a left-focus group. The authors stated that a right hemispheric lateralization of the sympathetic cardiac control can be assumed. These results were confirmed for children with TLE. Our results for the right-focus group confirm that the coupling effects between the HRV-LF and delta-envelope depends on the lateralization of the seizure.

Methodological aspects
The advanced analysis used in this study yields a clear improvement of our previous results, i.e. the coupling effects are more pronounced which allows a more confident interpretation of the results. Our previous results can be confirmed and new results have been added. The EEG of the interictal and preictal periods in TLE children is typically characterized by spikes or sharp-wave discharges and intermittent rhythmic delta-activity [56]. Our detailed IMF analysis reveals that the two IMFs, which envelopes show couplings with HRV-LF, correspond to delta-activity. An additional advantage of EMD approaches is that nonlinear characteristics of the signals components are preserved in the corresponding IMFs [25]. Yeh et al [57] used EMD to select low-and high-frequency oscillations from sleep EEG at different states and applied nonlinear measures for their quantification.
Tensor decomposition (we used the PARAFAC) provides a reorganization and dimension reduction of the results of tvCOH analysis (time-frequency coherence maps). Thus, a better interpretation of the results of time-frequency analyses is possible [24]. The topography-timefrequency-profile of single factors provides information with regard to the time-frequency characteristic of coupling and allows a feedback concerning the selection of ROIs (timefrequency ranges). Our results confirm the choice of ROI definition and of the coupling dynamic which was identified by statistical analysis.

Conclusion
It can be stated that this methodological concept can be beneficially applied for the analysis of tvCOH between HRV-LF and EEG-envelopes. By means of MEMD the individualized, seizurerelated EEG frequency bands can be extracted and by using tensor decomposition the topographytime-frequency profile of couplings for each patient of both groups can be revealed.
By means of this signal-adaptive analysis concept coherence patterns can be demonstrated which indicate interactions between the 'epileptic network' and the CAN before and after the seizure. These couplings speak for chronic and/or acute impairments (instabilities) of complex physiological regulation mechanisms caused by the epilepsy or the epileptic seizure. Additionally, our results confirm the hypothesis of a right hemispheric lateralization of sympathetic cardiac control of the HRV-LF.
Our study provides a methodological strategy which can also be beneficially applied to coupling analyses between other physiological networks. It is our hope that this work will provide the incentive.
For a multivariate K-component signal, , and a set of , , , 6. Extract the 'detail' d(t) = x(t) − m(t). If d(t) respects the stopping criterion for a multivariate IMF, the procedure is repeated for x(t) − d(t), otherwise for d(t). Since 'extrema' cannot be properly defined for multivariate signals, the stopping criterion doesn't impose the condition of same number of extrema and zero crossings. Even so, the stopping criterion for multivariate IMFs is similar to the one used in EMD. Mode-mixing is a problem also with MEMD. A graphical representation of a third-order tensor decomposition as a sum of two factors is given in figure A3.1.
The PARAFAC algorithm decomposes the tensor into a sum of multi-linear terms without imposing orthogonality constraints. The solution of the PARAFAC model is unique under very mild conditions [35,37]. One condition was proposed by Kruskal [61] for threeway arrays: The basic method used for fitting the PARAFAC model is alternating least squares. The loadings in two modes are successively assumed to be known, and the unknown set of parameters for the last mode are estimated. The algorithm start with the selection of the number of factors and the initialization of matrices B and C. The matrix A is estimated from X, B and C by least squares regression. The next steps is to estimate the matrix B likewise, and then matrix C in the same manner. The process is repeated until the model converges to a solution (little change in fit or loadings). A more detailed description of the algorithm and improvement suggestions can be found in [35,37]. Constrains like non-negativity, unimodality and orthogonality can improve the speed of the algorithm [62]. The aspects regarding the factor number selection is addressed in [35], where one suggestion is to compare the results with external knowledge about the data that are modelled (like spectral or temporal patterns). It was pointed out that when extracting too many components (factors), the noise is being increasingly modelled and the true factors are modelled by more correlated components.