Dynamic decomposition of spatiotemporal neural signals

Neural signals are characterized by rich temporal and spatiotemporal dynamics that reflect the organization of cortical networks. Theoretical research has shown how neural networks can operate at different dynamic ranges that correspond to specific types of information processing. Here we present a data analysis framework that uses a linearized model of these dynamic states in order to decompose the measured neural signal into a series of components that capture both rhythmic and non-rhythmic neural activity. The method is based on stochastic differential equations and Gaussian process regression. Through computer simulations and analysis of magnetoencephalographic data, we demonstrate the efficacy of the method in identifying meaningful modulations of oscillatory signals corrupted by structured temporal and spatiotemporal noise. These results suggest that the method is particularly suitable for the analysis and interpretation of complex temporal and spatiotemporal neural signals.

where the coefficients c k are chosen in a way to have stable solutions. An important tool for analyzing a linear differential equation is the impulse response function G(t). This function is defined as the response of the system to a unit-amplitude impulse δ(t): Using the impulse response function, a solution of the linear SDE driven by an arbitrary random input w(t) can be written as follows: This means that the stochastic process α(t) is an infinite linear superposition of responses to the random uncorrelated input w(s). Using Eq.
(3) we can derive the mean and covariance function of α(t). The mean function is defined as where the triangular brackets · denote the expectation with respect of the distribution of the random input w(s). Using (3) in (4), we obtain: Here, we used the fact that the order of expectation and integration can be interchanged and that the expectation of the white noise process is equal to zero. Analogously, we can obtain the covariance function as follows: Since w(s) is white, its covariance w(s)w(s ) is given by the delta function σ 2 α δ(s − s ), where σ 2 α is the variance of the random input. The integral over s can be solved by using the translation property of the delta function: Using this formula and introducing the new integration variable s * equal to t − s, the covariance function becomes Since the covariance function depends on t and t only though their difference τ = t − t , we denote it as k (τ ).

S2 Appendix: Spherical harmonics and spherical Fourier transform
Spherical harmonics are the generalization of sine and cosine on the surface of a sphere. They are parametrized by the integers l and m, of which l is a positive integer and m ∈ {−l, . . . , l}. These two parameters determine, respectively, the angular frequency and the spatial orientation. Spherical harmonics are defined by the following formula: where P |m| l is a Legendre polynomial [1]. Spherical harmonics form a set of orthonormal basis functions and, consequently, we can use them to define a spherical Fourier analysis [2]. Specifically, the spatiotemporal process α(x, t) can be expressed as a linear combination of spherical harmonics whereα(l, m; t) is the l, m-th spherical Fourier coefficient as a function of time, defined asα (l, m; t) = C α(l, m; t)H m l (x)dx.
Eqs. (10) and (11) are the equivalent of respectively inverse and direct Fourier transform for functions defined on the surface of a sphere.

S3 Appendix: Properties of the Kronecker product and GP regression with separable covariance matrices
In order to derive the posterior expectations of the spatiotemporal GP regression, it is useful to introduce some of the properties of the Kronecker product between matrices. The Kronecker product between two N × N matrices is defined by the block form: The following formula relates the regular matrix product with the Kronecker product: The inverse and transpose of a Kronecker product are respectively and The following formula relates the Kronecker product to the vectorization of a matrix: Using these formulas, we can now derive the posterior expectation of the spatiotemporal GP regression. Combining the spatiotemporal prior and the observation model using Bayes' theorem, we obtain the posterior This is the product of two multivariate Gaussian densities and it is therefore a multivariate Gaussian itself. Its expectation is given by Using (13) and (15), the expression simplifies to: This formula involves the inversion of a matrix that is the sum of two Kronecker product components. Inverting this matrix would be computationally impractical. We simplify the problem by imposing Σ = ΛDΛ T . In this case, Eq. (14) allows to invert the spatial and temporal covariance matrices separately: In most realistic cases, the MEG observation model Λ will not be full rank, therefore we introduced a Tikhonov regularization parameter λ.
Using Eq. (16), we finally arrive at the posterior expectation:

S4 Appendix: Modeling vector-valued sources using block matrices
In the main text, the source reconstruction formulae are expressed for fixed dipole directions v(x). The solution for the general case, in which the dipole direction is estimated from the data, is obtained by introducing an independent set of spherical