The detection of transient directional couplings based on phase synchronization

We extend recent approaches based on the concept of phase synchronization to enable the time-resolved investigation of directional relationships between coupled dynamical systems from short and transient noisy time series. For our approach, we consider an observed ensemble of a sufficiently large number of time series as multiple realizations of a process. We derive an index that quantifies the direction of transient interactions and assess its statistical significance using surrogate techniques. Analysing time series from noisy and chaotic systems, we demonstrate numerically the applicability and limitations of our approach. Our findings from an exemplary application to event-related brain activities underline the importance of our method for improving knowledge about the mechanisms underlying memory formation in humans.


Introduction
Characterizing the interactions between two or more (sub-)systems is an important but challenging issue in many scientific fields, ranging from physics to neuroscience [1]. Although interactions cannot be measured directly, they can nevertheless be quantified by applying appropriate time series analysis methods to system observables. In recent decades, a number of techniques that are based on information or synchronization theory have been proposed to characterize the strength or direction of an interaction (see [2]- [4] for an overview). More recently, these techniques were supplemented with methods based on the Fokker-Planck formalism [5,6]. Most approaches assume the underlying systems to be stationary, the coupling between them to be stable, and the time series to be composed of a sufficiently large number of data points. Thus, a short intermittent coupling cannot be resolved accurately with these techniques, and fast dynamical changes are usually regarded as an undesirable complication. Nevertheless, transient coupling might constitute an even more interesting aspect of interacting dynamical systems.
As a prominent example from neuroscience, we mention here the investigation of so-called event-related potentials (ERPs). The electroencephalogram (EEG) reflects electrical neural activity due to intrinsic dynamics and/or responses to external stimuli. To examine the pathways and time courses of information processing in the brain under specific conditions, numerous experiments have been developed controlling sensory or other inputs (stimuli) [7,8]. The neural activity induced by this kind of stimulation leads to potential changes in the EEG. These ERPs often exhibit a sequence of multi-phasic peak amplitudes after stimulus onset and extending over a time period of several 100 ms. ERPs reflect different stages of information processing in the brain, and insights into the underlying neural processes can already be achieved from the analysis of specific aspects of ERPs, such as peak amplitude or moment of occurrence (latency). Since individual ERPs possess very low peak amplitudes, as compared to the ongoing EEG, multiple realizations with respect to a given stimulus are commonly averaged, assuming phaselocked responses not correlated with the ongoing EEG. More recently, data analysis techniques were proposed that make use of only a few measurements instead of large preprocessed ensembles [9]- [13].
The human brain is a complex network of interacting subsystems, and it is now commonly accepted that synchronization plays an important role in brain function [14]- [17]. In order to further improve insights into the mechanisms of interactions subserving information processing in the brain there is a strong need for reliable analysis techniques that allow one to characterize both the strength and direction of an interaction from short and transient signals such as ERPs.
While the strength of an interaction only quantifies the degree of mutual influence, knowledge about the direction of an interaction is essential to investigate the influence of one system on another. For example, knowing the direction of influences between the neocortex and the hippocampus is crucial for a better understanding of mechanisms of transfer of memory contents between these structures. While a number of approaches allow one to investigate the strength of interactions [18]- [22], only a few methods have been proposed for the estimation of the direction of interactions from transient signals [23,24]. These methods are based on an embedding of the data into some high-dimensional state space and require the parameters embedding dimension as well as time delay to be chosen appropriately [25,26], which might be nontrivial when analysing time series from systems with only poorly understood properties. Moreover, since the phases of signals in particular are sensitive to weak couplings, a phase synchronization-based method for the detection of directed interactions between stationary systems, as proposed in [27], can be regarded as a more promising approach. We present here an extension of this approach that allows the time-resolved investigation of the direction of interactions from short and transient noisy time series.

Time-resolved detection of directional couplings
Given the time series of observables from systems X and Y , the phase time series [φ X,Y ( jδt) = φ X,Y (t j )] j=1, ... ,n can be extracted using either the Hilbert or a wavelet transform [28,29]. Following [27], an asymmetric influence of system X onto system Y can be estimated by modelling the underlying phase dynamics, which are considered to be generated by unknown with noisy perturbations ζ X,Y .
With the slope of the phases derived from the time series, the deterministic part F can be modelled using finite Fourier series F. In order to estimate the coefficient vector x, the so-called design matrix A can be derived, and the linear least square problem for approximating F can then be written as where b denotes the vector of phase slopes for all t j . Now, in order to resolve directional couplings on short time scales, we consider an observed ensemble of m time series as multiple realizations of some process of transient interactions induced by some event. We rewrite the phase time series as [φ X,Y (t r j )] r =1, ... ,m j=1, ... ,n with the ensemble index r , and j = 1 denotes the time index of the onset of coupling. For each time point, the Fourier series F X j for system X then reads using the following combination of coefficients: |α| 3 for β = 0, |β| 3 for α = 0, and |α| = |β| = 1 [30]. To simplify matters we define [F r jk ] r =1, ... ,m j=1, ... ,n, k=1, ... ,d for each Fourier term k, each realization r , and each time point j, and obtain a time-dependent version of the design matrix: 4 We define the time-dependent phase slope vector b j via the phase increments φ where τ denotes an appropriately chosen time delay. Following [31], we here choose τ = min(T X , T Y ) with the mean periods T X,Y of the phase time series. The timedependent phase slope vector then reads Given these definitions, (2) can now be solved for each time index j, and the influence of system Y onto system X is quantified by The sum on the rhs can be derived analytically [32]. With the influence of system X onto system Y , which is defined in complete analogy, we quantify a time-dependent directional relationship as D j (X, Y ) is not bounded, but by definition it will attain positive values if system X drives system Y and negative values for the opposite case. Note that C j can attain nonzero values even for uncoupled systems, e.g. due to a similar temporal evolution of the systems. In order to minimize false interpretations of directed couplings, we developed a surrogate-based testing scheme that is related to the methods proposed in [23,33,34]. Under the null hypothesis that the underlying phase dynamics of both systems are independent, the surrogates are designed such that they preserve the temporal structure of each time series and destroy the dependence between systems. To test if the process of system X is independent of that of system Y , we employed the following surrogate test. We randomly permute the realizations r of phase variable φ Y (t r j ) and calculate [C p j (X |Y )] p=1, ... ,s with the permuted realizations of system Y . With s = 39 random permutations and a significance level of 97.5%, C j (X |Y ) is set to zero only if C j (X |Y ) max p (C p j (X |Y )). The testing for C j (Y |X ) is defined in complete analogy. In the following, we refer to the surrogate-tested quantities asC j andD j .

Numerical examples
To illustrate our approach, we present findings obtained from an analysis of time series from coupled noisy phase oscillators [27] with natural frequencies ω X = 1.1 and ω Y = 0.9. Here X,Y denotes the coupling strength, and ξ X,Y are Gaussian δ-correlated noisy perturbations with amplitude h. The oscillators (8) were integrated using the Euler-Maruyama method [35] with a step size of 0.009π and sampling interval δt = 0.113. With initial conditions chosen randomly from [0, 2π) we generated m = 1000 realizations, each consisting of n = 1200 data points. The transient diffusive couplings between oscillators were realized such that X = 0.1 ∧ Y = 0 for In parts (a) and (b) of figure 1, we show C j (X |Y ) and C j (Y |X ), as well as the maximum and minimum values of the corresponding surrogates. As expected, C j (X |Y ) and C j (Y |X ) increase shortly after the onset of the coupling and remain well above the maximum values of the surrogates until the coupling is turned off. Likewise, the directional relationship reflects these coupling changes (cf figure 1(c)); we observe, however, D j (X, Y ) = 0 also for uncoupled oscillators. In contrast,D j = 0 except for statistical outliers (cf figure 1(d)). During the intervals with a constant nonzero coupling, both D j (X, Y ) andD j fluctuate around some level. From an analysis of 20 sets with m = 1000 realizations each we derive the means D j and standard deviations (cf figure 1(e)) and conclude that these fluctuations can-to a large extent-be related to the noisy perturbations. indicates system X to predominantly drive system Y and D j < 0 indicates system X to be driven by system Y . Time intervals with unidirectional couplings (W Y →X and W X →Y ) are highlighted in light grey.
Next, we proceed with analysing interactions between nonidentical nonlinear oscillators. We here chose Lorenz oscillators [36] as a prototypical example: with R X = 35, R Y = 39, and the nonzero couplings X = 4.0 for j ∈ W Y →X and Y = 4.0 for j ∈ W X →Y (cf [23]). The equations of motion were integrated with sampling rate δt = 0.01 using an algorithm based on the Livermore solver for ordinary differential equations [37,38].
To generate 20 sets of m = 1000 realizations each, we chose the initial conditions randomly in the state space near the Lorenz attractor and estimated phase time series from the x X,Y variables via the Hilbert transform using the analytical signal approach [39,40]. The temporal evolution of D j clearly reflects the coupling changes (cf figure 2) but in contrast to the noisy oscillators the range of absolute values ofD j differs for the intervals W Y →X and W X →Y , which can be related to the different control parameters. When comparing figures 1(e) and 2 we note that the range of absolute values of the directional relationship varies for the unidirectionally coupled oscillators investigated here. Thus, an appropriate normalization ofD j may be expedient when analysing empirical data. We also notice thatD j does not immediately indicate a change in the coupling, which can be related to the definition of the phase slope vector (cf (5)) and coupling properties of the oscillators.
In order to estimate the performance of our approach, particularly with respect to the analysis of empirical data, we investigated the impact of noise contaminations and the number of realizations m on the correct detection of directional relationships. We repeated the analyses for coupled Lorenz oscillators as outlined above for different numbers of realizations m and varied the coupling strength X,Y , as well as the amount of noise contamination (additive Gaussian δ-correlated noise; the signal-to-noise ratio (SNR) is σ 2 s /σ 2 n ). Then, we accumulated correct (γ C = j∈W − Y →X |D j | + j∈W + X →Y |D j |) and incorrect indications of directional relationships (γ I = j∈W Y ↔X |D j | + j∈W + Y →X |D j | + j∈W − X →Y |D j |). With γ = γ C /(γ C + γ I ), which is confined to the interval [0,1], we define a performance estimatorγ as the mean from 20 sets of oscillators. In the upper parts of figure 3, we show the dependence ofγ on the coupling strength X,Y and the number of realizations m. For coupling strengths X,Y > 2.5, we observeγ 0.75 even for a number of realizations as low as 200. For coupling strengths X,Y ≈ 1.0, an increase of the number of realizations (m > 1000) suffices to reliably detect directional relationships.
In the lower parts of figure 3 we show, for a fixed coupling strength X,Y = 4.0, the dependence ofγ on SNR and on the number of realizations m. Assumingγ 0.75 to indicate a sufficient performance, directional relationships can be resolved with m > 100 realizations for noise contaminations with SNR 10. As expected, even stronger noise contaminations require a larger amount of realizations (m 3000) in order to achieve a comparative performance.
Before closing this section, we conclude that our approach, which is based on the concept of phase synchronization, allows one to investigate-in a time-resolved manner-the direction of interaction from short and transient noisy time series. Our numerically obtained findings indicate that directional interactions can be detected reliably with only a few hundreds of realizations.

Detection of transient directional couplings in event-related potentials
In this section we present exemplary findings obtained from applying our method to ERPs that were recorded intracranially from an epilepsy patient who underwent presurgical evaluation ERPs recorded from the different brain structures during the encoding and the retrieval condition. We averaged the raw EEG time series in a time interval extending from 300 ms before stimulus onset to 1000 ms after stimulus onset, and corrected amplitudes with respect to the mean amplitude from the prestimulus interval, which is usually taken to be a period of inactivity. The dashed vertical line at j = 0 indicates stimulus onset.
of drug-resistant temporal lobe epilepsy. The patient had signed informed consent that his/her clinical data might be used and published for research purposes, and the study protocol had previously been approved by the local ethics committee. In order to define the zone of seizure origin for resective surgery, depth electrodes had been implanted bilaterally along the longitudinal axes of the hippocampi [41] to enable recording of brain electrical activities with high signal-to-noise ratio and with high spatial resolution. For this patient seizures proved to originate unilaterally from the left temporal lobe; thus electrodes located in the right brain hemisphere enabled recordings of activities unrelated to epilepsy. We analysed ERP data that were recorded while the patient performed two memory-related tasks (cf [42,43]). During an encoding phase 896 pictures (either houses or faces) were presented sequentially for 2500 ms on a computer screen with an inter-stimulus interval of 1500 ms. In order to ensure that the patient adequately attended to and processed each item, she/he was asked to rate each as pleasant or unpleasant by pressing one of two buttons. During the retrieval phase 576 pictures (either houses or faces) were again presented, but this time the sequence consisted of 384 old pictures that had already been presented, and another 192 new presentations. The patient was asked to react to all pictures by pressing one of four buttons, indicating whether she/he judged a given picture to be a first presentation (sure or unsure) or a repeat (sure or unsure). While the patient performed these tasks, brain electrical activities were recorded from multiple electrodes within the frequency band 0.1-300 Hz with a sampling rate of 1000 Hz and using a 16-bit analogue-to-digital converter. Linked mastoid electrodes served as reference. Discarding trials with pronounced artifacts resulted in 738 trials for the encoding phase and 434 trials for the retrieval phase.
For our exemplary analyses, we concentrated on ERP recordings from the hippocampus and the rhinal cortex. These brain structures are known to be crucially involved in memory processing [19,44,45]. Using magnetic resonance imaging scans that had been performed after implantation of the electrodes, we identified those electrode contacts that best sampled activities from these brain structures.
In figure 4, we show-for each brain structure-averaged ERPs that we derived for each condition separately in the time interval ranging from 300 ms before stimulus onset to 1000 ms 9 after stimulus onset. The time interval ranging from 0 to 1000 ms after stimulus onset is sufficient to capture the main ERP components (hippocampus: a positive peak at around 200 ms and a negative peak at around 500 ms; rhinal cortex: a negative peak at around 400 ms) that have been described in previous studies [46,47]. We observe here similar ERP components, and the clear amplitude differences seen for the different conditions already indicate that the rhinal cortex and the hippocampus are involved in the memory processing steps investigated here. Although some findings indicate a more dominant influence from rhinal cortex to hippocampus during encoding and an inverse influence during retrieval [44,48], averaged ERPs do not provide clear-cut information as to whether there are transient directed interactions between these brain structures. Moreover, we are not aware of any study that has investigated such interactions in a time-resolved manner.
We therefore performed a time-and frequency-resolved analysis of directional relationships between the rhinal cortex and the hippocampus using the raw EEG time series from the trials for each condition as described above. We concentrated on the well-known EEG-frequency bands θ Note that particularly activities in the θ -and γ -band are associated with encoding and retrieval of memory [19], [48]- [50]. After bandpass filtering (Butterworth characteristic) each trial forward and backward to account for phase shifts, the data were normalized to zero mean, which corresponds to setting the dc Fourier coefficient to zero. To avoid edge effects, we tapered each trial using a cosine half-wave (Hanning window) before deriving phase time series using the Hilbert transform. Since the calculation of the Hilbert transform requires integration over infinite time, which cannot be performed for a time series of finite length, we discarded 1/8 of the calculated instantaneous phase values on each side of every trial. We then calculated D j andD j as described in section 2. Since these analyses yield 15 600 values ofD j (6 frequency bands, 1300 points per time interval, 2 conditions), we applied a Bonferroni correction to account for multiple comparisons. We adjusted the significance threshold by dividing it by the total number of statistical tests and obtained a corrected probability value forC j as α B = 0.025/15 600 = 1.6 × 10 −6 . Assuming that directed interactions, if they exist, extend over at least some tens of milliseconds, we considered a sequence of at least four consecutive values ofC j not equal to zero to indicate statistically significant ((0.025) 4 < α B ) directed interactions. Following our findings in section 3, we rescaled each time series D j and (Bonferroni-corrected)D j such that it is confined to the interval [−1, 1], and denote them by D * j andD * j , respectively. In figure 5, we show the time-resolved directional relationships between the rhinal cortex and the hippocampus for the two experimental conditions and for the investigated frequency bands. Except for the α-band, we observe significant directional influences around 200 ms and in the period 400-1000 ms after stimulus onset in all other frequency bands. In the θ-band, we observe for the encoding condition that the rhinal cortex first influences the hippocampus (at about 130 ms), while the influence is then reversed later on (at about 500 ms). This sequence appears to be reversed for the retrieval condition; here the hippocampus first influences the rhinal cortex (at about 270 ms) and the influence is then reversed later on (at about 580 ms). Although our findings fit quite well to current hypotheses of memory formation in humans [51,52] we here refrain from providing an interpretation. Future studies on a larger group of subjects may elucidate the possible neurofunctional relevance of directed interactions between different brain structures for memory formation in humans.

Conclusion
We have proposed an approach based on the concept of phase synchronization which allows the time-resolved investigation of directional relationships between coupled dynamical systems from short and transient noisy time series. Given an ensemble of such time series as multiple realizations of some interaction process, we have extended the approach proposed by Rosenblum and Pikovsky [27] and derived an index that quantifies the direction of a transient interaction. We assessed the statistical significance of our index with a surrogate-based testing scheme [23,33,34]. With several numerical examples that are representative of noisy and chaotic systems that can be nonidentical, we have exemplified the applicability of our approach. Moreover, we have provided an estimate of the minimum amount of realizations needed to reliably identify directional relationships for different coupling strengths and for different contaminations with observational noise. Our preliminary findings obtained from a time-and frequency-resolved analysis of directional relationships between ERPs recorded from the human rhinal cortex and hippocampus underline the potential impact of our approach on improving the understanding of neural processes underlying memory formation and other cognitive processes.
As one of the directions for further study, we consider improvements from using more sophisticated methods for estimating the phase from raw time series [53,54], as well as the use of other weighting functions [55,56]. As another direction, we mention comparisons of our approach with other quantifying measures for the direction of a transient interaction [23,24]. We expect that our approach will also find applications in other scientific fields where transient couplings are considered an interesting aspect of interacting dynamical systems.