## Abstract

Quasi-stationarity is ubiquitous in complex dynamical systems. In brain dynamics, there is ample evidence that event-related potentials (ERPs) reflect such quasi-stationary states. In order to detect them from time series, several segmentation techniques have been proposed. In this study, we elaborate a recent approach for detecting quasi-stationary states as recurrence domains by means of recurrence analysis and subsequent symbolization methods. We address two pertinent problems of contemporary recurrence analysis: optimizing the size of recurrence neighbourhoods and identifying symbols from different realizations for sequence alignment. As possible solutions for these problems, we suggest a maximum entropy criterion and a Hausdorff clustering algorithm. The resulting recurrence domains for single-subject ERPs are obtained as partition cells reflecting quasi-stationary brain states.

## 1. Introduction

The brain is an open complex system that receives a broad range of sensational inputs, carries out some computations and triggers various sorts of output. The corresponding neural information is processed in time either sequentially, for example, as in low-level responses to sensory input [1–4], or in parallel in different brain areas as in most higher-level neural processes, for example, in motor action [5] and cognition [6,7]. In both sequential and parallel processes, electrophysiological recordings exhibit a sequence of quasi-stationary signal states that reflect the underlying neural processing steps [8–11]. Here, quasi-stationary states evolve on a large time scale, whereas transitions between them are much faster.

Sequential patterns of quasi-stationary states have been observed in various quite different natural systems. For instance, the chaotic Lorenz attractor exhibits two quasi-stationary states (the so-called *Lorenz wings*), which are visited recurrently over time. The Lorenz model was one of the first nonlinear atmospheric models applied in weather forecasting [12], while it also describes the dynamics in solid-state laser systems [13]. In neural systems, recurrent temporal activity patterns are known from the replay of memory in the hippocampus [14], in bird songs [9] and in electroencephalographic (EEG) data obtained during epileptic seizures [15–17]. Some more detailed analysis of the dynamics of quasi-stationary neural states reveals properties of metastable attractors, which attract and repel in different directions in phase space. Good examples for such metastable attractors are saddle nodes with an attracting and a repelling manifold. Hence, sequences of metastable attractors that are connected along their repelling and attracting directions could generate measurable signals of neural activity exhibiting sequences of quasi-stationary states. Mazor & Laurent [18] found such a sequence of metastable attractors in the activity of projection neurons in the insect olfactory bulb [10,18], whereas Freeman & Rogers [19] revealed sequential rapid phase transitions between spatially synchronized states in EEG signals. Most recently, Hudson *et al.* [20] discovered metastable transition networks in the recovery from anaesthesia.

During cognitive tasks, Lehmann and co-workers [8,21–24] observed segments of quasi-stationary EEG topographies, which they called *brain microstates*. The EEG related to cognitive tasks exhibits a sequence of signal components and is termed *event-related potential* (ERP). Lehmann *et al.* were able to identify the extracted microstates with ERP components, which, in turn, are signal components well established in neuropsychology as markers of certain neural cognitive processes. This link between quasi-stationary signal states and ERP components corroborates the working hypotheses of this study: ERP components reflect quasi-stable attractor states in brain dynamics [25–28]. This hypothesis indicates that the brain performs information processing in a sequence of steps [29] that is reflected in quasi-stationary signal states. Therefore, improving the detection of such states in neurophysiological data promises to improve the study of the corresponding states, their neural origin and essentially their meaning in neural information processing.

For the detection of quasi-stationary states from real-world data, several methods have been proposed. Obviously, one may first plot the signal trajectory in signal space and extract the time windows of apparent quasi-stationarity by visual inspection. To this end, the dimensionality of multivariate signals may be reduced by principal component analysis or related techniques, cf. [18,20]. These methods give first insights into the dynamics, but allow neither classification of the quasi-stationary states nor extraction of their durations and their transition phases. More advanced methods for segmenting quasi-stationary states have to consider (i) non-stationary signal features, i.e. applying temporally instantaneous methods, (ii) multivariate time series, i.e. multidimensional signals exhibiting heterogeneous signal characteristics in different dimensions, and (iii) low signal-to-noise ratio where the noise level is not known *a priori* [30,31].

Techniques applied in the context of brain–computer interfaces aim to extract signal features from multivariate signals [32,33]. Most of them perform a first supervised learning task in order to determine the specific features to be extracted from the EEG. However, this learning phase is not applicable in the detection of quasi-stationary states. This work aims to extract signal features that are unknown *a priori* and which may differ slightly between different trials in an experiment.

One of the first successful methods dealing with these signal characteristics is the microstate segmentation technique of Lehmann *et al.* [8,21,23,24,34], which extracts transient quasi-stationary states in multivariate EEG signals based on the similarity of their spatial scalp distributions. It considers the multivariate EEG as a temporal sequence of spatial activity maps and extracts the time windows of the microstates by computing the temporal difference between successive maps. This procedure allows one to compute the time windows of microstates from the signal and classifies them by the spatial averages over EEG electrodes in the extracted state time window. This method works well for high-dimensional EEG signals with good spatial resolution and high signal-to-noise ratio. Other similar methods extract signal features by clustering techniques [26–28] taking into account the similarity of the time series in different electrodes. Such methods are less sensitive to noise compared with the microstates method because they avoid the computation of temporal differences.

Although methods for experimental single-trial EEG data exist, they imply a missing statistical evaluation of the detected time windows, i.e. one cannot ensure that the states found reflect an underlying neural metastable attractor and not just noise. Hence, to statistically evaluate the gained results, it is necessary to perform a statistical evaluation study over several trials. However, in turn, this implies the need to map detected temporal segments of different trials to each other and, in the best case, align them to each other. This task is difficult to perform, because single trials are known to exhibit temporal jitter to each other [35,36], i.e. they are delayed and shortened to each other. Previous studies have attempted successfully to deal with this problem by applying dynamic time warping on EEG [37]. This technique considers a certain signal template and detects it in the dataset under study with respect to dilation and shrinking operations. The method is not well adjusted to our problem because it is based on a predefined signal template, which is not known in practice.

This study elaborates and further improves our recently developed recurrence domain analysis method [38] in two ways. First, we suggest an alternative method for optimizing the neighbourhood size in the recurrence plot (RP), because the earlier proposal turned out to exhibit some bias. We modify our maximum entropy principle by an additional recoding step that maps transients in the symbolic dynamics onto one distinguished symbol. Second, we suggest a new technique for identifying and aligning recurrence domains obtained from different realizations of a dynamic process (e.g. different trajectories starting from randomly prepared initial conditions or different trials in an EEG experiment). We introduce a Hausdorff partition clustering method achieving this alignment. Section 2 introduces the technique of recurrence domain analysis and applies a new method on single realizations, i.e. single datasets. Then, several extensions of this method are presented followed by a novel method based on the Hausdorff distance to align multiple symbolic sequences gained from multiple datasets. Section 3 concludes.

## 2. Symbolic recurrence domain analysis

Nonlinear dynamical systems generally obey Poincaré's famous *recurrence theorem*^{1} [39], which states that almost all trajectories starting in a ‘ball’ *B*_{ɛ}(*x*_{0}) of radius *ɛ*>0 centred at an initial condition *x*_{0}∈*X* (where denotes the system's phase space of dimension *d*) return infinitely often to *B*_{ɛ}(*x*_{0}) as time elapses. For time-discrete (or discretely sampled) dynamical systems, these recurrences can be visualized by means of Eckmann *et al.*'s [40] RP technique where the element
2.1
of the *recurrence matrix* ** R**=(

*R*

_{ij}) is unity if the state

*x*_{j}at time

*j*belongs to an

*ɛ*-ball centred at state

*x*_{i}at time

*i*, i.e.

*x*_{j}∈

*B*

_{ɛ}(

*x*_{i}), and zero otherwise [40,41].

Often, RPs exhibit a characteristic ‘chequerboard texture’ [40] indicating the system's *recurrence domains*, which are quasi-stable regions in phase space with relatively large dwell times that are connected by transients. Figure 1*c*, for example, depicts a typical RP for a so-called heteroclinic contour connecting saddle nodes (details are reported in §2*a*). Because states moving along such a trajectory slow down in the vicinity of the saddles, these can be regarded as quasi-stationary states. Hence, they appear as ‘blobs’ in the RP figure 1*c*. Every time the system returns into one of these recurrence domains, the RP displays such a texture, which is symmetric with respect to parallel translations along the horizontal and vertical time axes, which is the case for the ‘blobs’ in figure 1*c*. Other paradigmatic examples for recurrence domains are, for example, the Lorenz wings centred around the unstable foci of the Lorenz attractor [12], or saddle sets (such as saddle nodes or saddle tori) that are connected by *stable heteroclinic sequences* (SHS) as shown in figure 1*c* [10,42,43].

In a recent study [38], we have proposed a method for the detection of recurrence domains by means of symbolic dynamics [44–46]. Our approach is motivated by the fact that a recurrence *R*_{ij}=1 leads to overlapping *ɛ*-balls *B*_{ɛ}(*x*_{i})∩*B*_{ɛ}(*x*_{j})≠Ø that can be merged together into equivalence classes, which eventually partition, together with their complements and the still isolated balls, the phase space into its respective recurrence domains. Merging *A*_{j}=*B*_{ɛ}(*x*_{i})∪*B*_{ɛ}(*x*_{j}) into a set *A*_{j} if states *x*_{i} and *x*_{j} are recurrent (*R*_{ij}=1) and if *i*>*j*, we simply replace the larger time index *i* in the RP ** R** by the smaller one

*j*. Let be a finite partition of the phase space

*X*into

*n*disjoint sets

*A*

_{k}. Then, the discrete trajectory (

*x*_{i}) is mapped onto a symbolic sequence

*s*=(

*s*

_{i}) according to the cells of being visited by the states

*x*_{i}, i.e.

*s*

_{i}=

*k*when

*x*_{i}∈

*A*

_{k}.

In the study [38], we suggested creating the symbolic sequence *s*=(*s*_{i}) from an initial sequence of time indices *r*_{i}=*i* subjected to a rewriting grammar [47], which we refer to henceforth as the *recurrence grammar* because this grammar is simply another interpretation of the recurrence matrix ** R** obtained in equation (2.1). It contains

*rewriting rules*that replace large time indices

*i*by smaller ones

*j*(

*i*>

*j*) when the corresponding states

*x*_{i}and

*x*_{j}are recurrent,

*R*

_{ij}=1. Moreover, if three states

*x*_{i},

*x*_{j}and

*x*_{k}are recurrent,

*R*

_{ij}=1 and

*R*

_{ik}=1, for

*i*>

*j*>

*k*, the grammar replaces the two larger indices

*i*,

*j*by the smallest one

*k*. In other words, the recurrence grammar comprises rewriting rules 2.2 and 2.3

We illustrate the action of the recurrence grammar with a simple example. Let us assume that a trajectory consists of only five data points (*x*_{1},*x*_{2},…,*x*_{5}) giving rise to the recurrence matrix
2.4
The algorithm commences with the fifth row, detecting a recurrence *R*_{51}=1. Because 5>1, a rewriting rule 5→1 is generated. Because the next recurrence in row 5 is trivial, the algorithm continues with row 4, where *R*_{42}=*R*_{43}=1. Now, two rules 4→2 and 3→2 are created. The remaining rows 3, 2 and 1 can be neglected owing to the symmetry. Thus, we obtain the following recurrence grammar
2.5
Recursively applying this grammar to the series of time indices *r*=12345 yields *s*=12221, i.e. a system with two recurrence domains 1 and 2. Note that recurrence grammars could be ambiguous, i.e. one left-hand symbol could be expanded into different right-hand symbols by several rules. Yet, this is not a problem for the segmentation of the trajectory into recurrence domains when the grammar is recursively applied. After at most two iterations, all ambiguities become resolved through the smallest time index symbolizing the earliest recurrence domain.

Furthermore, the symbolic alphabet of the rewritten sequence may contain several gaps. We therefore augmented our algorithm by a search-and-replace mechanism recoding unused smaller indices by larger ones, thereby squeezing the symbolic repertoire together.

From the resulting symbolic sequence *s*, a symbolic recurrence plot [48,49] can be computed. Doing this for our example above, we obtain figure 1*d* where the ‘blobs’ have been converted into squares. Now, the chequerboard texture is clearly visible: each square corresponds to a quasi-stationary state captured by a recurrence domain in phase space.

### (a) Numerical simulations

For further illustration and validation of our algorithm, we carry out a numerical simulation study of a stable heteroclinic contour, i.e. a closed SHSs between three saddle nodes [10,42,43] that may be regarded as a model for event-related brain potentials, according to our working hypothesis [26,28].

The stable heteroclinic contour was obtained from a generalized Lotka–Volterra dynamics of *n*=3 competing populations
2.6
with growth rates *σ*_{k}>0, and interaction weights *ρ*_{kj}>0, *ρ*_{kk}=1. The particular parameter settings were *σ*_{1}=1, *σ*_{2}=1.2 and *σ*_{3}=1.6.

Lotka–Volterra systems are well studied in theoretical ecology where they describe the intra- and interspecific interactions, i.e. competition for limited resources and predator–prey networks, among *n* populations with abundances *x*_{k} in an ecosystem. Moreover, these systems have a long-standing tradition in computational neuroscience [50,51] as a suitable approximation [52] of the Wilson–Cowan neural mass model [53].

The Lotka–Volterra equations often describe limit cycles in appropriate parameter regimes. However, they also exhibit SHS [10,42,43] and stable heteroclinic channels (SHCs) [10,54] for asymmetric interactions, which is of particular importance for neural applications.

#### (i) Recurrence grammars

The heteroclinic contour in our simulation is governed by: *ρ*_{3,1}=*σ*_{3}/*σ*_{1}+0.5, *ρ*_{2,1}=*σ*_{2}/*σ*_{1}−0.5 for saddle 1; *ρ*_{1,2}=*σ*_{1}/*σ*_{2}+0.5, *ρ*_{3,2}=*σ*_{3}/*σ*_{2}−0.5 for saddle 2; and *ρ*_{2,3}=*σ*_{2}/*σ*_{3}+0.5, *ρ*_{1,3}=*σ*_{1}/*σ*_{3}−0.5 for saddle 3. The system was initialized with ** x**=(1,0.17,0.01)

^{T}and numerically integrated over [0,100] with

*Δt*=0.1429 and the Runge–Kutta algorithm

`ode45`of MATLAB. Figure 1 shows the resulting trajectory

**(**

*x**t*) in figure 1

*a*and the three coordinates,

*x*

_{1}(

*t*),

*x*

_{2}(

*t*) and

*x*

_{3}(

*t*) in the upper panel of figure 1

*b*.

Regarding the coordinates as ‘recording channels’ of an (idealized) EEG experiment, we can identify the traces of *x*_{1} (blue online), *x*_{2} (yellow online) and *x*_{3} (green online) as ‘ERP components’, peaking at the respective saddle nodes in phase space.

For illustration purposes only, we show a typical RP of the dynamics obtained for *ɛ*=0.15 in figure 1*c*. The RP exhibits a ‘skeleton’ of parallel diagonal lines which is characteristic for oscillatory dynamics: the distance of two lines reflects the oscillation's period. As already mentioned above, the saddles appear as ‘blobs’ owing to the system's slowing down around these quasi-stationary states. Applying our recurrence domain algorithm to the simulated data with a different *ɛ**=0.017 yields the symbolic segmentation displayed in the bottom panel of figure 1*b*. Each peak is now assigned to one recurrence domain symbolized by the same colour. The corresponding phase space partition is depicted by the coloured balls in figure 1*a*. Moreover, figure 1*d* depicts the symbolic RP resulting from this encoding, where the ‘blobs’ eventually became squares, thus indicating the partitioning of the system's phase space into equivalence classes.

One of the pertinent problems in recurrence analysis is optimizing the ball size *ɛ* [55]. In our previous study [38], we suggested maximizing the ratio of symbol entropy
2.7
over the cardinality of the symbolic repertoire *M*(*ɛ*), where *p*_{k} is the relative frequency of symbol *k* in the symbolic sequence *s*. However, this utility function strongly penalizes large alphabets through the denominator *M*(*ɛ*) such that the proposed procedure is biased towards binary partitions, which is not appropriate in most applications, such as in our SHS example as well. Therefore, here we remedy this bias by an additional recoding of *s*. The symbolic dynamics obtained from recurrence grammar encoding still contains subsequences (or ‘words’) of monotonically increasing time indices, resulting from isolated *ɛ*-balls. These isolated balls, reflecting the transients of the dynamics, entail large symbol alphabets that are punished by the utility function above.

For our new optimization method, we scan the sequence *s* for monotonically increasing time indices, such as 4567, for example, which are then replaced by a new symbol 0, resulting in the subsequence 0000. Thereby, all isolated balls receive the same label 0, now indicating ‘transient’ in the symbolic dynamics. This new symbol is plotted in red (online) in figure 1; red segments and partition cells therefore denote transients in our encoding.

Hence, the optimization algorithm performs a sequence of three encoding steps: (i) rewriting by the recurrence grammar, (ii) compressing the alphabet for closing gaps, and (iii) mapping transients onto 0. Figure 2 shows the intermediate results of this optimization procedure. In figure 2*a,* we plot the final symbolic sequences (called *s*, again) as a function of the ball size *ɛ*.

The symbolic dynamics plotted in figure 2*a* starts with *ɛ*=0.001 in the first row. Here, almost all balls remain isolated, which is shown by the symbol 0 printed in olive here. For increasing *ɛ*, the number of recurrences increases owing to the trajectory's slowing down in the vicinity of the saddle nodes. These saddles correspond to the ‘tongues’ widening up to the recurrence domains. For some range of *ɛ*, the results are quite robust, indicating the stability of our algorithm. At the bottom row of figure 2*a,* *ɛ*=0.032, eventually all balls merge together into one domain covering the entire trajectory. This trivial domain is printed in dark brown here.

Computing the Shannon entropy equation (2.7) for each row in figure 2*a* yields the distribution in figure 2*b*. For small and large *ɛ*, the symbolic dynamics is highly degenerate: for small values, it consists mainly of the transient symbol 0; for large values, there is essentially only one remaining symbol. In both cases, these symbols are highly predictable, which is reflected by low entropies. Because *p*_{k} is the relative frequency of symbol *k* in the symbolic sequence *s*, it can be regarded as the system's dwell time in partition cell *A*_{k} and therefore as the duration of either a quasi-stationary state or of a transient regime. Considering now an approximately uniform symbol distribution with maximal entropy, we arrive at an optimization constraint maximizing both the durations of quasi-stationary states and their separation through transient regimes of comparable duration. In our particular example of figure 2*b,* the entropy distribution is relatively broad and exhibits several local maxima. Therefore, we suggest looking not for one of these (probably) spurious maxima for optimizing *ɛ*, but rather for the median of a symmetric distribution, or likewise for another quantile in the case of an asymmetric entropy distribution. In this way, we obtained the optimal ball size *ɛ** for figure 1 through the median (the *q*=0.5 quantile).

#### (ii) Hausdorff partition clustering

Lotka–Volterra systems as in equation (2.6) do not only possess SHS. These trajectories are additionally captured by SHCs [54] that illustrate another important property of event-related brain potentials. Adopting our working hypothesis again, that ERP components correspond to saddle sets, or, more generally, to recurrence domains in the brain's phase space [26,28], we need to understand the emergence of ERP waveforms from ensemble averages over realizations of the system's dynamics starting from randomly distributed initial conditions [3,56].

To this end, we simulated an ensemble of 20 trajectories starting from randomly distributed initial conditions in the vicinity of the first saddle node (the blue recurrence domain in figure 1). Figure 3 shows the phase portrait of this ensemble (each realization plotted in a different colour) in figure 3*a* and the *x*_{3}-components in figure 3*b*. In addition, we present the ensemble averages as the fat black curves in both.

Figure 3 nicely demonstrates that all simulated trajectories are confined to a ‘tube’, namely the SHC, stretched across the saddle nodes. This channel is stable against small perturbations in the initial conditions and also against additive noise [54], i.e. for moderate noise levels, the dynamically evolving states never leave the SHC. Moreover, figure 3 shows the dispersion of the ensemble average which is inevitable in standard ERP analyses (yet see [36,56] for alternative approaches). This dispersion is due to the velocity differences dependent on the distance of the actual state from the saddles. The closer a state comes to the saddle, the larger is its acceleration along the unstable and its deceleration along the stable manifold. Hence, the velocities with which states explore their available phase space crucially depend on the initial conditions. As the system evolves, the dynamics accumulates phase differences, deteriorating the averaged ensemble trajectory. This finding is consistent with the phase-resetting hypothesis of ERPs [57] as phase-resetting takes place in the preparation of initial conditions.

Our example is also relevant for the larger scope of nonlinear data analysis problems. Computing RPs and subsequent recurrence quantification analysis [55,58] is numerically rather expensive and therefore restricted to relatively short time series. In practice, a long-lasting time series is therefore cut into many segments for which recurrence analyses are carried out individually. The pertinent problem is then comparison and alignment of the individual results.

In the framework of our symbolic recurrence domain analysis, we suggest here a solution based on a clustering algorithm whose results are presented in figure 4. In figure 4*a,* we plot the symbolic sequences of the 20 SHS realizations from figure 3, now encoded with the *q*=0.7 quantile.

Figure 4*a* shows different realizations using essentially the same colour palette that is due to the recurrence grammar and recoding algorithm. However, the same symbol (i.e. colour) does not necessarily refer to the same recurrence domain in the system's phase space (with one exception: transients are always encoded as 0, here plotted in dark blue). Therefore, we have to render the sequences in figure 4*a* comparable with respect to corresponding regions in phase space.

For that aim, we first gather all discrete sampling points belonging to the same phase space cell *A*_{k} that is labelled by the symbol *k*, i.e. *A*_{k}={*x*_{i}|*s*_{i}=*k*}⊂*X*. Because high-dimensional phase spaces are usually very sparsely populated with sampling points from experimental or simulated data, we then project the cells onto the *d*-dimensional hypersphere through *y*_{i}=*x*_{i}/∥*x*_{i}∥. This is also justified by the goal to clustering recurrent topographies of ERP data [8,22,34].^{2}

In the next step, the recurrence domain partitions and of two subsequent realizations are merged together into a set system by regarding all symbols (apart from 0) as being different. This set system is, in general, not a partition any more. From the members of we then calculate the pairwise Hausdorff distances^{3}
2.8
where
2.9
measures the ‘distance’ of the point ** x** from the compact set

*A*⊂

*X*. The Hausdorff distance between two compact sets vanishes when they are intersecting. The recurrence domains taken into account here are clearly compact sets as they are finite unions of intersecting

*ɛ*-balls.

From the pairwise Hausdorff distances (2.8), we compute the *θ*-similarity matrix ** S** with elements
2.10
which has essentially the same properties as the recurrence matrix

**from equation (2.1). Therefore, we merge the members of into new partition cells by interpreting**

*R***as another recurrence grammar for rewriting large indices of**

*S**A*

_{i}by smaller ones from

*A*

_{j}(

*i*>

*j*) when they are

*θ*-similar (

*S*

_{ij}=1).

For numerical implementation, we recursively decompose the ensemble of realizations into halves until only one or two trajectories are obtained. In the latter case, their recurrence domain partitions are merged together and subjected to the recurrence grammar given by the Hausdorff similarity matrix ** S**. The resulting partition contains the unions of similar recurrence domains that are passed to the next iteration of the algorithm.

The result of the Hausdorff partition clustering algorithm applied to our SHC simulation is shown in figure 4*b* for *θ*=0.7. Now, all symbols (colours) refer to the same recurrence domains in phase space (respectively to their topographic projections). Realizations 1 to 12 exhibit very similar behaviours, as they explore the same recurrence domains. Note also the differences in phase and duration in the individual realizations. Although the algorithm shows nice convergence, it is numerically very demanding using the MATLAB interpreter language. However, we have been able to substantially improve and speed up the Hausdorff clustering algorithm. The results will be published in a subsequent publication. Clustering the simulated recurrence domains took approximately 100 h computing time on a conventional PC. Considered as a model of ERP data, figure 4*b* shows great resemblance to statically encoded ERPs [56] where ERP components correspond to meandering vertical stripes.

### (b) Event-related brain potentials

The analysis of electroencephalographic data by means of symbolic dynamics techniques has a long-standing tradition. It could be traced back to Lehmann [60] and Callaway *et al.* [61], who encoded extremal EEG voltages in a binary fashion and estimated the relative frequencies of these symbols across trials. The former approach led later to the so-called half-wave encoding [62], whereas the latter was used through cylinder measures and word statistics [56]. Further developments in this field were the symbolic resonance analysis [63,64] and order pattern analyses [65–67]. In addition, the segmentation of the EEG into quasi-stable brain microstates can be subsumed under these techniques [8,22,34], which were recently combined with spectral clustering methods on approximate Markov chains [17,68,69].

In this section, we reanalyse an ERP experiment on the processing of ungrammaticalities in German [70] (see [71–74] for other studies on symbolic dynamics of language-related brain potentials). Frisch *et al.* [70] examined processing differences for different violations of lexical and grammatical rules. Here, we focus only on the contrast between a so-called *phrase structure violation* (11-b), indicated by the asterisk, in comparison with grammatical control sentences (11-a).

Sentences of type (11-b) are ungrammatical in German, because the preposition *am* is followed by a past participle instead of a noun. A correct continuation could be, for example, *im Garten wurde am Zaun gearbeitet* (‘work was going on in the garden at the fence’) with

*am Zaun*(‘at the fence’) comprising an admissible prepositional phrase.

The ERP study was conducted with 17 subjects in a visual word-by-word presentation paradigm. Subjects were presented with 40 structurally identical examples per condition (i.e. 40 trials). The critical word in all conditions was the past participle printed in bold font above. EEG and additional electrooculogram (EOG) for controlling eye movement were recorded with a 64-electrode montage from which 59 channels were measuring EEG proper. These were selected for recurrence domain analysis.

#### (i) Grand averages

First, we carry out the recurrence domain analysis for the grand averages over all 17 subjects, which are presented in the upper panels of figure 5. At the left-hand side, the multivariate ERP time series for the correct condition (11-a) are plotted as coloured traces, one for each EEG electrode.^{4} The same for the time series of the phrase structure violation condition (11-b) at the right-hand side. Both plots start 200 ms before stimulus presentation at *t*=0. In both conditions, similar N100/P200 ERP complexes are evoked. Moreover, condition (11-b) exhibits a large positive P600 ERP at central and posterior electrodes around 600 ms after stimulus presentation.

For the recurrence domain analysis, we regard the observation space [75] spanned by the instantaneous EEG voltages ** v**=(

*v*

_{i}), with 1≤

*i*≤

*d*as EEG electrode index and

*d*the number of recording channels, as the system's

*d*-dimensional phase space [76,77] and compute entropy distributions

*H*(

*ɛ*) (2.7) in the range

*ɛ*∈[0.2 μV,5 μV] that are very asymmetric. Therefore, we chose the

*q*=0.77 quantile for optimization by visual inspection of both segmentations, yielding

*ɛ**=2.2 μV for the control condition. We then compute the recurrence grammars and the resulting segmentation into

*d*-dimensional recurrence domains with this parameter for both conditions. The results are shown in the bottom panels of figure 5.

From figure 5 (bottom panel), we draw three important conclusions. (i) The pre-stimulus interval and also the first 100 ms are assigned to the same recurrence domain, reflecting some resting state brain activity. (ii) In the time window from 100 to 300 ms, a similar partitioning into successive recurrence domains is observed for both conditions. This indicates similarities in the early attentional processes assigned to the N100/P200 ERP complexes [4]. (iii) Both conditions differ crucially after 300 ms. In the correct control condition (11-a), there is only one large recurrence domain in this time window, whereas the phrase structure violation condition (11-b) exhibits two recurrence domains: a first one in the time window of lexical access and syntactic diagnoses processes around 400–700 ms, and another one in the time window of syntactic repair and reanalysis processes around 700–1200 ms, which is consistent with other findings about subcomponents of the late positivity [73].

#### (ii) Single subject clustering

Finally, we carry out the Hausdorff partition cluster analysis for the phrase structure violation condition (11-b). Here, we present only a proof of concept using the single-subject ERPs shown in figure 6 for one characteristic electrode Pz that is located above the brain's parietal lobe at the midline of the head. Now, each trace depicts the ERP average of one subject.

Single-subject ERPs as shown in figure 6 exhibit much more variation than grand averages from figure 5, which is reflected by the amplitude scales: 40 μV in figure 6 versus 20 μV in figure 5. Furthermore, large differences in the P600 waveforms around 800 ms are also visible between subjects. Thus, single-subject ERPs provide a good benchmark for ensembles of realizations of noisy dynamics.

All individual single-subject ERPs of condition (11-b) are subjected to the same recurrence analysis with the *q*=0.77 quantile of the entropy distribution (2.7). The results are shown in figure 7*a*.

Figure 7*a* displays large differences between individual single-subject ERPs, although some resemblances with the late grand average recurrence domains from figure 5 (bottom right) are present, for example, in subjects 1, 5, 10, 13 and 17. As every symbol must be treated differently across realizations, the cardinality of the alphabet for this representation is about 272.

In order to obtain a good clustering despite the apparent differences, we chose a similarity threshold *θ*=0.985 for our Hausdorff algorithm. The resulting clustering of recurrence domain topographies is depicted in figure 7*b*. The cardinality of the alphabet is drastically reduced to 54 different symbols. Obviously, some recurrence domains are shared by many realizations, as reflected by the repeating symbol appearing in medium blue in subjects 1–3, 7–9, as well as 11, 12 and 16 around 400–800 ms. However, the desired ‘meandering vertical stripes’ are not yet visible.

For further evaluation, the symbolic dynamics obtained from different experimental conditions, different subjects and eventually from single trials could be subjected to statistical analysis by computing word distributions, cylinder entropies and statistical hypothesis tests (e.g. *χ*^{2} tests on word distributions or permutation tests on symbolic ensembles) [56,64]. We leave this as well as further parameter search for optimizing encodings for a future publication.

## 3. Conclusion

Starting from the working hypothesis that event-related brain potentials (ERPs) reflect quasi-stationary states in brain dynamics, we have elaborated an earlier proposal for the detection of quasi-stationary states and for the segmentation of electroencephalographic time series into recurrence domains. This segmentation technique interprets the recurrence matrix as a rewriting grammar that applies to the time indices of an ERP dataset. After further recoding steps, the quasi-stationary states are detected as recurrence domains in phase space as indicated by a symbolic dynamics. Moreover, we have suggested a method for the alignment and unification of multiple EEG trials (i.e. single-subject ERPs) by means of Hausdorff partition clustering.

We think that the presented methods could be of significance for the greater community of nonlinear time series analysis, as we have addressed some pertinent problems in recurrence analysis and symbolic dynamics. We have suggested an optimization procedure for choosing the ball size of RPs by maximizing the entropy of the symbol distribution that is obtained by applying recurrence grammars and subsequent recoding of transients. In addition, the alignment and comparison of RPs of different time windows or different realizations was an unsolved problem so far. We solved this problem by merging together recurrence domains from different time series upon their similarity with respect to the Hausdorff distance either in phase space or in projection onto the unit sphere.

## Funding statement

A.H. and P.b.G. acknowledge funding from the European Research Council for support under the European Union's Seventh Framework Programme (FP7/2007-2013) ERC grant agreement no. 257253. In addition, P.b.G acknowledges support by a Heisenberg Fellowship of the German Research Foundation DFG (GR 3711/1-2).

## Acknowledgements

We thank the guest editors for their kind invitation to contribute to this Theme Issue of the *Philosophical Transactions A*. In this study, we reanalysed a language processing EEG dataset by courtesy of Stefan Frisch, Anja Hahne and Angela Friederici.

## Footnotes

One contribution of 12 to a theme issue ‘Enhancing dynamical signatures of complex systems through symbolic computation’.

↵1 Under the assumptions that the dynamics possesses an invariant measure and is restricted to a finite portion of phase space.

↵2 A classification method based on the estimation of phase space probability densities was reported by Wright & Schult [59]. However, because probability density functions are even harder to estimate for sparsely sampled data, we use the recurrence domain partitioning instead.

↵3 Note that Wackermann

*et al.*[8] used a very similar method for clustering brain microstate topographies.↵4 For a legend of electrode labels, see the electronic supplementary material.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.