We apply the method of frequency-selective excitation waves in excitable media to characterize synchronization phenomena in interacting complex dynamical systems by measuring coincidence rates of induced excitations. We relax the frequency-selectivity of excitable media and demonstrate two applications of the method to signals with broadband spectra. Findings obtained from analyzing time series of coupled chaotic oscillators as well as electroencephalographic (EEG) recordings from an epilepsy patient indicate that this method can provide an alternative and complementary way to estimate the degree of phase synchronization in noisy signals.
1. Introduction
Excitable media (EM) are spatially distributed excitable systems that can rapidly propagate excitations (suprathreshold impulses) over long distances without damping. Such traveling waves have been observed in many scientific contexts and—as it is now widely accepted—play an important role in information processing in many biological systems. EM are usually modeled as a collection of locally coupled units with excitable or threshold dynamics [1]. Neighboring units of the medium interact with each other via diffusion-like transport processes which, in many cases, results in a spatially homogeneous steady state. Recent studies extended the notion of excitability of EM to spatially distributed systems, whose unperturbed steady state resembles spatiotemporal chaos. This results in a steady state that actively destroys long-range spatiotemporal correlations and thus prevents the propagation of a single excitation. To generate and then to facilitate a wave-propagation phenomenon, the EM have to be driven with perturbations of appropriate amplitudes and certain (resonant) frequencies. Thus, the phenomenon of frequency-selective excitation waves in EM [2–5] provides a possible way to extend the conventional amplitude-selective notion of excitability in dynamical systems. Applying an unknown signal as a local perturbation to the first unit of EM and measuring the number of excitations in the last unit provide a possibility to detect the presence of rhythmic components in the applied perturbation. Moreover, by increasing the spatial size of the medium it is possible to suppress accidental spikes induced in the first unit and thus to enhance frequency-selectivity of this approach.
The concept of Cellular Neural Networks (CNNs) provides a solid framework to model and to study pattern formation and wave propagation phenomena in many spatially distributed dynamical systems [6–9]. An important subclass of CNN allows one to approximate the behavior of reaction-diffusion partial differential equations (PDEs) and is thus well suited as a framework to model various EM [10]. Probably the most tractable mathematical model of EM is the FitzHugh-Nagumo (FHN) PDE. The discretized version of this equation provides a generic mathematical model of EM and can in principle be implemented in form of reaction-diffusion CNN.
Despite the many advantages, applications of the method of frequency-selective excitation waves in EM are limited. A permanent propagation of excitations in the media requires—by definition—the presence of rhythmic components in the applied perturbation, which limits the applicability of the method to either monochromatic signals or signals with sharp peaks in the power spectrum well above the noise floor. In many field applications, however, one is often confronted with noisy broadband signals.
The human brain is a complex and nonstationary dynamical system. It is capable to generate a vast variety of spatiotemporal patterns that reflect physiologic as well as pathophysiologic states. Electroencephalography is an important tool in neuroscientific research and especially in clinical practice to measure the aforementioned patterns at a high temporal resolution. Electroencephalography is used for diagnostic purposes and in the presurgical evaluation of epilepsy patients. Epilepsy represents one of the most common neurological disorders, second only to stroke. The cardinal symptom of this disease is recurrent seizures that are usually characterized by an abnormal synchronized firing of neurons involved in the epileptic process. But even between seizures, the epileptic brain is different from normal and accompanying alterations seen on the EEG are often associated with complex characteristic spectral patterns and changes in the degree of synchronization between EEG signals. Since these alterations are usually not rigorously defined and, moreover, are patient specific it is rather difficult to detect them automatically in EEG recordings. There is now growing evidence that an improved understanding of the epileptic process can be achieved through the analysis of interactions in epileptic brain networks [11] and further improvements can be expected from time series analysis technique that allow one to estimate the degree of synchronization in the human EEG especially in cases where a high level of noise contaminations cannot be easily avoided.
In this work we apply the method of frequency-selective excitation waves in EM to characterize synchronization phenomena in time series from two coupled chaotic oscillators as well as in multichannel, multiday EEG recordings. For this purpose we relax the property of frequency-selectivity of EM by considering media that are composed of only a few excitable units. With two such media we estimate the coincidence rate of excitations induced in the last units of the media due to a local perturbation of the first units with broadband signals. Our preliminary findings indicate that this approach provides a reliable estimate of the degree of phase synchronization between signals.
2. Methods and Data
2.1. Frequency-Selective Induction of Excitation Waves in Excitable Media
We consider an EM as a diffusively coupled network of resonators, that is, oscillators exhibiting a transition from their resting states to an oscillatory state via a supercritical Hopf-type bifurcation. The FHN oscillator is a generic example of an excitable oscillator that resides (in its parameter space) either near a supercritical or a subcritical Hopf bifurcation [4, 12, 13]. We here model an EM as a chain of only L 3 diffusively coupled FHN oscillators with zero-flux boundary conditions. The equations of motion of the th unit read:
where is a time scaling coefficient, is the Kronecker delta, and denotes an external perturbation. Following [4] the equations were integrated using a fourth-order Runge-Kutta algorithm with step size dt 0.05, followed by a downsampling to dt 0.5. The parameters were set in the Canard region near a supercritical Hopf bifurcation (a 0.1, b 0.015, m 0.015, I 0.062). With this choice and without spatial diffusion (D 0) each unit exhibits low amplitude (subthreshold) oscillations with a characteristic frequency . The time scaling coefficient is used to adjust the characteristic frequency of the model with respect to the external perturbation. A supra-threshold perturbation leads to an excited state, which is defined by oscillations (i.e., spikes) that exhibit a significantly higher amplitude ( 0.9) than those of an unperturbed unit (0.3), as already shown in [4]. A nonzero diffusive coupling (D 0.5) leads to an active desynchronization of all units [14]. Such a spatiotemporally disordered state of the EM prevents wave propagation that is induced by a single short-lasting local supra-threshold perturbation . However, wave propagation is facilitated if the medium is perturbed locally with a frequency located near the characteristic frequency of an unperturbed FHN unit . The medium thus represents a sort of narrow band-pass frequency filter exhibiting frequency-selective firing of the last unit . The periodic firing of the last unit can thus be considered as a detection of rhythmic components in the perturbation signal applied to the first unit ; see (1). To quantify detection of rhythmic components in the local perturbation of a medium we define the spike rate in the dynamics of the last unit :
where is the Heaviside step function , and denotes a threshold for the spike amplitude, which we set to S = 0.9 in order to match the spike amplitude of an uncoupled FHN unit. counts the number of spikes induced in the last unit of the medium within a time window that consists of data points. The normalization factor defines the maximum possible number of spikes (i.e., during permanent firing of the last unit) within the given time window, that is, . Thus, indicates detection of rhythmic components in the local perturbation.
In Figure 1 we show the spike rate in the last unit of a medium as a function of amplitude and frequency of a sine wave that is applied as a local perturbation to the first unit. In contrast to the results obtained in the original works on this method [2–5] we observe a rather smooth response over the frequency range Hz. By increasing the number of units within the medium and/or decreasing the diffusive coupling strength between units it is possible to obtain a rather sharp resonant peak and thus increase the frequency-selectivity of the method [4].
Figure 1: Spike rate in the last unit of a medium composed of L 3 FHN units. The first unit was perturbed with sine waves of varying amplitudes and frequencies. The time scaling coefficient was set to 10 and time series consisted of 5000 data points.
Since we aim at approximating the degree of phase synchronization between two signals (consisting of data points each) we here consider two identical EM, say A and B, each tuned to the same characteristic frequency [15]. According to (1) only the first units (i 1) of the media are perturbed by the signals. Let us assume that the signals contain rhythmic components with frequency . A perturbation of the media will then lead to waves propagating throughout the media. The coincidence between spikes generated by the last units (i L) of media A and B per time step can be used to approximate the degree of phase synchronization between rhythmic components in the applied signals. We define
where again denotes a threshold for the spike amplitude, is the Dirac delta function, and is a coincidence threshold with 1 P 2. This definition implies that by setting the limiting case P 2 the coincidence rate will provide an estimate of zero-lag phase synchronization only. By decreasing the coincidence threshold it is possible, however, to cover the regime of lagged phase synchronization. In the following we set P 1.3 to allow a certain temporal mismatch between spikes in the last units of the media. Preliminary investigations using time series from different model systems indicate that with this setting it is indeed possible to detect a lagged phase synchronization in broad band signals. However, setting too close to 1 increases the risk of detecting accidental coincidences; that is, the situation where accidentally induced spikes will contribute to . This results in a certain bias in the regime of weak synchronization as will be shown later in Figure 3. The normalization factors and are defined as above. The coincidence rate is confined to the interval [] and approaches 1 for identical spike sequences. For large , C(A, B) 0 if the sequences of spikes are completely uncorrelated or if the perturbations (even identical perturbations) do not contain rhythmic components matching the characteristic frequencies.
2.2. Measuring Phase Synchronization
In order to demonstrate the applicability of our method to characterize the degree of synchronization between two signals we compare it with a conventional time series analysis technique. Traditionally, phase synchronization is defined as the locking of the phases of two oscillating systems 1 and 2 [16]:
In order to quantify the degree of phase synchronization we used the mean phase coherence [17, 18]
where is the sampling rate of the discrete time series of length . denotes the circular variance of an angular distribution obtained by transforming the differences in phase onto the unit circle in the complex plane [19]. By definition is confined to the interval [] where R 1(V 0) indicates fully synchronized systems. In order to determine the phases of two signals and , we followed the analytic signal approach [20, 21] which renders an unambiguous definition of the instantaneous phase for an arbitrary signal :
where
is the Hilbert transform of the signal ( denoting the Cauchy principal value). Application of the convolution theorem turns (7) into
where denotes Fourier transform and inverse Fourier transform, respectively.
An important property of the analytic signal approach is that the instantaneous phase always relates to the predominant frequency in the Fourier spectrum [22] of a broadband signal. A frequency-selective analysis of phase synchronization thus either necessitates prefiltering of signals before applying the Hilbert transform or an alternative definition of phases, for example, by the complex wavelet transform [23]. Both definitions of a phase are mathematically related to each other [24]. Quantifying the degree of phase synchronization via the coincidence rate is frequency-selective by definition. Note that in contrast to the coincidence rate the mean phase coherence always approaches 1 for two identical signals.
2.3. Model Systems
In order to test our approach we analyzed time series generated by two diffusively coupled Rössler systems:
with a small mismatch of the natural frequencies and . The parameter denotes the coupling strength. The differential equations were integrated using a fourth-order Runge-Kutta algorithm with step size dt 0.05 [25]. We generated twenty realizations (using different initial conditions, normally distributed with zero mean and unit variance) of time series consisting of N 105 data points. In order to eliminate transients we discarded the first points. The data were then normalized to zero mean and standard deviation and eventually contaminated with additive uniform white noise at a signal-to-noise ratio of SNR −6 dB. The amplitude of the noise is thus twice as high as the amplitude of the signal. In order to measure the degree of synchronization when varying the coupling strength we analyzed the -components of the Rössler systems using the coincidence rate and the mean phase coherence .
2.4. Field Data
In epilepsy monitoring units, EEG is recorded over an extended period of time for diagnostic purposes or for the presurgical evaluation of candidates for resective therapy. Successful surgical treatment of so-called focal epilepsies requires exact localization of the brain region that can generate seizures and its delineation from eloquent cortex that is indispensable for defined cortical functions (see [26] for a comprehensive overview). This task mostly relies on the observation of typical seizures on the video-EEG, which is currently regarded as the gold standard. Epileptic seizures, however, often represent a rather infrequent phenomenon, and thus the question arises as to what extent information obtained from EEG data recorded during the seizure-free interval can help to identify and to delineate the seizure-generating brain region.
We here retrospectively analyzed long-term, multichannel EEG recordings from an epilepsy patient suffering from temporal lobe epilepsy. For this patient the presurgical workup indicated an epileptic focus in the left mesial temporal lobe and after surgical resection of the seizure generating structure the patient is seizure-free. Intracranial EEG were recorded continuously over a longer period from electrodes implanted bilaterally along the longitudinal axis of the hippocampal formation (cf. Figure 2). The data were sampled at using a 16-bit analog-to-digital (A/D) converter and filtered within a frequency band of . The recording lasted 104 hours during which ten seizures occurred. We used a moving-window approach (nonoverlapping segments of N 214 data points) to measure synchronization between all pairs of EEG signals. Since we were interested in recordings from the seizure-free interval only, we excluded data recorded 90 minutes prior to seizures, during actual seizure activities, and from periods of 30 minutes duration immediately following seizures. Prior to estimating the coincidence rate signals were normalized to zero mean and standard deviation .
Figure 2: Schematic of intracranially implanted depth electrodes. (a) Axial view; (b) sagittal view. Each catheter-like, 1-mm-thick silastic electrode contains 10 cylindrical contacts of a nickel-chromium-alloy (2.5 mm) every 4 mm.
Figure 3: Means and standard deviations of the mean phase coherence (circles and triangles) and the coincidence rate (squares) estimated from the noisy (SNR −6 dB) -components of coupled Rössler systems for varying coupling strength . Triangles indicate values of from low-pass filtered data. To calculate the coincidence rate the internal frequency of FHN units in both media was chosen to match the natural frequencies of the Rössler oscillators ( by setting ).
The intrahippocampal EEG does not have a white spectrum but instead has a predominant spectral content in the so-called (0.5–4 Hz) and (4–7 Hz) range [27]. For a characterization of synchronization phenomena in EEG recordings we therefore concentrated on these frequency bands and adjusted the internal frequency of the FHN-units (by setting the time scaling coefficient to (cf. (1) and Figure 1)) such as to cover this frequency range.
3. Results
3.1. Model Systems
In Figure 3 we show the dependence of the coincidence rate and of the mean phase coherence on the coupling strength for the noise-contaminated diffusively coupled Rössler systems. In addition, we show the dependence of for time series that we low-pass filtered (3rd-order Butterworth characteristic; cutoff-frequency set to ) prior to analysis. Despite the strong noise contamination (SNR = −6 dB) allowed to detect the transition to the regime of phase synchronization () between the Rössler systems. Due to the noise contamination the resolvability of the different synchronization regimes using the mean phase coherence was drastically reduced but could be recovered with an appropriate filtering of the data.
Note that even in the regime of no or weak coupling ( we observed C 0.45, which resulted from accidental spikes propagating throughout the media. This bias can be reduced by increasing the spatial size of the media, however, at the expense of limiting the applicability of the method to narrow band signals.
3.2. Field Data
In Figure 4 we show, as an example, for all intra- and interhemispheric channel combinations (cf. Figure 2) the averaged coincidence rate , obtained from averaging -values from the seizure-free periods only. We performed the same analysis with the mean phase coherence and findings are presented in Figure 5. Despite apparent limitations of our method in the regime of weak phase synchronization already observed with model systems, the coincidence rate captures spatial patterns of synchronization similar to those observed with the mean phase coherence . With both approaches we observe that the average degree of synchronization from intrahemispheric channel combinations exceeds the one from interhemispheric channel combinations. Interestingly, even a crude temporal averaging (i.e., from data covering more than four days of recording) and spatial averaging (i.e., from data from intrahemispheric channel combinations only) the degree of synchronization within the left brain hemisphere containing the epileptic focus (ipsilateral hemisphere; ) was slightly higher than in the right (contralateral) brain hemisphere . Analyses using the mean phase coherence revealed qualitatively similar results (cf. Figure 5) with and .
Figure 4: Averaged coincidence rates estimated for all channel combinations of EEG data recorded during the seizure-free intervals. Note the difference in the degree of intrahemispheric synchronization between the focal (left) hemisphere and the nonfocal (right) hemisphere.
Figure 5: Same as Figure
4 but for the averaged mean phase coherence
. EEG data were low-pass filtered with a 3rd-order Butterworth filter and the cutoff-frequency set to 7 Hz.
4. Conclusion
We presented a biologically inspired approach to estimate the degree of phase synchronization in noisy signals. It is based on the method of frequency-selective excitation waves in excitable media whose coincident outputs are used as an estimate for the degree of phase synchronization in locally applied perturbations. By considering EM that are composed of only a few excitable units and thus by relaxing the frequency-selectivity of EM it was indeed possible to apply the method to signals with broadband spectra. From our preliminary investigations using model systems with well-known properties we could show that the method appears sensitive to intermediate and strong couplings even in cases where other approaches (e.g., the mean phase coherence ) require prefiltering of signals. However, our investigations also indicate that the method appears to be insensitive to detect weak couplings due to the propagation of accidental spikes. Moreover, a careful selection of several important parameters—such as the number of excitable units constituting the medium, the diffusive coupling or the coincidence threshold —is required when analyzing broadband signals from dynamical systems with only poorly understood properties.
Despite these limitations, preliminary analyses of synchronization phenomena in multichannel, multiday electroencephalographic recordings from an epilepsy patient indicate that our method provides results that are qualitatively similar to those obtained with the mean phase coherence . In line with findings reported in [17, 27–31] we here observed a higher degree of phase synchronization in brain electrical activity recorded from the hemisphere containing the epileptic focus as compared to the nonaffected hemisphere. Although promising, these findings need to be validated on the data from a larger group of patients. A more detailed study including other frequency bands () on a larger data base is currently underway and will be presented elsewhere.
Extending our previous findings [14, 15] and alternative approaches [32–35] we consider a future implementation of excitable media as analog or digital electrical circuits (e.g., as a FitzHugh-Nagumo cellular neural network) as an attractive perspective for the development of miniaturized detectors of synchronous activities in noisy field data.
Acknowledgments
The authors would like to thank Ronald Tetzlaff and his group for close support. This work was supported by the Deutsche Forschungsgemeinschaft (Grant no. LE 660/2-4).