Let us consider a set of discretely sampled trajectories {(xt(c))1≤t≤T|t,T,c∈ℕ} in the phase space X of a dynamical system with sampling time t. The system depends on a (discretized) control parameter c ∈ ℕ (hence taken from natural numbers, here), indicating several experimental conditions. T denotes the number of samples and X is a subset of ℝⁿ. An example might be a multivariate measured signal with dimension n. Note that for univariate time series, X can be reconstructed through phase space embedding [28] or, more recently, through spectral embedding techniques [29, 30]. For a given c ∈ ℕ the system's realization (xt(c)) is a function from the index set T = {t|1 ≤ t ≤ T} into X, i.e., (xt(c)) ∈ X^T.

Starting point of the (symbolic) recurrence structure analysis (RSA) is Poincaré's famous recurrence theorem [31] stating that almost all trajectories starting in a “ball” B_ε(x₀) of radius ε > 0 centered at an initial condition x₀ ∈ X return infinitely often to B_ε(x₀) as time elapses, when the dynamics possesses an invariant measure and is restricted to a finite portion of phase space. For time-discrete dynamical systems, these recurrences can be visualized by means of Eckmann et al.'s [32] recurrence plot (RP) technique where the element

of the recurrence matrix R = (R_ij) is unity if the state x_j at time j belongs to an ε-ball

centered at state x_i at time i and zero otherwise [32, 33]. Here d: X × X → ℝ⁺ denotes some distance function, which is a metric in general.

Recurrent events (R_ij = 1) of the dynamics lead to intersecting ε-balls B_ε(x_i) ∩ B_ε(x_j) ≠ ∅ which can be merged together into equivalence classes of phase space X. By merging balls together into a set A_j = B_ε(x_i) ∪ B_ε(x_j) when states x_i and x_j are recurrent and when i > j, we simply replace the larger time index i in the recurrence plot R by the smaller one j, symbolized as a rewriting rule i → j of a recurrence grammar [12, 26, 27].

Applying the recurrence grammar that corresponds to a recurrence plot R for given ball size ε recursively to the sequence of sampling times r_i = i with 1 ≤ i ≤ T yields a segmentation of the system's multivariate time series into discrete states s_i = k, i.e., a symbolic dynamics based on a finite partition P={Ak⊂X|0≤k≤ma} of the phase space X into m_a + 1 disjoint sets A_k, such that x_i ∈ A_k. In this representation the distinguished state 0 captures all transients, while m_a is the number of detected MS [12].

The results of the RSA and the obtained symbolic dynamics depend heavily on the parameter ε that defines the ball size of recurrent events. For finding an optimal partition and hence an optimal value for ε, beim Graben and Hutt [12, 26] and beim Graben et al. [27] have suggested several entropy-based criteria. To gain an optimal segmentation, the idea is to assume an underlying stochastic model for the symbolic sequences. Then the optimal segmentation is as close as possible to the model under consideration. For instance, assuming a uniform probability distribution of symbols, the optimal value ε maximizes the entropy of the states [26]. Here, we employ the most recently suggested Markov optimization criterion [27], maximizing a utility function

for a Markov transition matrix P estimated from the bi-gram distribution of the symbolic sequence s. Here

are renormalized entropies with p0j′=p0j/∑j=1n-1p0j for the first row and pi0′=pi0/∑i=1n-1pi0 for the first column of P. The transition rates p_ij in P are estimated for each value of ε. The parameter ε^* maximizing u(ε) is then the optimal threshold for which the recurrence structure detected from the time series mostly resembles an unbiased Markov chain model.

Applying the RSA to the time series (xi(c)) ∈ X^T for several experimental conditions separately yields a symbolic sequence s^(c) for each condition. Each sequence is based on a unique phase space partition P(c)={Ak(c)⊂X|0≤k≤m(c)}, c ∈ ℕ. In order to unify such different descriptions for all experimental conditions within a common picture based on a single partition, beim Graben and Hutt [12] suggested a recursive Hausdorff clustering method. As this method turned out to be rather time consuming in practical applications, we present a substantially improved non-recursive algorithm for the alignment of multiple realizations of dynamical system's trajectories in the sequel.

For this aim let us assume two experimental conditions, i.e., c = 1, 2, without loss of generality. Our new, more parsimonious, approach is introduced by a recoding of the symbolic sequences s^(c) for c > 1 based on the symbolic sequence s⁽¹⁾. Thus we leave the symbols for the first condition (c = 1) unchanged but alter the symbols for the second condition through an additive constant, i.e.,

with T as the length of each sequence s^(c). This is possible as symbols are simply represented as integer numbers in our framework¹. Note that the “transient symbol” 0 is not affected by the recoding, thereby merging together all transients across conditions.

Afterwards, both the original time series (xi(c)) ∈ X^T and their recoded segmentations q^(c) are concatenated into two long series

such that the concatenation products are now functions from 2T = {t|1 ≤ t ≤ 2T} into X, i.e., (ξ_i) ∈ X^2T.

From these data we gather all sampling points that belong to the same phase space cell B_k labeled by the symbol k, i.e., B_k = {ξ_i|η_i = k} ⊂ X. The family of these m_b sets is in general not a partition, but certainly a covering Q={Bk⊂X|0≤k≤mb} of X. From the members of Q we calculate the pairwise Hausdorff distances [34]

where

measures the “distance” of the point x from the compact set A ⊂ X. The Hausdorff distance of two overlapping compact sets vanishes.

Now we proceed as in the previous approach by beim Graben and Hutt [12]: From the pairwise Hausdorff distances (7) we compute a θ-similarity matrix S with elements

which consists of zeros and ones as the recurrence matrix R from Equation (1). Therefore, S can also be regarded as a recurrence grammar, for merging the members of Q into new partition cells by rewriting large indices of B_i through smaller ones from B_j (i > j) when they are θ-similar (S_ij = 1). The result of the Hausdorff clustering is a unique new partition R={Ck⊂X|0≤k≤mc} entailing a symbolic dynamics p^(c) for both conditions based on the same symbolic repertoire.

For the application to the ERP data shown below, we chose a rather small similarity distance of θ = 0.25 in order to identify at least the pre-stimulus domain across conditions.

In the view of our subsequent modeling in Section 2.2, recall that the found segments of the RSA are MS in phase space. Therefore, we compute their centers of gravity as time-averaged topographies

where x_i ∈ C_k and T_k is the number of samples in C_k.

The recurrence structure analysis, Section 2.1, of the measured data yields a sequential dynamics for each experimental condition. From our present analysis Section 3.1 let us assume that we obtain a total number of m_c MS after Hausdorff clustering. For illustration reasons (cf. Section 3), we assume m_c = 8 and we gain the following pattern sequences

forming n × l matrices for the n-dimensional EEG observation space with l = 6 MS per condition. The spatial patterns result from an underlying system dynamics, that exhibits a HO. To gain deeper insight into the system dynamics that generates the SHS gained experimentally, a second step aims to derive the underlying system model dynamics. To this end, we identify the MS as equilibria in the underlying dynamical system model. Each pattern, i.e., column of V^(c) can be seen as a spatial discretization of a continuous spatial neural field v(x) [35] serving as a saddle in a heteroclinic sequence of field activations [36–38].

In the following we use the findings of beim Graben and Hutt [36] and Schwappach et al. [38] for the modeling of the SHS dynamics by means of heterogeneous neural fields. Our starting point is the Amari equation

describing the evolution of neural activity u(x, t) at site x ∈ Ω ⊂ ℝ^b and time t [39]. Here, Ω is a b-dimensional manifold, representing neural tissue. Furthermore, w(x, y) is the synaptic weight kernel, and f is a sigmoidal activation function, often taken as f(u) = 1/(1 + exp(− β (u − θ))), with gain β > 0, and threshold θ > 0. The characteristic time constant τ will be deliberately absorbed by the kernel w(x, y) in the sequel.

The metastable segmentation patterns of the RSA (11) are interpreted as stationary states, v(x), of the Amari equation that are connected along their stable and unstable directions, thereby forming a stable heteroclinic sequence (SHS: 3, 4). Such SHS are examples of heteroclinic orbits connecting different equilibria. Let {v_k(x)}, 1 ≤ k ≤ n be a collection of MS which we assume to be linearly independent. Then, this collection possesses a biorthogonal system of adjoints {vk+(x)} obeying

For the particular case of Lotka-Volterra neural populations, described by activities ξ_k(t),

with growth rates σ_k > 0, interaction weights ρ_kj > 0 and ρ_kk = 1 that are tuned according to the algorithm of Afraimovich et al. [3] and Rabinovich et al. [4], the population amplitude

recruits its corresponding MS v_k(x), leading to an order parameter expansion

of the neural field.

Under these assumptions, beim Graben and Potthast [40] and beim Graben and Hutt [36] have explicitly constructed the kernel w(x, y) through a power series expansion of the right-hand-side of the Amari Equation (12),

with Pincherle-Goursat [41] kernels

The kernel w₁(x, y) describes a Hebbian synapse between sites y and x whereas the three-point kernel w₂(x, y, z) further generalizes Hebbian learning to interactions between three sites x, y, z of neural tissue.

The kernels w₁(x, y) and w₂(x, y, z) allow the identification of the given system, characterized by the real kernel w(x, y) and the activation function f(u) through a Taylor series expansion of f around u = 0:

The linear term in this expansion is simply f′(0)w(x, y) which equals w₁(x, y) above.

Applying this method for each condition separately, yields two kernels w^(c)(x, y) for c = 1, 2. Interestingly, their convex combination

for μ ∈ [0, 1] entails a dynamical system depending on a continuous control parameter μ with limiting cases w^(c)(x, y) for μ = 1, 0.

Summarizing the modeling approach, we assume a Lotka-Volterra dynamics of the underlying system while identifying its fixed points with the MS gained experimentally from the data. To derive a spatio-temporal neural model that evolves according to the Lotka-Volterra dynamics, and hence exhibits the right number of HOs, we have employed a neural field model. The final underlying neural model (17) evolves similarly to the experimental data.

In our later brain signal simulation, we construct a b = 2 dimensional neural field with l = 6 MS per condition that are discretized with a spatial grid of n = 59 sites. Adjoint patterns are obtained as Moore-Penrose pseudoinverses [42, 43] of the matrices (11), i.e.,

For the temporal dynamics we prepare the HO solving Equation (14) with growth rates σ1(1)=0.15,σ2(1)=0.21,σ3(1)=0.27,σ4(1)=0.33,σ5(1)=0.39,σ6(1)=0.45 for condition (22-a) and σ1(2)=0.15,σ2(2)=0.21,σ7(2)=0.27,σ3(2)=0.33,σ6(2)=0.39,σ8(2)=0.45 for condition (22-b), respectively. The interaction matrices have been tuned according to the algorithm of Afraimovich et al. [3] and Rabinovich et al. [4] with a competition bias of ρ₀ = 3.

For the simulations with the neural field toolbox of Schwappach et al. [38], we prepare as initial conditions the first MS v₁ for two control parameters c and integrate the Amari Equation (12) over the pre-stimulus interval [−200, 0] ms. Since v₁ is an (unstable) fixed point, the system remains in this stationary state during this time. At time t = 0, we introduce a slight perturbation of the system mimicking the stimulation by the upcoming word. This kicks the state out of equilibrium and triggers the evolution along the prescribed heteroclinic sequence.

Running simulations for different conditions with their respective parameters allows finally the computation of simulated brain signal data and comparison of their functional connectivities.

In order to validate our neural field model of HO, we subject the simulated data to the recurrence structure analysis as well. Since the amplitude of the simulated data is slightly diminished, we use a Hausdorff similarity threshold θ = 0.1 in case of the experimental data. Other analysis details are the same as above in Section 2.1.

We reanalyze an ERP experiment on the processing of ungrammaticalities in German [44, Exp. 1] for easy comparison with our presentation in beim Graben and Hutt [12]. Frisch et al. [44] reported processing differences for different violations of lexical and grammatical rules. Here we focus on the contrast between a so-called phrase structure violation (22-b), indicated by the asterisk, in comparison to grammatical control sentences (22-a).

In German, sentences of type (22-b) are ungrammatical because the preposition am is followed by a past participle instead of a noun. A correct continuation would be, e.g., im Garten wurde am Zaun gearbeitet (“work at the fence was going on in the garden”) with am Zaun (“at the fence”) as an admissible prepositional phrase.

The ERP study was carried out in a visual word-by-word presentation paradigm with 17 subjects. Subjects were presented with 40 trials per condition, each trial comprising one sentence example that was structurally identical to either (22-a) or (22-b). The critical word was the past participle printed in bold font across all conditions. EEG and additional electro-ocologram (EOG) for controlling eye-movement were recorded with 64 electrodes; EEG was measured with n = 59 channels which spanned the observation space [45] of our multivariate analysis.

For preprocessing, continuous EEG data were cut into [−200, +1, 000]ms epochs, baseline aligned along the prestimulus interval [−200, 0]ms and averaged over trials per condition per subject after artifact rejection. Subsequently, those single-subject ERP averages were averaged over all 17 subjects per condition in order to obtain the grand average ERPs as the bases of our recurrence structure analysis Section 2.1 and neural field modeling Section 2.2.

First we report the grand average ERPs as snapshot sequences of the respective scalp topographies. Figure 1 shows the results for the correct condition (22-a) where every snapshot presents the EEG activation pattern at the indicated point in time. The color scale range is [−15, 15]μV.

The sequence exhibits only one prominent pattern, a frontal positivity around 232 ms, known as the attentional P200 component.

Accordingly, we present in Figure 2 the grand average for the phrase structure violation condition (22-b). Again the color scale ranges between [−15, 15]μV.

Also in this condition the P200 effect is visible as a frontally pronounced positivity between 232 and 280 ms. Moreover, an earlier effect with reversed polarity can be recognized around 136 ms, the tentative N100 component. Most important is a parietal positivity starting at 472 ms until the end of the epoch window at 952 ms. This late positivity is commonly regarded as a neurophysiological correlate of syntactic violations [46], reanalysis [47, 48], ambiguity [49] and general integration problems [50].

However, a proper interpretation of ERP effects relies on considering condition differences. Thus, we plot the ERP difference between violation condition (22-b) and correct condition (22-a) in the range [−8, 8]μV in Figure 3.

The difference patterns in Figure 3 clearly indicate the posteriorly distributed P600 ERP component that evolves between 500 and 1,000 ms post-stimulus.

During cognitive tasks, Lehmann et al. [51, 53] and Wackermann et al. [52] observed segments of quasistationary EEG topographies, which they called brain microstates. The microstate analysis extracts MS 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 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.

In order to detect brain microstates by means of the RSA method in Section 2.1 we compute recurrence plots based on the cosine distance function

as we are interested in detecting recurrent scalp topographies. This choice has also the advantage that the sparse 59-dimensional observation space is projected onto the unit sphere, resulting into a denser representation of ERP trajectories.

For the recurrence structure analysis we optimize the Markov utility function (3) for the grand averages of both conditions (22-a) and (22-b) separately. The resulting utility functions are depicted in Figure 4.

Figure 4 indicates that both conditions lead to optimal segmentations for ε^* ≈ 0.014.

The results of the recurrence structure analysis (RSA) for this optimal value of ε are shown in Figure 5. Figure 5A displays the grand average ERPs for the correct condition (22-a) and Figure 5B for the violation condition (22-b) where each colored trace denotes one recording electrode. The language-related P600 ERP component is clearly visible in Figure 5B as a positive going half-wave across many recording sites. The resulting segmentations are shown in Figure 5C for condition (22-a) and Figure 5D for condition (22-b) after additional Hausdorff clustering (θ = 0.25) for alignment between conditions. Moreover, we present the centers of gravity Equation (10) in Figure 5E for condition (22-a) and Figure 5F for condition (22-b) of the corresponding segments. The ERPs of both conditions start in a metastable baseline state 1 in the pre-stimulus interval. After a first transient (dark blue) still in the pre-stimulus interval both ERPs proceed into another baseline state 2. The P200 attentional component is realized as segment 3 in both conditions. Yet in condition (22-b) also the earlier N100 is detectable as MS 7 that does not exist in condition (22-a). The remaining segments 4, 5, 6 in the control condition (22-a) exhibit rather spatially flat topographies (cf. Figure 5E,F). While MS 6 is common to both conditions, the final state 8 in condition (22-b) reflects the crucial difference, namely the posteriorly distributed positive P600 component. These results are in full agreement with the understanding of language processing.

After identification of the MS, we construct the neural field model based on results of Section 2.2. First we present the functional connectivity kernels w^(c)(x, y) for conditions (22-a) (c = 1) and (22-b) (c = 2) in Figures 6A,B. Figure 6C additionally displays the kernel difference w⁽²⁾(x, y) − w⁽¹⁾(x, y). For better visualization, we have ordered the recording electrodes according to their hemispheric topography, thus creating nine “regions of interest” (ROI): LF: left-frontal, LT: left-temporal, LP: left-parietal, LO: left-occipital, C: central (midline electrodes), RF: right-frontal, RT: right-temporal, RP: right-parietal, RO: right-occipital².

Figure 6 reveals a kind of checkerboard texture for the kernels w^(c)(x, y). Mostly obvious is the strength along the main diagonal, indicating mainly self-connections within brain areas, such as in ROIs LO and RO. Since primary visual cortex V1 is situated at occipital brain areas, these patterns may be interpreted as reflecting the visual presentation paradigm of the experiment. Interestingly, also the kernel difference exhibits largest differences in LO-LO and RO-RO areas. There is also a hemispheric asymmetry between reciprocal connections of left-parietal (LP) with left-occipital (LO) and right-parietal (RP) with right-occipital (RO) areas. It is quite tempting to speculate that this asymmetry is due to the brain's asymmetry regarding language processing: The largest difference between conditions is in the RF area, which is strong for the correct condition (22-a), but rather weak for the phrase structure violation condition (22-b). This might be seen as a neural correlate of semantic and pragmatic integration processes that are supported by the right hemisphere [50]. These processes fail in the case of the violation condition.

Next we show the simulated spatiotemporal dynamics of the neural field simulation from Section 2.2. To this end, the dynamics is illustrated as a temporal snapshot sequence of spatial topographies.

Figure 7 shows the simulation results for the correct condition (22-a) where every snapshot presents the neural field activation pattern at the indicated point in time. The color scale range is [−10, 10]μV.

The sequence exhibits one prominent pattern, a frontal positivity around 235 ms, reflecting the attentional P200 component. Hence it resembles well the original time series shown in Figure 1.

Accordingly, we present in Figure 8 the simulated neural field for the phrase structure violation condition (22-b). Again the color scale ranges between [−10, 10]μV.

In this condition the P200 effect is visible as a frontally pronounced positivity between 331 and 426 ms, however that is delayed compared to the original ERP data (cf. Figure 2). Hence our simulation entails a desynchronization in comparison with the experimentally observed ERP dynamics. The reason for this deviation is the presence of only one time scale in the underlying dynamic model, reflected by the growth rates of neural populations in the Lotka-Volterra Equation (14) that are all of the same order of magnitude. Hence, the phasic MS 7 is not appropriately captured by our phenomenological model.

Moreover, the earlier N100 component with reversed polarity is now shifted toward the time window between 139 and 235 ms. The final parietal positivity starting off at about 470 ms in the ERP data is already present in our simulation, now starting at 618 ms until the end of the epoch window at 952 ms.

As above, a proper interpretation of ERP effect requires the computation of condition differences. We plot the simulated neural field difference between violation condition (22-b) and correct condition (22-a) in the range [−10, 10]μV in Figure 9 similar to the ERP difference plot shown in Figure 3.

A strong artifact effect from 139 until 235 ms is visible that is due to the misalignment of N100 and P200 in both conditions. However, the difference patterns clearly indicate the posteriorly distributed P600 ERP component that evolves between 618 and 1,000 ms post-stimulus.

For the recurrence structure analysis of the neural field simulation we optimize the Markov utility function (3) for the grand averages of both conditions (22-a) and (22-b) separately. The resulting utility functions are depicted in Figure 10.

Figure 10 indicates that both conditions lead to optimal segmentations for ε^* ≈ 0.0051.

The results of the recurrence structure analysis (RSA) are shown in Figure 11.

Figure 11A displays the simulated neural fields for the correct condition (22-a) and Figure 11B for the violation condition (22-b) where each colored trace denotes one simulated electrode. In both simulations the initial conditions were prepared as the stationary base line state 1 that does not change during the pre-stimulus interval. At stimulation time t = 0 a slight perturbation kicks the state out of equilibrium, triggering the sequential heteroclinic dynamics. The resulting segmentations are shown in Figure 11C for condition (22-a) and Figure 11D for condition (22-b) after additional Hausdorff clustering (θ = 0.1) for alignment between conditions. Despite the timing differences the pattern sequence is essentially the same as in Figure 5. This is confirmed by the centers of gravity Equation (10) in Figure 11E for condition (22-a) and Figure 11F for condition (22-b) of the corresponding segments. The simulated ERPs of both conditions start in a metastable baseline state 1 in the pre-stimulus interval. After a first transient (dark blue) still in the pre-stimulus interval both ERPs proceed into another baseline state 2. The P200 attentional component is realized as segment 3 in both conditions. In condition (22-b) also the earlier N100 is detectable as MS 7 yet with much longer duration than in the experimental ERP analysis. The remaining segments 4, 5, 6 in the control condition (22-a) exhibit rather flat topographies. While MS 6 is common to both conditions, the final state 8 in condition (22-b) reflects the crucial difference, namely the posteriorly distributed positive P600 component.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.