Data-driven Koopman operator approach for computational neuroscience

This article presents a novel, nonlinear, data-driven signal processing method, which can help neuroscience researchers visualize and understand complex dynamical patterns in both time and space. Specifically, we present applications of a Koopman operator approach for eigendecomposition of electrophysiological signals into orthogonal, coherent components and examine their associated spatiotemporal dynamics. This approach thus provides enhanced capabilities over conventional computational neuroscience tools restricted to analyzing signals in either the time or space domains. This is achieved via machine learning and kernel methods for data-driven approximation of skew-product dynamical systems. The approximations successfully converge to theoretical values in the limit of long embedding windows. First, we describe the method, then using electrocorticographic (ECoG) data from a mismatch negativity experiment, we extract time-separable frequencies without bandpass filtering or prior selection of wavelet features. Finally, we discuss in detail two of the extracted components, Beta (∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ \sim $\end{document} 13 Hz) and high Gamma (∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ \sim $\end{document} 50 Hz) frequencies, and explore the spatiotemporal dynamics of high- and low- frequency components.

One challenge stems from the complexity of brain activity patterns and the different frequencies characterizing these patterns. Most common methods of analysis rely on predefined parameters: e.g., the width of the band in bandpass filtering, the wavelets shape in wavelet analysis, or the linear operators in principal component analysis. Once necessary parameters are selected, data are projected onto them, and the output -a subset of data that fits our parameters -is examined to capture spatial and temporal dependencies within the brain's activity [6]. In other words, most available approaches require the researcher to impose conditions extrinsic to the dynamical system governing the data, and as a result, risk inaccurate representations of the states of that system.
The risks entailed in relying on bandpass filtering to separate various frequencies have been well documented [3]. Nonetheless, most analyses of EEG and ECoG data begin with filtering available datasets, a process that inevitably leads to some degree of departure from the original signal. As a result, the use of different filter parameters may result in very different outcomes [63].
Following the filtering process, eigenfunctions of linear operators are commonly used for performing dimension reduction to explore brain activity across channels or time. One example of such -principal component analysis (PCA) -employs eigenvectors of the empirical signal covariance matrix to recover independent patterns. The dynamical process underlying the observed data does not constrain the operators and can lead to the extraction of patterns with limited ties to the underlying physical mechanisms. For example, it has been established that PCA has limited ability to recover distinct frequencies from broadband signals, an especially important shortcoming in the presence of components that can account for a large portion of observed variance [4]. Moreover, by their construction, PCA patterns are linear functions of the input data, which constrain the generality of patterns that the method can recover, particularly when the observations do not provide access to all states of the system (typically the case in neuroscience applications).
Another class of techniques takes advantage of various wavelets (e.g., Hilbert-Huang, Morlet wavelets) to transform the signal into frequency domain (e.g., [11]) and select dominant modes. This approach allows researchers to move beyond the assumption that the signal is a sine wave. Yet, this too requires a prior choice of the wavelet transformation. Furthermore, additional tools are needed to allow for predictions regarding the behavior of the selected frequencies.
Another shortcoming of the above approaches is related to the character of electrophysiological data. While EEG and ECoG records patterns that evolve in space and time, most available computational methods study temporal and spatial changes separately [6].
Despite their shortcomings, these approaches have been widely used to study one of the central questions of contemporary neuroscience: What mechanisms underlie spatial and temporal coordination of neural activities across various brain areas and during different steps of cognitive processing? And while several models of spatiotemporal coordination were successfully developed [37], there is an urgent need for computational methods that overcome the shortcomings of standard analytical methods [6].
Recent years have seen increased interest in analysis techniques treating the data as "observables" -functions of the state of the underlying dynamical system [49]. This allows researchers to combine ideas from the spectral theory of dynamical systems with modern machine learning approaches to construct algorithms for signal decomposition that make minimal ad hoc assumptions of the input data. They target dynamically intrinsic quantities, making minimal assumptions on the structure of the data.
In this work, we introduce the applications of a data-driven approach that can successfully remedy the above limitations. Specifically, we focus on a Koopman operator approach that: (1) is data-driven, and the extracted components are intrinsic to the underlying dynamical system; (2) does not require bandpass filtering; (3) captures high-dimensional data across both spatial and temporal domains; (4) provides means to visualize the spatial and temporal coordination of neural activity; (5) decomposes complex dynamics into a hierarchy of periodic and quasi-periodic patterns even if the underlying dynamics is highly non-periodic.
The plan of this paper is as follows: Section 2 begins with a description of Koopman operator formalism and the rationale behind data-driven approximation of the eigenfunctions of Koopman operators using kernel methods for machine learning (Sections 2.2 and 2.3). Section 2.4 provides algorithms applied to reconstruct the spatiotemporal patterns of the Koopman eigenfunctions. Section 3 applies the Koopman operator approach to ECoG data and provides a description of the neural data (from a marmoset monkey) and the experimental design. Section 4 lists the parameters used in the analyses (number of nearest neighbors, number of examined samples, kernel's density distribution, etc.) and discusses the results for two of the decomposed patterns with frequency bands matching Beta and Gamma. This section also discusses potential limitations. Section 5 provides a short summary and briefly describes an extension of the current approach called vector-valued spectral analysis.

Data-driven approach to computational neuroscience
We apply a recently developed framework [15,23,24,32] inspired from the operatortheoretic formulation of ergodic theory that addresses many of the limitations of PCA and related approaches outlined above. Instead of estimating a covariance operator, our approach is based on data-driven approximations of the Koopman operator: the fundamental operator governing the evolution of observables (functions of the state) of a dynamical system.
It is a remarkable fact, realized in the work of Koopman in the 1930s [43], that the action of a general nonlinear dynamical system can be characterized without approximation by intrinsically linear operators acting on vector spaces of observables. These operators provide a natural and theoretically rigorous framework for spectral decomposition and statistical prediction in complex systems. In particular, the eigenfunctions of Koopman operators factor the dynamics, which may be chaotic, into periodic or quasi-periodic modes. These modes can be thought of as analogs of Fourier modes tailored to the dynamical system generating the data.
While this operator-theoretic framework has long been known, it has only fairly recently been employed in data-driven approaches [16,49]. A special feature of the techniques employed here is the use of kernel methods for machine learning [5,10,27] to approximate the eigenvalues and eigenfunctions of the Koopman operator with rigorous convergence properties and minimal assumptions on the system and measurement modality.
Elsewhere [30,59], we have applied this approach to identify a multi-scale hierarchy of modes of variability from sea surface temperature data in the Earth's climate system, spanning months to decades. There are several important similarities between spatiotemporal climatic data and multi-channel ECoG recording of brain potentials, including a two-dimensional surface for sampling the system, as well as meaningful signals that vary on multiple spatial and temporal scales. Here, we apply this approach and describe two of the extracted modes, previously shown to characterize neural responses to surprising sensory events [47,48].

Koopman operators formalism
We approach the problem of extracting spatiotemporal patterns of neuronal dynamics from the point of view of the operator-theoretic formulation of ergodic theory [19]. In particular, we treat electrophysiologic data as observables of a skew-product dynamical system. By that, we mean that there is a dynamical system, with state space X and flow map Φ t : X → X, t ∈ R, driving another dynamical system on a state space Y . Here, (X, Φ t ) describes the evolution of the external stimuli, and Y is the state space for the brain. The dynamics of the latter are modeled via a map Ψ t : X × Y → Y such that y t = Ψ t (x, y) corresponds to the state of the brain at time t, given that at time 0 the external stimulus had the state x ∈ X, and the brain the state y ∈ Y . Note that in this representation the brain dynamics are manifestly non-autonomous since y t depends on both y and the state x of the driving system. On the other hand, the product dynamical system on M = X × Y is autonomous, and can be described by the evolution map Ω t : M → M, where Ω t (x, y) = (Φ t (x), Ψ t (x, y)).
The dynamical system (M, Ω t ) is said to have a skew-product structure as component X affects the dynamical evolution of component Y , but Y does not influence X. In this description, the ECoG signal h ∈ R d recorded at d sensors on the marmoset cerebral cortex when the brain state is y ∈ Y can be described via an observation function F : Under mild assumptions on the state spaces X and Y , and the evolution maps Φ t and Ψ t , the product dynamical system (M, Ω t ) will possess ergodic invariant measures. That is, given a dynamical trajectory (x n , y n ), where is the state in M reached at time t n = nτ for a fixed sampling interval τ > 0 and initial conditions (x 0 , y 0 ) ∈ M, there exists a probability measure μ on M, invariant under the flow map Ω t , such that for every μ-integrable function f : M → C, the time averagef N = N −1 N−1 n=0 f (x n , y n ) converges almost surely to the expectation valuef = M f dμ. Note that the system can have multiple ergodic components, but data collected from a given experiment sample a specific ergodic component of the system. Associated with the triplet (M, Ω t , μ) is a Hilbert space H = L 2 (M, μ) of squareintegrable observables with respect to μ and a group of unitary Koopman operators U t : H → H governing the evolution of observables under Ω t . That is, given f ∈ H , g = U t f is defined as the observable satisfying g(x, y) = f (Ω t (x, y)) for μ-almost every (x, y) ∈ M. For our purposes, the space H will be the space of admissible temporal patterns generated by the marmoset neuronal dynamics. In particular, for every f ∈ H , the mapping t → U t f (x, y) defines a temporal pattern associated with a sampling along the dynamical trajectory starting at (x, y) ∈ M. Our approach is to identify dynamically intrinsic temporal patterns through approximate eigenfunctions of the Koopman operator. This will lead to a decomposition of the observation map F into a corresponding set of spatiotemporal ECoG patterns.

Koopman eigenfunctions
Significant temporal modulation frequencies can be identified by computing the Koopman eigenfunctions of a dynamical system and this can provide a better representation of the system than standard approaches such as discrete Fourier transform and PCA. In the setting of interest here, an observable z j ∈ H is a Koopman eigenfunction if it satisfies the eigenvalue equation for all t ∈ R. In the above, ω j is a real-valued frequency associated with eigenfunction z j . In particular, ω j is real due to the fact the dynamics Ω t preserves the measure μ.
A key property of Koopman eigenfunctions of measure-preserving systems is that, along a given dynamical orbit, they vary according to a multiplication by a periodic phase factor. Specifically, it is a direct consequence of (2) that for μ-a.e. point (x, y) ∈ M, U t z j (x, y) = z j (Ω t (x, y)) = e iω j t z j (x, y). This means that, whenever they exist, Koopman eigenfunctions give rise to temporal patterns with a highly coherent and predictable temporal evolution, each recovering a timescale 2π/ω j intrinsic to the dynamical system. This behavior is markedly different from that of the eigenfunctions of covariance operators approximated by PCA, which in general have no direct connection to the spectral properties of the dynamical system. In addition, compared to spectral estimation techniques based on the discrete Fourier transform, which require identification of the "true" frequency peaks from a set of candidate frequencies, an advantage of pattern extraction via (2) is that an estimate of the eigenfrequencies ω j is a direct output of the method. It can further be shown that (1) Koopman eigenfunctions corresponding to distinct eigenfrequencies are orthogonal on H (that is, z i , z j H = M z * i z j dμ = δ ij , where ·, · H denotes the Hilbert space inner product); (2) every Koopman eigenspace is one-dimensional; and (3) the spectrum {ω j } of Koopman eigenfrequencies is countable.
The closed span of all Koopman eigenfunctions identifies a subspace D = span{z j } ⊆ H , invariant under U t , leading to an associated orthogonal decomposition H = D ⊕ D ⊥ . Observables lying in D have an associated Koopman eigenfunction expansion, which can be thought of as a generalization of the Fourier transform of time-periodic signals. Specifically, every f ∈ D admits the decomposition f = jf j z j , wheref j = z j , f H are complexvalued expansion coefficients. Moreover, the dynamical evolution of f can be computed in closed form via It should be noted that the evolution of observables in the orthogonal subspace D ⊥ cannot be determined on the basis of a countable set of eigenfrequencies as in (3). Instead, the action of Koopman operators on such observables can be determined via an expansion involving its continuous spectrum [49]. In general, systems encountered in real-world applications are of the mixed spectrum type, meaning that they exhibit both eigenvalues and continuous spectrum (i.e., D contains nontrivial observables, but it is a strict subspace of H ). Due to their desirable temporal coherence and predictability properties, Koopman eigenfunctions and their corresponding eigenfrequencies are natural objects to identify for the purpose of analyzing complex systems.

Data-driven approximations
In order to compute solutions to the Koopman eigenvalue problem in (2), we first construct a data-driven basis of H using a class of kernel algorithms, called nonlinear Laplacian spectral analysis (NLSA), operating on delay-embedded data [26][27][28]. Selecting an integer parameter Q (the number of delays), we first define the augmented observation map F Q : and then define a kernel function k Q : M × M → R + , which provides a pairwise measure of similarity of points in M, based on the observation function F Q . For the purposes of this discussion, we consider a radial Gaussian kernel, as a concrete example, but in practice we work with anisotropic Gaussian kernels [22,[26][27][28] that also take into account time tendencies of the data. In the above, is a positive bandwidth parameter controlling the rate of decay of the kernel, and · the canonical Euclidean norm on R Qd . Given a finite dataset consisting of N + Q snapshots h −Q+1 , . . . , h N−1 , with h n = F (a n ) taken at the dynamical states a n = (x n , y n ) from (1), the kernel values k Q (a m , a n ) with m, n ∈ {0, . . . , N − 1} can be empirically computed by first forming the delay-embedded data vectors h n,Q = (h n , h n−1 , . . . , h n−Q+1 ) ∈ R Qd , and then evaluating k Q (a m , a n ) = e − h Q,m −h Q,n 2 / . Using these pairwise kernel values, we then compute the Markov kernel p Q,N (a m , a n ) = k Q (a m , a n ) l Q,N (a m )r Q,N (a n ) , introduced in the diffusion maps algorithm [10], where r Q,N and l Q,N are normalization functions given by k Q (a m , a n ), l Q,N (a n ) = N−1 n=0 k Q (a m , a n ) r Q (a n ) .
Note that p Q,N has the Markov property, N−1 n=0 p Q,N (a m , a n ) = 1, by construction. We then form the N × N ergodic Markov matrix P Q,N = [p Q,N (a m , a n )], and compute its eigenvalues and eigenvectors, For p Q,N defined as in (6), the eigenvalues of P Q,N are real and thus lie in the interval [−1, 1] by Markovianity. By convention, we order the eigenvalues in decreasing order, The eigenvector φ 0 has all components equal to 1 by standard properties of Markov matrices. As is standard practice, we take advantage of the exponential decay of the Gaussian kernel to sparsify the matrix P Q,N , retaining approximately k nn N nonzero entries per row based on the k-nearest neighborhood of the corresponding data point. This step significantly increases the scalability of the approach to large N , since apart from data with special structure, it would be difficult to carry out using kernels without a localizing behavior (e.g., covariance and polynomial kernels).
Each eigenvector φ j = (φ 0j , . . . , φ N−1,j ) ∈ R N can be thought of as sampling a continuous function ϕ j on M, such that ϕ j (x n , y n ) = φ nj . It can be shown that in the limit of large data, N → ∞, these functions converge to eigenfunctions of a Markov operator P Q : H → H , where and p Q is a Markov kernel defined analogously to p Q,N in (6) [15]. Importantly, as Q increases, every eigenspace of P Q corresponding to nonzero eigenvalue converges to a finite union of Koopman spaces [15,23]. This property makes the eigenfunctions produced by NLSA a highly efficient basis for solving the Koopman eigenvalue problem. Note that the ability of NLSA to approximate Koopman eigenspaces in the Q → ∞ limit does not depend on the kernel being Gaussian; it is, rather, a consequence of delay embedding that would also hold for other kernels, including covariance kernels constructed from delayembedded data as in singular spectrum analysis (SSA) algorithms [34]. An advantage of the nonlinear Gaussian kernels over linear covariance kernels is that, in the former case, the subspace spanned by eigenfunctions corresponding to nonzero eigenvalues (which amounts to the subspace that can be consistently approximated from data) is infinite-dimensional, irrespective of Q being finite or infinite. On the other hand, in the latter case, the subspace is finite-dimensional (in effect, the Gaussian kernel can be viewed as an "infinite-degree" polynomial kernel). For our purposes, this means that Gaussian kernels can provide a richer set of basis functions to approximate the Koopman operator than covariance kernels. Following [15,23,32], we employ the data-driven basis {φ j } from NLSA to compute approximations z j ∈ C N , andω j ∈ C to the Koopman eigenfunctions z j and the corresponding eigenfrequencies ω j , respectively, using a Galerkin method. See [25, Table 3] for pseudocode. In all cases, we order Koopman eigenfunctions in order of increasing Dirichlet energy: a non-negative quantity that intuitively measures their "roughness" as functions on M. The rationale for this ordering (which does not necessarily correspond to ordering with respect to frequency) is that the functions with low Dirichlet energy are more amenable to robust approximation from finite datasets, reducing the risk of sampling errors. Details on the Koopman Galerkin framework can be found in the references just cited. By virtue of the spectral convergence results on the (λ j , φ j ) stated above, the solutions of this data-driven Galerkin method can be shown to converge to eigenfunctions of an approximate Koopman operator, with a small amount of diffusion added to render the Koopman spectrum discrete. An aspect of the present analysis, which has not been discussed in earlier introductions of the method (e.g., [23]), is that we consider a skew-product dynamical system appropriate for non-autonomous neuronal dynamics. The source code used here is available for download at the third author's personal website (http://cims.nyu.edu/ ∼ dimitris), while pseudocode can be found in [23].

Spatial and temporal reconstruction
We now describe how the Koopman eigenfunctions z j can be used to decompose the observed ECoG signal h n into coherent spatiotemporal patterns. The reconstruction method is similar to the procedure used in SSA [34] and NLSA [25][26][27][28].
First, we compute the time-lagged projections of the observable F onto the Koopman eigenfunctions z j , given by where N = min{N, N + q}, and z nj is the n-th component of eigenvector z j . That is, A j (q) is an ECoG pattern in R d approximating the projection of F , shifted by q sampling time intervals, onto Koopman eigenfunction z j . In the Koopman operator literature, A j (0) is referred to as the Koopman mode associated with the observable F [7,49], though note that here we will use A j (q) with several values of q to perform reconstruction. In particular, following [25][26][27][28]34], we define the spatiotemporal pattern associated with the pair (z j , A j ) as the function F j : M → R d such that For eigenfunctions forming complex-conjugate pairs (z j , z j ) (which is is usually the case), we compute the sum F j + F j . The latter is real since the observation map F is real.

Mismatch negativity
Mismatch negativity (MMN) is a drop in the recorded potential observed when a change in the external environment takes place. Some examples of such changes are the duration of, or time between, sounds in the surroundings [21,51,52]. To better explain the external conditions, let us consider an example of an animal listening to short tones occurring every 40 seconds. In this example, we would expect to observe MMN in the case when the tone unexpectedly appears after only 20 seconds. MMN is characterized by a strengthening of the negative peak following the deviation from the preceding sequence [1]. Depending on the animal, following the deviant stimuli, the wave of MMN peaks anywhere between 30 ms (e.g., in rodents [54]) to 250 ms (in humans) [1,14,39,42,62]. While MMN is one of the best described neural markers of unexpected stimuli and changes in the environment [35], explanations for the mechanisms driving it can differ drastically [46,62]. This lack of agreement might be driven by the variety of brain areas contributing to MMN, among them the auditory cortex (AUD) [1,50], the frontal cortex [1,45], and the dorsolateral prefrontal cortex [2,45]. Furthermore, numerous studies suggest a significant temporal breadth of MMN sub-components that appear in different stages of MMN [1,33]. Though MMN measures can indicate normal attentional states, there is growing appreciation for the use of MMN to identify endophenotypes for neurological disorders, including psychiatric conditions [53,56]. For example, MMN anomalies have been observed among individuals with schizophrenia [40,41,56], as well as sleep-deprived individuals [36].

Data description
To test the applicability of the Koopman operator approach described in Section 2 for neural data, we used a subset of ECoG data from an experiment designed to evoke MMN in marmoset monkeys through changes in the frequency of sound waves. Data were obtained from http://neurotycho.org/. This open-source database contains recordings from a study designed and conducted by Komatsu, Takaura, and Fujii [42] for one of two subjects in the experiment. Along with ECoG recordings, the authors provide sound stimulus parameters, and other necessary information to replicate their analysis (we refer readers to [42] for a more detailed description of the study's design and data acquisition protocol). ECoG data consist of recordings from 32 electrodes placed across the monkey's left hemisphere (Fig. 1) with a sampling rate of 1000 Hz. The total duration of the experiment was 15 minutes. During this period, the subject was passively listening to trains of tones. Each tone lasted for 64 ms. The time between tones was constant (439 ms), so each interval lasted for 503 ms. In total, the subject was exposed to 240 sound trains. Each train consisted of repetitions of the same tone. The trains differed in tones' frequency (from 250 Hz to 6727 Hz) and the number of times the tone was repeated (3, 5, or 11). Consistent with [42], we considered the first tone "deviant" and the last tone "standard". To validate the Koopman operator framework, we replicated and described findings presented by Komatsu and colleagues [42] in [47]. In this paper, consistent with past research highlighting the importance of both Beta and Gamma frequency oscillations in response to unexpected auditory stimuli [37], we extend our earlier work [47] describing two of the extracted modes that fall in the Beta and Gamma frequency range.

Default setup and robustness of results
Here, we examine 45,266 samples (t 0 = 147,359 ms to t 1 = 192,625 ms) that include 13 trains with a sampling rate of 500 Hz. This dataset contains a total of 89 tones with tonal frequency ranging form 250 Hz to 4000 Hz. Similarly to [47], this selection is arbitrary but allows us to obtain a balance between a sufficient number of samples and computational cost related to subsequent analysis.
We computed Koopman eigenfunctions as described in Section 2 using the anisotropic Gaussian kernels introduced in NLSA and the family of cone kernels introduced in [22]. Ultimately, we found our results to be robust against the type of kernel employed. However, we selected the NLSA kernel as our primary focus here to evaluate faster (i.e., corresponding to MMN duration) dynamics. The choice of the NLSA kernel bandwidth parameter was = 25, and the k-nearest neighborhood for each data point has a size of k nn = 5000. The specific value of is rather large versus other applications [30,59], and the possible reason is the more intermittent character of the ECoG signal as compared to fluid flows. These, in turn, can be related to under-sampling of the MMN signal over the analyzed period of the experiment. The same logic applies to a small number of nearest neighbors. For largeenough and small-enough k nn , however, the results presented here seem to sample the local neighborhood robustly and sufficiently when it comes to extraction of modes needed to capture response to deviant sound.
A range of tests led us to conclude that 250 time-steps (equal to 500 ms) is an appropriate choice for the embedding window given the short timescale of the modes of interest (tens of milliseconds). An insufficiently large embedding window does not allow clear separation of different frequencies, as expected theoretically, while an overly long embedding window introduces low-frequency modes (not of interest here, thus not displayed). Long embedding windows also split high-frequency modes, introducing redundancy (given the purpose of our study). As such, we show that a theoretically meaningful response to "deviant" sounds appears robust within a suitable range of the embedding window (see Figs. 2 and 3). 3) in the hierarchy of modes retrieved by Koopman analysis; specifically, they correspond to the complex-conjugate pairs {z 5 , z 6 }, {z 9 , z 10 }, and {z 13 , z 14 } (henceforth abbreviated z 5,6 , z 9,10 , and z 13,14 , respectively). These response frequencies correspond to the Beta band in several primate species including the marmoset [8,37,61]. Our result thus confirms prior studies that find a change in Beta frequency associated with sound deviants or natural dynamically modulated sounds [8,37]. However, here Fig. 2 Real part (left) and its frequency spectrum (right) for Koopman eigenfunctions z 5,6 (blue), z 9,10 (red), and z 12,13 (green) we identify multiple peaks within the Beta range without having to select an arbitrarily defined range of modulation frequencies. Among the pairs z 5,6 , z 9,10 , and z 12,13 examined here, eigenfunction pair z 5,6 is the most periodic in nature, exhibiting a prominent 13 Hz (77 rad/s) spectral peak. Eigenfunctions z 9,10 and z 12,13 have frequency peaks at 10 Hz and 16 Hz, respectively. Their temporal patterns are more irregular and have a less prominent low-frequency envelope. Such decomposition of complex signals into periodic and irregular patterns has been observed in the case of other dynamical systems [57][58][59], but it can be achieved with sufficiently long embedding windows (this is also true for traditional methods such as SSA).

Temporal coordinates
In [42], the differences between evoked-response potentials following the last and first tone of each train were averaged across all trains for each electrode, and the differences in averages were tested for statistical significance. Figure 3 depicts reconstructed patterns attained via the Koopman operator approach, aligned similarly for electrode 25 that corresponds to the primary auditory cortex, A1, according to [42]. The extracted pattern was averaged for 13 consecutive onsets of deviant tones and 13 preceding standard tones. The figure shows the reconstructed signal for electrode 25 obtained via the procedure described in Section 2.4 using eigenfunctions z 5,6 , z 9,10 , and z 13,14 . To extract an average response to deviant sound, the reconstructions are averaged over the 13 response trials analyzed here. The average corresponds to the more periodic conjugate pair of the temporal pattern z 5,6 after subtracting the overlap with the more chaotic z 9,10 and z 12,13 . It is worth noting that this projection captures amplification of the signal shortly after the appearance of the "deviant" sound. In particular, a large positive peak is observed at approximately 10 ms followed by a large negative peak at 25-50 ms. A second large negative peak is observed around 125 ms after the "deviant" sound, and a third follows at approximately 200 ms. The  c At 30 ms, the negative valence signal strengthens and expands in the core and surrounding auditory cortex (electrodes [22][23][24][25][26][27][28][29], and the negative signal travels towards the parietal cortex (electrode 1). d At 50 ms, the negative valence signal in the parietal cortex is enhanced and propagates from electrode 1 to the rest of the brain. e Then at 92 ms, a second pair of high-amplitude, positive valence peaks appears in about the same positions as seen for 10 ms. f A second negative valence peak follows at approximately 130 ms temporal dynamic of the mode suggests distinct onset and terminating responses following deviant sound onset that matches well with the dynamics observed in core or primary auditory cortices in other studies [20,38,44,55]. The reader should also notice that all signs, amplitudes, and time relative to the signal onset of the deviant match well with those depicted in Fig. 2 of [42]. Subsequent negative and positive events are followed by the particular signal weakening (with slight strengthening observed at t = 503 ms), which corresponds to the standard tone occurrence (Fig. 3, grey line, 503 ms).
To summarize, we observe that spatiotemporal reduction of a multidimensional signal into three coordinates derived from approximate Koopman eigenfunctions permits reconstruction of an intermittent neural response and its temporal emergence following the deviant signal (Fig. 4).

Statistical tests
We tested statistical significance for 12 electrodes identified by the authors of the experiment [42] following the guidelines for the extracted mode. Specifically, we tested the average difference between ERPs of the standard and deviant stimuli following the onset of the tone, adjusting p-values with the false discovery rate (FDR) and using Wilcox ranked sum tests. The results replicate the findings presented in Fig. 3a of [42]. Notably, the pvalues presented in Table 1 are smaller than the p-values in [42]. This reduction in the probability of type I error is most likely driven by the reduction in the error term when comparing amplitude differences of a single decomposed frequency. As mentioned above, another challenge in analyzing EEG and ECoG data is the limitation faced by variance-based approaches (e.g., PCA) when capturing low-amplitude signals. These approaches often fail to adequately collect such low-amplitude components in the presence of a strong signal that accounts for a large amount of variance in the observed data [4]. This limitation is particularly problematic when studying high-frequency, lowamplitude oscillations, such as Gamma signals. We previously found that the data-driven Koopman operator approach is a robust tool for identifying low-amplitude events that are masked by larger-amplitude events [27,48]. Here, to further evaluate the usefulness of this approach for analysis of neural data, we examine sound evoked activity in the Gamma frequency band, which typically has low amplitude.
It has been established that Gamma is another frequency band that oscillates in response to an unexpected auditory stimulus [18,37]. Synchronous cortical activity in higherfrequency bands like Gamma can be driven by a single sensory modality and are restricted spatially to a single cortical area and temporally to short periods of time [60]. These properties make it difficult to localize Gamma related activity. In contrast, synchronous cortical activity within slower frequency bands can be driven by multiple sensory modalities, is observed over a large spatial area (i.e., multiple cortical fields) and lasts for a longer period of time [18,60].
Here, we illustrate the time evolution of a Koopman eigenfunction pair z 141,142 and its frequency spectrum (Fig. 5). Note that this mode has a higher frequency (50 Hz) than the Beta band activity associated with the lower-frequency Beta modes (Fig. 2, z 5,6 , z 9,10 , and z 13,14 ).
As per [61], we expected the spatiotemporal dynamics of the reconstructed signal for z 141,142 to be localized. That is, a plot of the signal displayed over the electrodes matrix should exhibit more diversity in space than slower frequencies (e.g., Beta shown in Fig. 4). Indeed, as shown in Fig. 6, the reconstructed signal often exhibits positive amplitudes in some areas of the cortex while exhibiting negative amplitudes in others. Additionally, the low amplitude of the signal aligns our results with the predictions based on the existing body of scientific literature.
In this work we did not cross-validate the approach with PCA or wavelet analysis methods as applied in neuroscience [17]. A more extensive exploration of the Koopman operator approach to ECoG data could benefit from cross-validating the results with other methods (e.g. wavelet analysis) and larger datasets. We thus recommend caution when approaching the results of the statistical tests summarized in Table 1. Here we have applied the method to data collected from a single marmoset over a subset of trials. And while the significance testing (Table 1) confirms that reported in [42] over a larger number of trials for the same subject, a dataset from a single subject carries a higher risk of false positives and limitations to external validity of the findings.

Summary
In this paper, we have demonstrated the feasibility of using a data-driven Koopman operator approach, implemented with kernel methods from machine learning (NLSA), in neuroscience to quantify spatial and temporal changes in cortical potentials in response to time-and-tone varying oddball sound sequences. Our study replicates results reported previously for the same cortical recordings [42] and uses the same traditional method of evoked potential averaging and significance testing to cross-validate sound-evoked effects on frequency bands isolated with NLSA and the Koopman operator.
Our results suggest that MMN can be captured by a reduced number of coordinates (temporal modes) derived from the eigenfunctions of an operator governing the evolution of observables under neuronal dynamics. These patterns describe brain dynamic shifts between quasiperiodic and nonperiodic regimes operating over a wide range of the frequency spectrum and can be operationally employed for predictive purposes [57]. Such intermittent dynamics are known to be poorly captured by Fourier methods and are better handled with wavelet transforms [11,12]. It is worth highlighting that the current approach can separate temporally coherent, and dynamically meaningful, frequency components with complex wave forms (also of potentially low mean amplitude) from high-dimensional time series. A recent extension of the method, called vector-valued spectral analysis [29,31,58], does not impose any spatiotemporal separability, which could allow for further refinement of the current approach.
In conclusion, the consistency between the current results and past research shows the capability of the data-driven Koopman operator approach to analyze neuroscientific data. In particular, the successful signal decomposition demonstrated here can help address several questions related to the underlying nonlinear dynamics and the puzzling statistics associated with them (asymmetry and skewness of dominant modes [30]). For example, using the remaining variance after subtracting the overlap with the more chaotic Koopman eigenfunctions z 9,10 , and z 13,14 , we may now extract orthogonal modes that may be a part of an overarching dynamics, even if their frequency spectra lie in very close proximity. Furthermore, the results suggest new biological dynamics associated with MMN response in two frequency bands.