Brain–heart interactions considering complex physiological data: processing schemes for time-variant, frequency-dependent, topographical and statistical examination of directed interactions by convergent cross mapping

Background: A multitude of complex methods is available to quantify interactions in highly complex physiological systems. Brain–heart interactions play an important role in identifying couplings between the central nervous system and the autonomic nervous system during defined physiological states or specific diseases. The crucial point of those interaction analyses is adequate pre-processing taking into account nonlinearity of data, intuitive graphical representation and suitable statistical evaluation of the achieved results. Objective: The aim of this study is to provide generalized processing schemes for such investigations taking into account pre-processing, graphical representation and statistical analysis. Approach: Two defined data sets were used to develop these processing schemes. Brain–heart interactions in children with temporal lobe epilepsy during the pre-ictal, ictal and post-ictal periods as well as in patients with paranoid schizophrenia and healthy control subjects during the resting state period were investigated by nonlinear convergent cross mapping (CCM). Surrogate data, bootstrapping and linear mixed-effects model approaches were utilized for statistical analyses. Main results: CCM was able to reveal specific and statistically significant time- and frequency-dependent patterns of brain–heart interactions for children with temporal lobe epilepsy and provide a statistically significant pattern of topographic- and frequency-dependent brain–heart interactions for schizophrenic patients, as well as to show the differences from healthy control subjects. Suitable statistical models were found to quantify group differences. Significance: Generalized processing schemes and crucial points of pre-processing, adapted interaction analysis and performed statistical analysis are provided. The general concept of analyses is transferable also to other methods of interactions analysis and data representing even more complex physiological systems.

Brain-heart interactions, that is interactions between the cortical activity of the central nervous system (CNS; e.g. represented by the electroencephalogram (EEG)) and the autonomic nervous system (ANS; e.g. represented by specific characteristics of heart rate variability (HRV)), have become another emerging field of analysis. In general, evaluations of brain-heart interactions are able to provide complementary information beyond the important and still relevant results of univariate analysis of the central and the autonomic nervous system itself. Information achieved by univariate analyses of both systems form the basis of interaction analyses in terms of the pre-definition of signal characteristics like e.g. relevant frequency content and time patterns. Brain-heart interactions have been extensively analyzed, e.g. by means of polysomnography/during sleep (Faes et al 2014. Beissner et al (2013) conducted a meta-analysis of human neuroimaging experiments evaluating central autonomic processing. The major goal of this analysis was to localize cortical and subcortical areas involved in autonomic processing, to find potential subsystems for the sympathetic and parasympathetic divisions of the ANS, and to find potential subsystems for specific ANS responses to different stimuli/tasks.
Clearly, directed interaction measures taking into account the complex and possibly nonlinear behavior of data to be analyzed will necessarily persist in this area. There is a wide range of different approaches to define directed interactions in complex systems. The most prominent methods are based on the concept of Granger causality (GC); multivariate approaches (Milde et al 2010) including time-variant and frequency-selective adaptions like partial directed coherence (PDC) are available (Leistritz et al 2013). Furthermore, model-free measures of information transfer like transfer entropy (TE) (Schreiber 2000, Vicente et al 2011 or the concept of timedelay stability (TDS)  exist to quantify linear and nonlinear interactions in complex physiological systems. Furthermore, other metrics like predictability-based measures exit to perform between-systems analyses, like e.g. predictability decomposition taking into account self, causal and interaction predictability (Faes et al 2016) or within-systems analyses like the evaluation of complexity and causality by comparing model-based and model-free approaches (Porta et al 2014).
Convergent cross mapping (CCM), exemplarily adapted in our study, was introduced by Sugihara et al (2012) and defines interactions between time series in terms of nonlinear stability and predictability. Independent of the actual approach adapted for the quantification of (linear/non-linear) interactions, it has become increasingly essential for all of these investigations (1) to properly pre-process the data to be analyzed, (2) to statistically evaluate the complex pattern achieved by the interaction analysis, and (3) to be able to visualize complex interactions in an adequate manner, also taking into account statistical results.
Therefore, the aim of this study is to show defined and general processing schemes for interaction analyses including those crucial points. Development of such processing schemes are demonstrated by two complex investigations of interaction analysis by means of CCM: quantification of brain-heart interactions in children with temporal lobe epilepsy during the pre-ictal, ictal and post-ictal periods of seizure, and quantification of brain-heart interactions in a group of patients with paranoid schizophrenia and healthy control subjects during the resting state period. Both applications include the analysis of interactions between HRV and specific EEG components.
The main advantage of looking at both application is that (1) from the technical point of view different features of analysis like time variance, frequency selectivity, topography and/or group comparison are necessary to proceed, and (2) from the methodological point of view advanced information exist about complementary linear/non-linear and directed/non-directed interaction analyses in the data , Schiecke et al 2016, Schulz et al 2016, 2018.
Analysis of brain-heart interactions in children with TLE is an extension of a previous study investigating non-linear interactions only in a limited frequency range (Schiecke et al 2016). Non-directed, linear brain-heart interactions in the same group of children were investigated in Piper et al (2014). Interactions between the CNS and ANS during epilepsy are possibly partly also induced by changes of each of the involved single systems. A general overview of cardio-respiratory variability during TLE in children is given in Varon et al (2015); a timevariant, frequency-selective HRV analysis taking into account linear and non-linear measures in the specific group of children with TLE was performed in Schiecke et al (2014). Characterizations of a specific EEG pattern during TLE can be found in Williamson et al (1993). In general, the main goal of interaction analyses for epileptic patients is (1) to get more information concerning the causes of the sudden unexpected death in epileptic patients (SUDEP), (2) to predict epileptic seizures and (3) to detect the focus of a seizure.
Cardiac autonomic dysregulation in patients with acute schizophrenia which is possibly associated with an increased risk of cardiac arrhythmias and sudden cardiac death is well described, e.g. in Bar et al (2005), Zahn and Pickar (2005) and Bar et al (2010). Former analysis of brain-heart interactions in schizophrenic patients included linear and nonlinear analysis of short-term central-cardiorespiratory couplings (Schulz et al 2018) and central-autonomic couplings achieved by an adaptation of PDC and methods of joint symbolic dynamics (Schulz et al 2016).
In section 2 proposed processing schemes and main methodological considerations can be found; in section 3 both investigated data sets were introduced and application-adapted processing schemes are given. An extensive description of results including appropriate visualizations and statistical evaluations can be found in section 4, followed by a discussion of achieved concepts in section 5.

Methods
First, in this section a general processing scheme was developed. Accordingly, all methods used for pre-processing, data analysis, statistical evaluation and graphical representation are introduced.

General processing scheme
Our general processing scheme developed for the analysis of directed brain-heart interactions performed on complex physiological data is presented in figure 1, left. A more detailed schematic view of all relevant information for general pre-processing by multivariate empirical mode decomposition (MEMD) and the accordingly adapted CCM approach can be found in figure 1, right.
While (I) the performed pre-processing in case of brain-heart interactions will be nearly the same for a wide range of applications, there are clear differences in (II) the adaptation of CCM analysis, (III) statistical evaluation and (IV) graphical representation depending on the complexity of investigated physiological system (figure 1 right, e.g. need of time variance, frequency selectivity, topological views, group comparison…).
Adapted processing scheme for CCM analysis, statistical evaluation and graphical representation of both applications will be given in section 3.

Pre-processing by means of MEMD
MEMD was introduced by Rehman and Mandic (2010). MEMD extends the standard EMD procedure to a procedure applicable to multivariate data. In general, EMD is a signal-adaptive algorithm for multiscale decomposition and time-frequency analysis of multicomponent signals (Huang et al 1998) and has, in comparison to linear filtering, the clear benefit of preserving the nonlinearity of the investigated data. EMD decomposes a multicomponent signal into a finite number of amplitude-and frequency-modulated zero mean oscillatory signals, called intrinsic mode functions (IMFs), and a monotonic residual (Flandrin et al 2004).
MEMD decomposes multivariate multicomponent signals (e.g. a specific number of EEG channels) into IMFs and results in the same number of IMFs for each of the specific set of investigated channels. Thus, there is no problem of assignment of corresponding IMFs over channels, but for an assignment of corresponding IMFs over different subjects an automated approach is preferable. Standard criteria for the investigation of the properties of specific IMFs are consideration of power spectra. Usually, time-invariant spectra of all IMFs are computed with the Fourier transform per subject and per channel and compared according to the frequency content. Adaptation of standard Kuhn-Munkres algorithm can be used for such kind of assignments, details and examples of its application can be found in Schiecke et al (2015b).

Interaction analysis by means of convergent cross mapping (CCM)
The CCM approach was introduced by Sugihara et al (2012). The basic idea of CCM is to test for causation between time series X and Y in a nonlinear way. For that, the correspondence between 'shadow manifolds' constructed from lagged coordinates of the time series values of X and Y are used. Looking at the causality concept of CCM, the main issue is the contrast to intuition and particularly the concept of GC: when causation is unilateral (e.g. 'X drives Y'), then it is possible to estimate X from Y, but not Y from X. Practically, interaction from time series X to Y are quantified by CCM correlation, defined by the absolute value of the Pearson correlation coefficient between the original time series X and an estimation of X using the CCM with Y (and interactions from Y to X are quantified by the one between the original time series Y and an estimation of Y using the CCM with X, respectively).
For the determination of bivariate, directed CCM both values CCM X→Y and CCM Y→X are compared for any pair of time series X and Y. Practically, CCM X→Y and CCM Y→X are examined using data sets X = x(t) and Y = y (t) of increasing data length L (the 'library length'). In the bidirectional case ('X drives Y stronger than vice versa') CCM X→Y will converge faster/reach a higher plateau than CCM Y→X respectively. For an investigation of time-varying directed interactions between time series X and Y, an interval-based estimation of CCM can be performed by using a sliding window of an appropriate library length L. Details of adapted estimation routine, practical details of the basic algorithm as well as performed simulations can be found in Sugihara et al (2012), Schiecke et al (2015a) and Schiecke et al (2016).

Surrogate data approach
A surrogate data approach was used to find a threshold to quantify statistically significant interaction values. The null hypothesis of the surrogate data approach is that there is no directed interaction from time series X to Y, and from time series Y to X, respectively ('small directed values of CCM'). For 'causal surrogates', the causal interaction is destroyed only in one direction, leaving the degree of association in the opposite direction untouched (Faes et al 2010). Surrogate data were obtained by destroying the Fourier phases but keeping the power spectra of investigated signals by means of phase randomization for each direction of interaction (Theiler et al 1992).
Accordingly, CCM estimations was performed on the causal surrogates. Based on this, statistical thresholds for significant interactions values result.

Bootstrapping approach
A bootstrapping approach was used to determine statistical properties of CCM estimated for a defined number of realizations.
Statistical properties of CCM values are of interest regarding (1) CCM values achieved by time-varying interval-based approaches and (2) for CCM values achieved over all investigated EEG channels. For the estimation of tubes of confidence intervals of the extracted CCM courses without any particular distribution assumption, 1000 bootstrap samples of size K were drawn from the estimated CCM, respectively. A case resampling with replacement was used.
One thousand bootstrap replications of each extracted parameters were computed. Based on these replications, the lower limit (the 2.5% quantile) defined the lower bound, and the upper limit (the 97.5% quantile) defined the upper bound of the tubes of 95% confidence intervals of estimated CCM. Theoretical details of bootstrapping approaches can be found in Efron and Tibshirani (1994)

Mixed-effects model approach
For statistical evaluation of group differences with the CCM as dependent variable a linear mixed-effects model was performed (Pinero and Bates 2000). Linear mixed-effects models offer the advantage of allowing the investigation of variability between subjects (heterogeneity) and simultaneously adjusting for the withinsubject correlation due to the existence of several measurements/results per subject. 'EEG electrode', 'direction of interaction' and 'EEG frequency ranges' are examples of relevant properties for CCM analyses. Finally, and based on the linear mixed model approach, the estimated means for the respective categories together with their 95% confidence intervals result. Figure 1. General processing scheme for directed analysis of brain-heart interactions. Necessary steps within (I) pre-processing, (II) data analysis, (III) statistical evaluation, and (IV) graphical representation are given on the left side; a more detailed schematic view of the general pre-processing by MEMD and adapted CCM analysis is given on the right hand side.

Graphical representation of achieved results
Adapted graphical representations are necessary to reflect different features of analysis like time variance, frequency selectivity, topography and/or group comparison.
In case of a time-variant implementation of CCM and a small number of investigated electrode positions results can be represented in the form of (1) time-frequency maps of mean CCM values per position and direction of interaction (e.g. using MATLAB, color code: CCM values between 0 and 1; no inclusion of statistical results possible) and (2) time courses of mean CCM values, e.g. per position and per frequency range of interest achieved by MEMD (inclusion of tubes of confidence intervals and statistical thresholds possible).
In case of a time-invariant implementation of CCM and a large number of investigated electrode positions results can be presented in the form of (1) topographic maps of mean CCM values per group, per direction of interaction and per frequency range of interest achieved by MEMD (e.g. using Fieldtrip, color code: CCM values between 0 and 1; statistical results can be shown by separate topographic maps using a value of 1 for significant and value of 0 for non-significant CCM) and (2) electrode-based views of mean CCM values per group investigated, per direction of interaction and per pre-defined frequency range (inclusion of tubes of confidence intervals and statistical thresholds possible).
Any further complexity of data to be analyzed (e.g. time-variant implementation of CCM, large number of investigated electrodes, multiple frequency ranges and group-related investigation) would imply an advanced summarization and graphical presentation of results, like e.g. tensor decompositions (Pester et al 2015).
Graphical view of results of linear mixed-effects model approach is based on the representation of mean values over categories together with the corresponding confidence intervals (e.g. using MATLAB).

Subjects and data analysis
Two different data sets were used to demonstrate the adapted processing routines, performed statistical evaluations and adapted graphical representations for the analysis of brain-heart interactions taking into account different features like time variance, frequency selectivity, topography and/or group comparison of achieved results.

Recordings and data material
The first data set originates from a group of 18 children with temporal lobe epilepsy (11 female, seven male, median age 9.4 years, range from 6.5-18 years, median seizure length 88 s, range 52-177 s) obtained during presurgical evaluation performed at the Vienna pediatric epilepsy center following a standard protocol. Information concerning the classification of seizure type, onset and termination can be found in Mayer et al (2004). The protocol was approved by the local ethical committee of the University Hospital Vienna. For each child, EEGs (23 channels) were recorded referentially from gold disc electrodes placed according to the extended 10-20 system with additional temporal electrodes. One-channel ECGs were recorded from an electrode placed under the left clavicle. EEG and ECG data were recorded referentially against electrode position CPZ (sampling frequency 256 Hz) and filtered (1-70 Hz). Recordings were performed during the pre-ictal, ictal and post-ictal periods.

Processing scheme and performed analyses
In general, one recording per child was analyzed, including 5 min before and 5 min after seizure onset, representing pre-ictal, ictal and post-ictal periods.
Pre-processing of EEG data included a down sampling of the 23 EEG channels from 256 Hz to 64 Hz and an artefact removal procedure using independent component analysis (ICA). Twenty channels of EEG data were used for further analysis. MEMD was performed for all EEG channels to be able to differentiate between defined and individually adapted frequency ranges resulting in sets of IMFs per channel and child. Four specific IMFs which are of interest for our analysis were selected (IMF(β-EEG); IMF(α-EEG); IMF(θ-EEG); IMF(δ-EEG)). After selection of the IMFs, their envelopes were computed using the Hilbert transform. The sampling frequency of the envelopes was 8 Hz. Details of the adapted MEMD approach and calculation of the envelopes of EEG-IMFs can be found in Piper et al (2014).
Pre-processing of ECG data included QRS detection 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. The series of events was low-pass filtered by means of an FFT filter (cutoff frequency ⩽ half of the mean HR, French-Holden algorithm; see (French and Holden 1971)). Multiplication of the low-pass filtered series of events with the sampling rate and with 60 beats per minute (bpm) and a downsampling to 8 Hz resulted in the final HRV representation.
The main characteristics of the CCM analysis, statistical evaluations and graphical representation of data set 1 (children with TLE) can be found in figure 2.
An interval-based CCM approach of the whole recording (10 min; pre-ictal, ictal and post-ictal period; sliding windows of L = 256 data points = 30 s) was performed and resulted in a bivariate, time-variant and frequency-selective view of directed interactions between HRV and EEG envelopes.
The statistical analysis of the achieved CCM results included the following: • a surrogate data approach to find a threshold to quantify statistically significant interactions; • a bootstrapping approach to quantify statistical properties of directed CCM time courses estimated for K = 18 realizations of the investigated data.
For each investigation performed (see processing scheme in figure 2: two EEG channels, four IMFs, two directions), 50 surrogate data sets per realization (K = 18 children) were obtained resulting in 900 repetitions per situation. CCM estimations were carried out for the library length L used for interval-based CCM estimation (L = 256 data points). The 95th percentile of the surrogate data's CCM was computed and served as the threshold for statistically significant CCM values per investigated situation.
Statistical properties of CCM values are of interest regarding CCM values achieved by the time-varying interval-based approach. For the estimation of tubes of confidence intervals of the extracted CCM courses without any particular distribution assumption, 1000 bootstrap samples of size K = 18 were drawn from the estimated time courses of CCM achieved by the interval-based approach.
The graphical representation of the achieved CCM results included: • time-frequency maps of CCM means per electrode position (ipsilateral and contralateral) and per direction of interaction; • time courses of CCM means with tubes of confidence intervals and statistical threshold per electrode position (ipsilateral and contralateral) and per pre-defined frequency range.   (First 1997). Psychotic symptoms were quantified using the Positive and Negative Symptom Scale (PANSS) (Kay et al 1987). Schizophrenic patients were treated with antipsychotic medication (77% being atypical neuroleptics and 23% being a mixture of antidepressant and atypical neuroleptics). For healthy control subjects, in-depth interviews and clinical investigations were performed to exclude any potential psychiatric or other diseases and to check for any interfering medication. The patients and control subjects did not receive any medication that might have influenced the cardiovascular system. All participants gave written informed consent to a protocol approved by the Ethics Committee of the University Hospital, Jena, Germany. Patients were only included after a psychiatrist certified their ability to give full consent to the study protocol. For each patient and control subject ECG (three channel, sampling rate 500 Hz) and EEG (64 channel, sampling rate 500 Hz) were recorded synchronously for 15 min. The EEG was obtained using active Ag/AgCl electrodes, and transmitted via a BrainAmp amplifier (Brain Products, Germany, ground AFZ, reference FCZ). Electrodes were positioned by using an electrode cap and according to the extended 10-20 system. Recordings were started after a resting period of 10 min. Participants were asked to close their eyes, relax, and breathe normally to avoid hyperventilation. No further breathing instructions were given.

Processing scheme and performed analyses
In general, one 15 min recording per patient and per control was analyzed during the resting state period. The main characteristics of the CCM analysis and statistical evaluations of patients with schizophrenia and healthy control subjects can be found in figure 3. The adapted processing scheme and visualization of the structure of the analyzed data is given.
Pre-processing of EEG data included MEMD of all 64 channels per subject followed by an automated assignment of corresponding IMFs across all subjects using Kuhn-Munkres algorithms (Schiecke et al 2015b). Only corresponding IMFs across subjects were inspected according to their frequency content. Specific IMFs of interest to our analysis were selected (IMF(γ-EEG); IMF(β-EEG); IMF(α-EEG); IMF(θ-EEG); IMF(between θ and δ); IMF(δ-EEG)). After the selection of the IMFs, their envelopes were computed using the Hilbert transform. The sampling frequency of the envelopes was 5 Hz. The pre-processing of ECG data included QRS detection performed after digital band pass filtering (10-50 Hz) of the ECG to detect the time of the maximum amplitude of each R wave. The resulting series of events was used for the HR computation. The series of events was low-pass filtered by means of an FFT filter (cutoff frequency ⩽ half of the mean HR, French-Holden algorithm; see (French and Holden 1971)). Multiplication of the low-pass filtered series of events with the sampling rate and with 60 beats per minute (bpm) and a downsampling to 5 Hz resulted in the final HRV representation.
The main characteristics of the CCM analysis, statistical evaluations and graphical representation of data set 2 (patients with schizophrenia and healthy control subjects) can be found in figure 3.
A time-invariant, bivariate CCM analysis of the whole recording (15 min; during the resting state period) was performed and resulted in a bivariate, topographic and frequency-selective view of the directed interactions between HRV and EEG envelopes.
The statistical analysis of the achieved CCM results included the following: • a surrogate data approach to find a threshold to quantify statistically significant interactions; • a bootstrapping approach to quantify statistical properties of directed CCM topographies estimated for K (K = 17 and K = 21) realizations of the investigated data; • a linear mixed-effects model approach to quantify group differences between schizophrenic patients and healthy control subjects.
For each investigation performed (see scheme in figure 3: 64 EEG channels, six IMFs, two directions), 30 surrogate data sets per realization of patients (K = 17) and 24 surrogate data sets per realization of controls (K = 21) were obtained, resulting in 510 repetitions per situation for the patient group and 504 repetitions per situation for the control group, respectively. CCM estimations were carried out for the library length used for time-invariant CCM estimation (L = 4500 data points). The 95th percentile of the surrogate data's CCM was computed and served as the threshold for statistically significant CCM values per investigated situation.
Statistical properties of CCM values are of interest regarding CCM values achieved over all investigated EEG channels. For the estimation of tubes of confidence intervals of the extracted CCM courses without any par-  Concerning the mixed-effects model approach, an intercept only model for the random effects based on individual subjects was applied. Inclusion of random effects for EEG electrodes did not further improve the fit of the model. For this reason, the intercept-only model was kept. The covariates 'direction' and 'EEG frequency bands' and their statistical interactions were included.
The graphical representation of the achieved CCM results included the following: • topographic maps of CCM means per group (control and patient) and per direction of interaction and per pre-defined frequency range; • topographic maps of significant CCM means per group (control and patient) and per direction of interaction and per pre-defined frequency range; • electrode-based views of CCM means with tubes of confidence intervals and statistical threshold per group (control and patient), per direction of interaction and per pre-defined frequency range; • estimated means for the respective categories together with their 95% confidence intervals per group (control and patient), per direction of interaction and per pre-defined frequency range.

Results
The results are presented for both investigated data sets. Estimated CCM values as well as results of performed statistical evaluations are given in figures 4 and 5 for the children with TLE and in figures 6-9 for patients with schizophrenia and healthy control subjects.

Children with TLE
The adapted specific processing scheme for the evaluation of brain-heart interactions by means of non-linear CCM in children with TLE was able to show specific time and frequency patterns of those interactions. Time-frequency maps of mean CCM values per electrode position and direction of interaction were able to clearly present an overview of results of interval-based, bivariate CCM analysis between HRV and pre-defined EEG-IMFs (figure 4).  The results of linear mixed-effects model approach. Estimated means for the respective categories together with their 95% confidence intervals (CI) are presented for patients (round marker) and controls (squared marker) per six pre-defined frequency ranges (top down: gamma to delta EEG range) and per direction of interactions (blue: CCM EEG-IMF→HRV , red: CCM HRV→EEG-IMF ). Significant differences between patients and controls are denoted by an asterisk (corresponding non-overlapping CI are marked by a black elliptical frame).
There was not much difference between ipsilateral (figure 4, left) and contralateral (figure 4, right) positions, but there was a clear difference between both directions of interaction (figure 4, upper and lower line) and also between pre-and post-ictal periods (figure 4, all subplots, more pronounced for higher frequencies of interest like β and α). The statistical evaluation of achieved results is not possible by means of time-frequency maps like the one presented in figure 4.
Time courses of mean CCM values achieved by interval-based, bivariate CCM analysis between HRV and per pre-defined EEG-IMFs (four frequency ranges, figure 5) were also able to include results of statistical evaluation by means of tubes of confidence intervals achieved by bootstrapping and statistical thresholds for significant interactions achieved by causal surrogates.
In general, the time courses of interactions between HRV and all pre-defined EEG-IMFs (figure 5) reveal higher values of CCM EEG-IMF→HRV (blue for ipsilateral and green for contralateral location) than vice versa (CCM HRV→EEG-IMF , red for ipsilateral and black for contralateral location, accordingly) during the pre-ictal period. The values of CCM HRV→EEG-IMF increase sharply for all pre-defined EEG-IMFs with seizure onset, mostly reach the values of CCM EEG-IMF→HRV , stay increased during the ictal period and decrease again in the post-ictal period. No clear difference in the time courses is visible by looking at the ipsilateral and contralateral locations.
Looking at the absolute values of CCM during the pre-ictal period, the greatest interaction are found between EEG-IMF(δ) and HRV, followed by EEG-IMF(θ). Less interactions are found for EEG-IMF(α) and EEG-IMF(β) having a higher variance over time and subjects (see time courses and CI-tubes achieved by bootstrapping).
Other differences between the four investigated frequency ranges are visible during the ictal and post-ictal periods. While the sharp increase during the ictal period is very short for interactions between EEG-IMF(α) and HRV as well as EEG-IMF(β) and HRV, there are longer lasting changes for interactions between EEG-IMF(δ) and HRV as well as EEG-IMF(θ) and HRV. Furthermore, CCM values of interactions between EEG-IMF(δ) and HRV as well as EEG-IMF(β) and HRV are clearly higher during the post-ictal period than during the pre-ictal period. In particular CCM HRV→EEG-IMF(δ) and CCM HRV→EEG-IMF(β) do not return to the original values estimated during the pre-ictal period.
The statistical significance of differences in the time courses of directed interactions (CCM EEG-IMF→HRV versus CCM HRV→EEG-IMF ) could be proven during the pre-and post-ictal periods only, and for EEG-IMF(δ), EEG-IMF(θ) and EEG-IMF(α). These significant differences are evaluated by means of non-overlapping tubes of 95% CI achieved by bootstrapping (figure 5). No statistically significant differences in the time courses of directed interactions could be shown for EEG-IMF(β) and during the ictal period (overlapping tubes of 95% CI, figure 5).
The statistical significance of the interaction (threshold by surrogate data) could be shown for all CCM EEG-IMF→HRV and for CCM HRV→EEG-IMF(δ) and CCM HRV→EEG-IMF(θ) . The statistical significance of the interaction is given by tubes of 95% CI not covering the statistical threshold achieved by surrogates. No statistical significance of the interaction could be shown for CCM HRV→EEG-IMF(α) and CCM HRV→EEG-IMF(β) (tubes of 95% CI cover statistical threshold, figure 5).

Patients with schizophrenia and healthy controls
The adapted specific processing scheme for the evaluation of brain-heart interactions by means of non-linear CCM in patients with schizophrenia and healthy controls was able to show specific group differences, frequency patterns and topographic differences of those interactions. Topographic maps of time-invariant, bivariate CCM analysis between HRV and pre-defined EEG-IMFs (six frequency ranges) for patients and controls as well as the corresponding results of the applied mixed-effects model approach to quantify group differences are given in figures 6 and 7. The tubes of confidence intervals achieved by bootstrapping, the statistical thresholds for significant interactions and a topographic view of quanti fied significant interactions are given in figures 8 and 9.
CCM analysis (figure 6) reveals stronger values for CCM EEG-IMF→HRV than CCM HRV→EEG-IMF over all predefined EEG-IMFs and for controls as well as patients. All CCM values are higher (and reveal stronger interactions) for controls than for patients with the exception of CCM HRV→EEG-IMF(θ) and CCM HRV→EEG-IMF(δ). There, interactions are stronger for patients than for controls. A topographic view of interactions shows a clear pattern appearing in the different investigated frequency ranges, which is different for patients and controls. For EEG-IMF(γ), there are peaks of interactions in the temporal areas (more pronounced for the controls); for EEG-IMF(β) a similar but weaker pattern appears (only left-sided for the patients); for EEG-IMF(α) the assumed occipital peaks of interactions can be found (nearly not visible for patients); no clear pattern but an increasing strength of interaction is visible for EEG-IMF(θ) to EEG-IMF(δ) (less pronounced again for patients).
The mixed-effects model approach reveals significant group differences between patients and controls for CCM EEG-IMF(β)→HRV , for CCM EEG-IMF(α)→HRV and for CCM HRV→EEG-IMF(δ) (figure 7;. In all other categories, no significant difference between patients and controls could be proven (figure 7; overlapping CIs).
The statistical properties of CCM achieved by bootstrapping in combination with statistical threshold for significant interactions over all investigated EEG channels (figure 8) also reveal clear differences between patients and controls.
Less-pronounced differences in topographic pattern of CCM values of patients are visible, but on the other side, patients show in most of the investigated situations a higher variance of absolute CCM values (figure 8; broader tube of 95% CI achieved by bootstrapping for patients).  Results achieved by surrogated data analysis over all 64 channels for controls (upper lines) and patients (lower lines), respectively, are shown per investigated frequency range (left to right: gamma to delta EEG range) and both directions of interaction including all 64 channels. The black color designates the significance of interaction with a value of 1; the white color shows no significance of interaction, with a value of 0.
The statistical significance of interaction in controls (figure 8: threshold by surrogate data; figure 9: topographic view of significant interactions achieved by surrogated data) could be shown for all CCM EEG-IMF→HRV and for CCM HRV→EEG-IMF(δ) (figure 8: tubes of 95% CI do not cover statistical threshold; figure 9: black color) but not for all channels for all other CCM HRV→EEG-IMF (figure 8: tubes of 95% CI sometimes cover statistical thresholds; figure 9: some white color).
The statistical significance of the interaction in patients (figure 8: threshold by surrogate data; figure 9: topographic view of significant interactions achieved by surrogate data) could be shown again for all CCM EEG-IMF→HRV and for CCM HRV→EEG-IMF(δ) (figure 8: tubes of 95% CI do not cover statistical threshold; figure 9: black color = value of 1) but not for all channels for all other CCM HRV→EEG-IMF . In particular, CCM HRV→EEG-IMF(β) , CCM HRV→EEG-IMF(α) and CCM HRV→EEG-IMF(θ) does not reveal statistically significant interactions for most of the channels (figure 8: tubes of 95% CI cover almost all statistical thresholds; figure 9: mostly white color = value of 0).

Discussion and conclusion
The main focus of this study was to present a generalized processing scheme for the investigation of brain-heart interactions of complex physiological systems including appropriate pre-processing, statistical evaluation necessary for such complex situations and adequate graphical presentation of results. Two example data sets were employed to achieve this aim.
From the methodological point of view, crucial point of adequate pre-processing is to avoid all linear filtering. Florin et al (2010) showed that the application of filters may significantly affect connectivity patterns and should thus be avoided as pre-processing steps, e.g. for GC analysis. Therefore, the application of MEMD, especially preserving the nonlinearity of data is an important solution for this problem and should be used to achieve frequency-selective results taking into account the nonlinear properties of the investigated data. If using MEMD in a large number of subjects and/or situations, then an automated assignment of resulting IMFs is necessary. An adaptation of the Kuhn-Munkres algorithm is available to solve this problem (Schiecke et al 2015b).
Time-variate and multivariate generalizations of most of linear and non-linear interaction methods are available. CCM has the disadvantage of having only a bivariate application but it works well in a time-varying manner, whereas e.g. TE is multivariate in its original implementation but is afflicted with the problem of the data length needed for estimation (Wollstadt et al 2014). Nevertheless, one important generalization of this investigation is that introduced processing schemes can independently be applied to any other linear or non-linear interactions measures like GC-based approaches or TE. Time-variant and frequency-selective views of the achieved results (data set 1) as well as topographic and frequency-selective views of achieved results sorted by group definition (data set 2) were given. Any further complexity of data to be analyzed and thus the results achieved would imply other summarizing and presenting of results, like tensor decompositions. The use of tensor decomposition to present complex results of PDC analyses is given in Pester et al (2015) and Mierau et al (2017).
Bootstrapping and surrogate data approaches have proven to be essential for statistical evaluation of complex results of interaction analysis and are necessary to be able to statistical quantify interactions in highly datadependent situations and are able to easily present statistical significance. An adequate null hypothesis for a surrogate data approach should involve the use of 'causal surrogates' (Faes et al 2010). In particular group comparisons taking into account various (dependent) covariates require complex statistical methods like e.g. linear mixed-effect models.
The physiological interpretations of achieved results for both investigated data sets were not the focus of this study and would possibly interfere with the clear methodological aim of this work. Nevertheless we were able to show that the investigation of brain-heart interactions in children with TLE was able to clearly demonstrate a statistically significant time-and frequency-dependent pattern of interactions. The extension to a broader frequency range of investigation in comparison to former study (Schiecke et al 2016) allowed for further differentiation between influencing factors and show the high complexity of underlying processes. Again, no differences between ipsilateral and contralateral location could be proven. The 'epileptic network' rather than 'epileptic focus' concept discussed in Lehnertz and Dickten (2015) may be a reason for that finding.
Investigations of brain-heart interactions in schizophrenic patients were able to provide a topographic view using all 64 EEG channels by differentiating among relevant frequency ranges. Physiologically reasonable topographic patterns and relevant statistically significant differences in this pattern between patients and controls could be identified. Further investigations should follow possibly using a larger group of patients or subgroups of patients looking at specific aspects of CNS-ANS couplings like specific HRV values and pattern of patients with schizophrenia and possible influence of that on brain-heart interactions.
In general, any linear and non-linear analysis of brain-heart interactions should be accompanied by a previous, adequate univariate view, e.g. of the HRV or EEG pattern, to be able to quantify characteristics of the central and autonomic nervous system during a specific physiological or pathophysiological state.