The Role of Spectral Complexity in Connectivity Estimation

Abstract

The study of functional connectivity from magnetoecenphalographic (MEG) data consists of quantifying the statistical dependencies among time series describing the activity of different neural sources from the magnetic field recorded outside the scalp. This problem can be addressed by utilizing connectivity measures whose computation in the frequency domain often relies on the evaluation of the cross-power spectrum of the neural time series estimated by solving the MEG inverse problem. Recent studies have focused on the optimal determination of the cross-power spectrum in the framework of regularization theory for ill-posed inverse problems, providing indications that, rather surprisingly, the regularization process that leads to the optimal estimate of the neural activity does not lead to the optimal estimate of the corresponding functional connectivity. Along these lines, the present paper utilizes synthetic time series simulating the neural activity recorded by an MEG device to show that the regularization of the cross-power spectrum is significantly correlated with the signal-to-noise ratio of the measurements and that, as a consequence, this regularization correspondingly depends on the spectral complexity of the neural activity.

1. Introduction

Magnetoencephalography (MEG) provides high temporal resolution measurements of the magnetic field associated to neural currents. The MEG device relies on superconducting sensors, named SQUIDs, organized in a helmet array close and around the scalp. MEG experimental time series can be used essentially to address two neuroscientific problems, whose solution requires both an accurate mathematical modelization based on Maxwell’s equations, and the numerical reduction of such formal models [1].

The first problem is concerned with the dynamical ill-posed inverse problem of estimating parameters associated with the neural sources inducing the magnetic field signal [2,3,4,5,6,7,8]. The second problem is concerned with the quantification of the interactions among neural sources located in different cortical areas and intertwined by means of either anatomical or functional connectivity [9,10,11,12,13,14].

In particular, the connectivity problem can be addressed by either computing proper connectivity metrics directly from the experimental time series provided by the MEG sensors or searching for connections in the source space, i.e., among the neural time series estimated as solutions of the inversion process. This second approach has the advantages of reducing the impact of volume conduction and providing results that can be more easily interpreted in the framework of neuroscientific models [15,16,17]. Several approaches for identifying connectivity paths rely on physiological models assuming that the functional communication between different brain areas is regulated by the synchronization of their activity at specific temporal frequencies [18,19]. This implies that, for these models, the frequency domain represents the natural computational framework where to perform the connectivity analysis. This is the reason why, in the present paper, we focus on the analysis of the cross-power spectrum, which is the mathematical quantity of reference for the computation of most frequency-domain connectivity measures [11,20,21]. From an operational viewpoint, the computation of the cross-power spectrum in the source space typically relies on a two-step procedure: first the neural activity is estimated by applying a regularized inversion method on the recorded time series and then the cross-power spectrum is computed from the Fourier transform of the estimated neural time series [22].

This paper investigates how to optimally select the regularization parameter in the inversion procedure in order to obtain the best possible estimate of the neural cross-power spectrum. In fact, we consider the Tikhonov method (better known as Minimum Norm Estimation (MNE) in the MEG world [2]) as it is one of the most commonly used inverse methods in connectivity studies [23,24,25]; we study the interplay between the regularization parameter providing the reconstructed neural time series minimizing the relative error in ${\ell }_2$-norm, and the one that allows the optimal estimate of the cross-power spectrum according to the normalized Frobenius norm. The conceptual motivation of this problem is illustrated in Figure 1, which tentatively sketches the result of recent investigations in MEG-based connectivity research, i.e., that the regularization parameter leading to the optimal estimate of the neural activity may not lead to the optimal estimate of the cross-power spectrum and vice versa. In fact, in [26] the authors used numerical simulations to compare the parameter that provides the best estimate of the power spectrum with the one that provides the best estimate of coherence and showed that the latter is in general two orders of magnitude smaller than the former. More recently, Vallarino et al. [27] addressed an analogous problem via analytical computations, considering a simplified model. Specifically, under the assumption that the neural time series are realizations of white Gaussian processes, the authors proved that the parameter providing the best neural activity estimate is more than twice as large as the one providing the best estimate of the cross-power spectrum.

The present paper focuses on an analysis of the impact of spectral complexity of the actual neural signal on the value of the two regularization parameters. Specifically, we simulate synthetic MEG signals and discuss how the optimal parameter for the reconstruction of the cross-power spectrum depends on its signal-to-noise ratio and how this latter quantity is related to the spectral richness of the neural sources. To this aim, we considered a simulation setting in which the signal is modeled as a multivariate autoregressive process.

This paper is structured as follows. Section 2 introduces the problem in a formal way. Section 3 describes how the synthetic data are simulated and analyzed. Section 4 presents the results of the analysis. Our conclusions are offered in Section 5.

2. Definition of the Problem

2.1. Forward Model

Let $\mathbf{X}\left(t\right)={(X_1\left(t\right),\dots ,X_N\left(t\right))}^{\top }\in R^N$ be a multivariate stationary stochastic process whose realizations $\mathbf{x}\left(t\right)$ can not be observed and let $\mathbf{Y}\left(t\right)={(Y_1\left(t\right),\dots ,Y_M\left(t\right))}^{\top }\in R^M$ be the process whose realizations $\mathbf{y}\left(t\right)$ are used to infer information on $\mathbf{x}\left(t\right)$. Let $\mathbf{Y}\left(t\right)$ and $\mathbf{X}\left(t\right)$ be related by the following equation

$$\mathbf{Y}\left(t\right)=\mathbf{G}\mathbf{X}\left(t\right)+\mathbf{N}\left(t\right),$$

(1)

where $\mathbf{G}\in R^{M\times N}$ is the forward matrix and $\mathbf{N}\left(t\right)={(N_1\left(t\right),\dots ,N_M\left(t\right))}^{\top }\in R^M$ is the measurement noise, that is here assumed to be a white Gaussian process with zero mean and covariance matrix ${\alpha }^2\mathbf{I}$, i.e., $\mathbf{N}\left(t\right)\sim \mathcal{N}(0,{\alpha }^2\mathbf{I})$, independent from $\mathbf{X}\left(t\right)$.

2.2. Cross-Power Spectrum

We are interested in reconstructing the cross-power spectrum of $\mathbf{X}\left(t\right)$, which describes the statistical dependencies between each pair of time series ${(X_j\left(t\right),X_k\left(t\right))}_{j,k\in \{1,\dots ,N\}}$. The cross-power spectrum is a one parameter family of $N\times N$ matrices ${\mathbf{S}}^{\mathbf{X}}\left(f\right)$, whose $(j,k)-th$ element is defined as

$$S_{j,k}^{\mathbf{X}}\left(f\right)=\underset{T\to \infty }{lim}\frac{1}{T}E\left[{\widehat{X}}_j(f,T){\widehat{X}}_k{(f,T)}^H\right],$$

(2)

where ${\widehat{X}}_j(f,T)$ is the Fourier transform of $X_j\left(t\right)$ over the interval $[0,T]$, defined as

$${\widehat{X}}_j(f,T)={\int }_0^T{X}_j\left(t\right)e^{-2\pi ift}dt$$

(3)

and $X^H$ is the Hermitian transpose of X [28].

Given a realization $\mathbf{x}\left(t\right)$ of the process $\mathbf{X}\left(t\right)$, the cross-power spectrum ${\mathbf{S}}^{\mathbf{X}}\left(f\right)$ can be estimated via the Welch method [29], which consists in partitioning the data in P overlapping segments multiplied by a window function, ${\left\{w\left(t\right){\mathbf{x}}^p\left(t\right)\right\}}_{p=1}^P$, computing their discrete Fourier transform ${\widehat{\mathbf{x}}}^p\left(f\right)=\frac{1}{L}{\sum }_{t=0}^{L-1}{\mathbf{x}}^p\left(t\right)w\left(t\right)e^{\frac{-2\pi itf}{L}}$ and averaging:

$${\mathbf{S}}^{\mathbf{x}}\left(f\right)=\frac{L}{PW}\sum _{p=1}^P{\widehat{\mathbf{x}}}^p\left(f\right){\widehat{\mathbf{x}}}^p{\left(f\right)}^H,f=0,\dots ,L-1,$$

(4)

where L is the length of each segment and $W=\frac{1}{L}{\sum }_{t=0}^{L-1}w{\left(t\right)}^2$.

It is often the case that the data reach a high dimension, and visual inspection of the cross-power spectrum is not doable. In such cases a metric that describes the spectral properties of the signals would be useful. Here we use the spectral complexity coefficient, defined as follows.

Definition 1.

Given a realization $\mathbf{x}\left(t\right)$ of the process $\mathbf{X}\left(t\right)$, and the corresponding cross-power spectrum ${\mathbf{S}}^{\mathbf{x}}\left(f\right)$, we define the spectral complexity coefficient as the average of the elements of the upper triangular part of the matrix obtained by computing the squared ${\ell }_2-$norm over the frequencies of ${\mathbf{S}}_{j,k}^{\mathbf{x}}\left(f\right)$, $j,k=1,\dots ,N$, that is

$$c=\frac{2}{N(N+1)}\sum _{j=1}^N\sum _{k=j}^N\sum _f{\left|{\mathbf{S}}_{j,k}^{\mathbf{x}}\left(f\right)\right|}^2.$$

(5)

The spectral complexity coefficient assumes small values if the elements of the cross-power spectrum are flat, that is when time series do not present any periodic trend and no dependencies among the pairs of time series are present. On the contrary, it assumes large values if the elements of the cross-power spectrum are peaked, that is when time series present periodic trends and complex relations among them. Finally, we observe that in Definition 1 only the elements on the upper triangular part of ${\mathbf{S}}^{\mathbf{x}}\left(f\right)$ are considered because ${\mathbf{S}}^{\mathbf{x}}\left(f\right)$ is Hermitian.

2.3. Two-Step Approach for Cross-Power Spectrum Estimation

Let us now consider a realization of Equation (1). Further than an estimate of the hidden data $\mathbf{x}\left(t\right)$, an estimate of the cross-power spectrum can be obtained from $\mathbf{y}\left(t\right)$. Such estimate can be achieved through a two-step process [22]:

• First, a regularized estimate ${\mathbf{x}}_{\lambda }\left(t\right)$ of $\mathbf{x}\left(t\right)$ is obtained by solving the inverse problem associated to Equation (1). Here we consider the Tikhonov regularized solution [30] of the problem which is defined as

  $${\mathbf{x}}_{\lambda }\left(t\right)=\arg \underset{\mathbf{x}\left(t\right)}{min}\left\{{∥\mathbf{Gx}\left(t\right)-\mathbf{y}\left(t\right)∥}_2^2+\lambda {∥\mathbf{x}\left(t\right)∥}_2^2\right\};$$

  (6)

  where $\lambda$ is a proper regularization parameter and ${\parallel ·\parallel }_2$ is the ${\ell }_2$-norm.
• Then, the corresponding estimate of the cross-power spectrum ${\mathbf{S}}^{{\mathbf{x}}_{\lambda }}\left(f\right)$ is computed from the reconstructed time series using the Welch method, as described in the previous section.

Remark 1.

In many applied fields, Tikhonov regularization with an ${\ell }_2$ penalty term has been outdated by more modern techniques that use sparsity-inducing penalization terms such as ${\ell }_1$ or ${\ell }_p$ with $0<p<1$. Indeed, also in the M/EEG literature there has been considerable effort in developing ${\ell }_1$ solutions [31,32], and mixed norm solutions [33]; both these approaches have proved to provide superior performances in terms of localization of neural activity. However, these newer methods are seldom used in connectivity studies, for good reasons: ${\ell }_1$ solutions computed independently at each time point produce extremely jittering reconstructions, resulting in highly sparse time courses that are not suitable for computing connectivity metrics. Mixed norms, that have been developed precisely to overcome this jittering problem, are computationally very expensive, and this actually prevents their use with the large datasets typically involved in connectivity studies.

When applying the described two-step process, the regularization parameter $\lambda$ in Equation (6) has to be set for the computation of ${\mathbf{x}}_{\lambda }\left(t\right)$. Thus, the problem naturally arises of the choice of such parameter, which can be set in order to optimally reconstruct either ${\mathbf{x}}_{\lambda }\left(t\right)$ or ${\mathbf{S}}^{{\mathbf{x}}_{\lambda }}\left(f\right)$. We define optimality through the minimization of the normalized norm of the discrepancy between the true and the reconstructed time series and cross-power spectra as follows.

Definition 2.

Given the regularized solution (6) and the cross-power spectrum (4), we define the optimal regularization parameter for the reconstruction of $\mathbf{x}\left(t\right)$ as

$${\lambda }_{\mathbf{x}}^{*}=arg\underset{\lambda }{min}{\epsilon }_{\mathbf{x}}\left(\lambda \right)\mathrm{with}{\epsilon }_{\mathbf{x}}\left(\lambda \right)=\frac{{\sum }_t{∥{\mathbf{x}}_{\lambda }\left(t\right)-\mathbf{x}\left(t\right)∥}_2^2}{{\sum }_t{∥{\mathbf{x}}_{\lambda }\left(t\right)∥}_2^2+{\sum }_t{∥\mathbf{x}\left(t\right)∥}_2^2};$$

(7)

and the optimal parameter for the reconstruction of ${\mathbf{S}}^{\mathbf{x}}\left(f\right)$ as

$${\lambda }_{\mathbf{S}}^{*}=arg\underset{\lambda }{min}{\epsilon }_{\mathbf{S}}\left(\lambda \right)with{\epsilon }_{\mathbf{S}}\left(\lambda \right)=\frac{{\sum }_f{∥{\mathbf{S}}^{{\mathbf{x}}_{\lambda }}\left(f\right)-{\mathbf{S}}^{\mathbf{x}}\left(f\right)∥}_F^2}{{\sum }_f{∥{\mathbf{S}}^{{\mathbf{x}}_{\lambda }}\left(f\right)∥}_F^2+{\sum }_f{∥{\mathbf{S}}^{\mathbf{x}}\left(f\right)∥}_F^2};$$

(8)

where ${\parallel ·\parallel }_F$ is the Frobenius norm; ${\epsilon }_{\mathbf{x}}\left(\lambda \right)$ and ${\epsilon }_{\mathbf{S}}\left(\lambda \right)$ will be called reconstruction errors.

The reconstruction errors range from 0 to 1 and penalize both a too small and a too large value of $\lambda$. In fact, they assume their maximum value when either $\lambda$ is very high and thus ${\mathbf{x}}_{\lambda }\left(t\right)$ is negligible with respect to $\mathbf{x}\left(t\right)$, or when $\lambda$ is too small and thus, vice versa, $\mathbf{x}\left(t\right)$ is negligible with respect to ${\mathbf{x}}_{\lambda }\left(t\right)$. This definition may appear overly complex compared to, e.g., a mere ${\ell }_2$-norm of the difference; however, in the presence of sparse data where only a few time series are non-zero, the simple ${\ell }_2$-norm would prefer a very high regularization parameter in order to minimize the error on the null time series, at the expense of the error on the non-zero ones; our definition aims to cope with this limitation of the ${\ell }_2$-norm. A similar definition has been introduced in [34].

In experimental contexts, where $\mathbf{x}\left(t\right)$ is not known, the choice of the optimal regularization parameter is crucial. This matter is widely discussed in literature [35,36,37,38], and many criteria have been proposed. Such criteria apply to Equation (1) and can be used to set the regularization parameter ${\lambda }_{\mathbf{x}}$. A possibility is to set the regularization parameter as a function of the signal-to-noise ratio (SNR), which describes the level of the desired signal with respect to that of the measurement noise; for Equation (1) the SNR is defined as follows.

Definition 3.

Consider the linear model (1). We define the signal-to-noise ratio of $\mathbf{X}\left(t\right)$ related to such model as

$${\mathrm{SNR}}^{\mathbf{X}}=10{log}_{10}\left(\frac{{\sum }_t{∥\mathbf{G}\mathbf{X}\left(t\right)∥}_2^2}{{\sum }_t{∥\mathbf{N}\left(t\right)∥}_2^2}\right).$$

(9)

To the best of our knowledge, the choice of the optimal regularization parameter for the reconstruction of the cross-power spectrum has never been related to the signal-to-noise ratio. This relation is presented in Section 4; however we first need to relate the cross-power spectrum of the unknown ${\mathbf{S}}^{\mathbf{X}}\left(f\right)$ with that of the data ${\mathbf{S}}^{\mathbf{Y}}\left(f\right)$.

By computing the cross-power spectrum of both sides of Equation (1) and from the linearity of the Fourier transform it follows that

$${\mathbf{S}}^{\mathbf{Y}}\left(f\right)=\mathbf{G}{\mathbf{S}}^{\mathbf{X}}\left(f\right){\mathbf{G}}^{\top }+{\mathbf{S}}^{\mathbf{N}}\left(f\right),$$

(10)

where the mixed terms ${\mathbf{S}}^{\mathbf{XN}}\left(f\right)$ and ${\mathbf{S}}^{\mathbf{NX}}\left(f\right)$ are negligible thanks to the independence between $\mathbf{X}\left(t\right)$ and $\mathbf{N}\left(t\right)$. Just like for Equation (1), we can define the signal-to-noise ratio for Equation (10) as follows.

Definition 4.

Consider the linear model (10). We define the signal-to-noise ratio of ${\mathbf{S}}^{\mathbf{X}}\left(f\right)$ related to such model as

$${\mathrm{SNR}}^{\mathbf{S}}=10{log}_{10}\left(\frac{{\sum }_f{∥\mathbf{G}{\mathbf{S}}^{\mathbf{X}}\left(f\right){\mathbf{G}}^{\top }∥}_F^2}{{\sum }_f{∥{\mathbf{S}}^{\mathbf{N}}\left(f\right)∥}_F^2}\right).$$

(11)

This definition is in line with the definition of ${\mathrm{SNR}}^{\mathbf{X}}$ for the signal, the main difference being in the use of the Frobenius norm rather than the ${\ell }_2$-norm, motivated by the fact that we are working with matrices rather than vectors.

2.4. Multivariate Autoregressive Models

To model the statistical relationships between the different components of the stochastic process $\mathbf{X}\left(t\right)$, in this work the latter is assumed to follow a stable multivariate autoregressive model of order P [39].

Definition 5.

A zero-mean stochastic process $\mathbf{X}\left(t\right)\in {\mathbb{R}}^N$ is said to follow a multivariate autoregressive (MVAR) model of order P if

$$\mathbf{X}\left(t\right)=\sum _{k=1}^P\mathbf{A}\left(k\right)\mathbf{X}(t-k)+\mathit{\epsilon }\left(t\right)\forall t,$$

(12)

where $\mathbf{A}\left(k\right)\in {\mathbb{R}}^{N\times N}$ are fixed coefficient matrices, and $\mathit{\epsilon }\left(t\right)\in {\mathbb{R}}^N$ is a white Gaussian noise process.

Moreover, the MVAR model described by Equation (12) is said to be stable if

$$det(\mathbf{I}-\sum _{k=1}^P\mathbf{A}\left(k\right)z^k)\ne 0\forall z\in \mathbb{C}s.t.|z|\le 1,$$

(13)

where $\mathbf{I}$ is the identity matrix of size N.

Remark 2.

From Equation (12) it can be easily seen that the process $\mathbf{X}\left(t\right)$ is uniquely determined by the process $\mathit{\epsilon }\left(t\right)$ and by the first P time points, $\mathbf{X}\left(0\right)$, ⋯, $\mathbf{X}(P-1)$. Indeed, consider for example an MVAR model of order 1 (a similar proof holds for the general case $P>1$ and can be found in [39]); then, for each time point t

$$\begin{array}{cc}\hfill \mathbf{X}\left(t\right)& =\mathbf{A}\left(1\right)\mathbf{X}(t-1)+\mathit{\epsilon }\left(t\right)\hfill \\ & =\mathbf{A}{\left(1\right)}^2\mathbf{X}(t-2)+\mathbf{A}\left(1\right)\mathit{\epsilon }(t-1)+\mathit{\epsilon }\left(t\right)\hfill \\ & =\mathbf{A}{\left(1\right)}^t\mathbf{X}\left(0\right)+\sum _{k=0}^{t-1}\mathbf{A}{\left(1\right)}^k\mathit{\epsilon }(t-k).\hfill \end{array}$$

Such a model satisfies the stability condition defined in Equation (13) if all the eigenvalues of the coefficient matrix $\mathbf{A}\left(1\right)$ have a modulus of less then one, the condition that guarantees the sequence of exponential matrices ${\left\{\mathbf{A}{\left(1\right)}^k\right\}}_k$ is absolutely summable.

According to Equation (12), if the process $\mathbf{X}\left(t\right)$ follows an MVAR model, then at each time point the value of $\mathbf{X}\left(t\right)$ can be derived as a weighted sum of the values of the process at the previous P time points, $\mathbf{X}(t-1)$, ..., $\mathbf{X}(t-P)$, plus a random perturbation $\mathit{\epsilon }\left(t\right)$. In particular, the $(i,j)$-th elements of the coefficient matrices, $a_{ij}\left(1\right)$, ..., $a_{ij}\left(P\right)$, describe how the value of the i-th component of the process depends on the past of the j-th component. Different connectivity patterns, with various levels of complexity, can thus be obtained by tuning the off-diagonal values of the coefficient matrices. Due to their flexibility and simplicity, MVAR models have been used by various authors in the framework of MEG functional connectivity estimation as a benchmark for testing and comparing different connectivity metrics [21,23,34,40,41,42,43]. Other models have been proposed to simulate different connectivity patterns, such as coherent sinusoidal time series [26], neural mass models [44,45], or Kuramoto models [46,47]. However a comprehensive comparison of all possible generative models is beyond the scope of this work.

3. Generation and Analysis Pipeline of the MEG Simulated Data

In this section we describe the numerical simulation that led to the main results of our study. First we introduce the continuous MEG forward problem and its discretized version, then we describe how we generated the data and, finally, we describe the inverse model and how we numerically computed the optimal regularization parameters.

3.1. MEG Forward Model

The MEG forward problem aims at computing the magnetic field produced outside the head by an electric current that flows inside the brain. The quasi static approximation of Maxwell’s equations provides the local relationship between the recorded magnetic field and the neural currents [1,48,49]. The two equations that are of interest here read as

$$\begin{array}{cc}& \nabla \times \mathbf{E}(\mathbf{r},t)=0\end{array}$$

(14)

$$\begin{array}{cc}& \nabla \times \mathbf{B}(\mathbf{r},t)={\mu }_0\mathbf{J}(\mathbf{r},t);\end{array}$$

(15)

where $\mathbf{E}(\mathbf{r},t)$ and $\mathbf{B}(\mathbf{r},t)$ are the electric and magnetic fields at location $\mathbf{r}$ and time t, ${\mu }_0$ is the magnetic permeability in vacuum and $\mathbf{J}(\mathbf{r},t)$ is the total electric current that flows inside the brain. The latter is the sum of two contributions

$$\mathbf{J}(\mathbf{r},t)={\mathbf{J}}^p(\mathbf{r},t)+{\mathbf{J}}^v(\mathbf{r},t),$$

(16)

${\mathbf{J}}^p(\mathbf{r},t)$ being the primary current directly related to the brain activity, while ${\mathbf{J}}^v(\mathbf{r},t)=-\sigma \left(\mathbf{r}\right)\nabla V(\mathbf{r},t)$ is the induced volume current due the non-null conductivity $\sigma \left(\mathbf{r}\right)$ of the brain, $V(\mathbf{r},t)$ being the electric scalar potential.

The manipulation of Maxwell’s Equations leads to the Biot–Savart Equation

$$\mathbf{B}(\mathbf{r},t)={\mathbf{B}}_0(\mathbf{r},t)-\frac{{\mu }_0}{4\pi }{\int }_{\Omega }\sigma \left({\mathbf{r}}^{\prime }\right){\nabla }^{\prime }V({\mathbf{r}}^{\prime },t)\times \frac{\mathbf{r}-{\mathbf{r}}^{\prime }}{|r-r^{\prime }{|}^3}\mathrm{d}{v}^{\prime },$$

(17)

where $\Omega$ is the volume occupied by the brain, the first term ${\mathbf{B}}_0(\mathbf{r},t)=\frac{{\mu }_0}{4\pi }{\int }_{\Omega }\mathbf{J}({\mathbf{r}}^{\prime },t)\frac{\mathbf{r}-{\mathbf{r}}^{\prime }}{|r-r^{\prime }{|}^3}\mathrm{d}{v}^{\prime }$ is the magnetic field induced by the primary current, whereas the second term is related to the volume current.

Solving the forward problem requires the computation of these two contributions knowing the primary current. Although for the first one straightforward numerical integration is feasible, for the second one it is common to model the head as the union of nested homogeneous volumes ${\left\{{\Omega }_j\right\}}_{j=1,\dots ,J}$ and to replace volume integration with surface integration. In this way, the Biot–Savart equation becomes

$$\mathbf{B}(\mathbf{r},t)={\mathbf{B}}_0(\mathbf{r},t)+\frac{{\mu }_0}{4\pi }\sum _{i,j}({\sigma }_i-{\sigma }_j){\int }_{\partial {\Omega }_{i,j}}V({\mathbf{r}}^{\prime },t)\frac{\mathbf{r}-{\mathbf{r}}^{\prime }}{|r-r^{\prime }{|}^3}\times {\mathbf{n}}_{i,j}\left({\mathbf{r}}^{\prime }\right)\mathrm{d}{s}^{\prime },$$

(18)

where $\partial {\Omega }_{i,j}$ is the contact surface between regions ${\Omega }_i$ and ${\Omega }_j$ and ${\mathbf{n}}_{i,j}\left({\mathbf{r}}^{\prime }\right)$ is the unit vector normal to the surface $\partial {\Omega }_{i,j}$ at ${\mathbf{r}}^{\prime }$ from region i to region j.

The forward problem can now be solved by computing the second term at the right hand side of Equation (18) after having computed $V(\mathbf{r},t)$ by solving the equation

$$\nabla ·{\mathbf{J}}^p(\mathbf{r},t)-\nabla ·\left(\sigma \left(\mathbf{r}\right)\nabla V(\mathbf{r},t)\right)=0,$$

(19)

which follows from Equation (15) by applying the divergence.

For further details on the MEG forward problem, we refer the reader to [1].

3.2. The Leadfield Matrix

Experimental contexts require the discretization of the forward problem. This involves a discretization of both the volume occupied by the brain and the volume outside the head.

When using a distributed model for the primary current ${\mathbf{J}}^p$, the brain volume is uniformly divided in N small parcels. If N is sufficiently big and thus each parcel has a sufficient small area, the activity in each brain parcel is approximated by a point-like source, henceforth denoted as dipole. From a mathematical point of view, each dipole is a vector whose strength and direction represent the intensity and orientation of the primary current in the corresponding brain area [1].

As for the volume outside the brain, it is natural to discretize it in correspondence of the MEG sensors. Let us denote the measured magnetic field as $\mathbf{y}\left(t\right)=(y_1\left(t\right),\dots ,y_M\left(t\right))$. Now, observing that the magnetic field $\mathbf{B}$ depends linearly on the primary current ${\mathbf{J}}^p$, the magnetic field in correspondence of the sensors of the instrument is

$$\mathbf{y}\left(t\right)=\sum _{k=1}^{N}G\left({\mathbf{r}}_k\right){\mathbf{q}}_k\left(t\right)+\mathbf{n}\left(t\right),$$

(20)

where ${\mathbf{r}}_k$, $k=1,\dots ,N$, is the location of the k-th brain parcel, $G\left({\mathbf{r}}_k\right)\in R^{M\times 3}$ is the corresponding leadfield matrix, and ${\left\{{\mathbf{q}}_k\left(t\right)\right\}}_{k=1,\dots ,N}$ are the electric current intensities along the three orthogonal direction of the N dipoles within the brain at time t and $\mathbf{n}\left(t\right)$ is the measurement noise. The l-th column of $G\left({\mathbf{r}}_k\right)$ contains the measurement at a sensor level when a unit current dipole is placed at location ${\mathbf{r}}_k$ and oriented along the l-th orthogonal direction.

In this work, we assume dipoles to be located only on the brain cortical mantle and their orientation to be normal to the local cortical surface [50]. In this case, the electric current intensities are scalar quantities (we refer to them as ${\left\{q_k\right\}}_{k=1,\dots ,N}$) and the leadfield matrices are column vectors (we refer to them as ${\left\{G_k\right\}}_{k=1,\dots ,N}$).

Let us define

$$\mathbf{x}\left(t\right):=(q_1\left(t\right),\dots ,q_N\left(t\right))$$

(21)

and

$$\mathbf{G}:=[G_1,\dots ,G_N]\in R^{M\times N};$$

(22)

reassembling Equations (21) and (22) in to Equation (20), we get

$$\mathbf{y}\left(t\right)=\mathbf{G}\mathbf{x}\left(t\right)+\mathbf{n}\left(t\right),$$

(23)

which can be interpreted as a realization of Equation (1). From now on we refer to $\mathbf{G}$ as to the leadfield matrix.

For the simulation presented in this work, we used the leadfield matrix available in the sample dataset of MNE Python [51]. We selected magnetometers and set a fixed orientation. For computational reasons, the available source space, containing 1884 sources, was uniformly down-sampled to obtain 274 sources. Thus, our model has $M=102$ sensors and $N=274$ dipole sources.

3.3. Data Generation

We simulated $N_{mod}=10$ pairs of active sources, ${(z_1\left(t\right),z_2\left(t\right))}^{\top }$, with unidirectional coupling from the first to the second; their time series follow a multivariate autoregressive (MVAR) model of order $P=5$ [39,43]

$$\left(\begin{array}{c}{z}_1\left(t\right)\\ z_2\left(t\right)\end{array}\right)=\sum _{k=1}^P\left(\begin{array}{cc}{a}_{1,1}\left(k\right)& 0\\ a_{2,1}\left(k\right)& a_{2,2}\left(k\right)\end{array}\right)\left(\begin{array}{c}{z}_1(t-k)\\ z_2(t-k)\end{array}\right)+\left(\begin{array}{c}{\mathit{\epsilon }}_1\left(t\right)\\ {\mathit{\epsilon }}_2\left(t\right)\end{array}\right),t=P,\dots ,T.$$

(24)

The non-zero elements $a_{i,j}\left(k\right)$ of the coefficient matrices were drawn from a normal distribution of zero mean and standard deviation $\gamma$, and T = 10,000. We retained only coefficient matrices providing (i) a stable MVAR model [39] and (ii) pairs of signals ${(z_1\left(t\right),z_2\left(t\right))}^{\top }$ such that the ${\ell }_2$-norm of the strongest one was less than three times the ${\ell }_2$-norm of the weakest one. In order to obtain time series with different spectral complexity coefficients we set $\gamma$ to $N_{mod}$ different values randomly drawn in the interval $[0.1,1]$. The values of the spectral complexity coefficient of the $N_{mod}$ simulated time series are reported in

Table 1. The table reports the values of the spectral complexity coefficients, $c_j$, associated to each simulated multivariate autoregressive (MVAR) model, $m_j$, $j=1,\dots ,N_{mod}$.

Model	1	2	3	4	5	6	7	8	9	10
Spectral complexity coefficient	1.41	1.96	2.14	3.10	3.44	4.17	4.64	5.67	6.69	8.67

. Finally, the resulting time series ${(z_1\left(t\right),z_2\left(t\right))}^{\top }$ were normalized by the mean of their standard deviations over time, so that pairs of time series drawn from different models had similar magnitude. Figure 2 shows a sample of the the cross-power spectra among the simulated pairs of time series. The figure shows that for increasing values of the spectral complexity coefficient the cross-power spectrum of the corresponding time series becomes more peaked. Each pair of simulated time series was then assigned to $N_{loc}=20$ pairs of point like sources randomly chosen in the source space, so that the ratio of the norms of the corresponding columns of the leadfield matrix was close to one, i.e., they had similar intensity at a sensor level, and their distance was grater than $7\mathrm{cm}$. The remaining $N-2$ sources were set to have null activity.

Source space activity was then projected to sensor level by multiplying the simulated source activity by the leadfield matrix and white Gaussian noise was added to obtain $N_{snr}=6$ levels of ${\mathrm{SNR}}^{\mathbf{X}}$ evenly spaced in the interval $[-20\mathrm{dB},5\mathrm{dB}]$.

Summarizing, we generated $N_{mod}·N_{loc}·N_{snr}=1200$ different sensor level configurations. The green box in Figure 3 shows a visual representation of the simulation pipeline.

3.4. Inverse Model

Source space time series were reconstructed using the Tikhonov method, also known as the minimum norm estimate (MNE) [2] within the MEG community. For each combination of source time series, source locations, and ${\mathrm{SNR}}^{\mathbf{X}}$ level, we computed the optimal regularization parameters ${\lambda }_{\mathbf{x}}^{*}$ and ${\lambda }_{\mathbf{S}}^{*}$ by minimizing the reconstruction errors ${\epsilon }_{\mathbf{x}}\left(\lambda \right)$ and ${\epsilon }_{\mathbf{S}}\left(\lambda \right)$, defined in Definition 2. The minimization procedure was achieved by using the Matlab built in function fminsearch that implements an iterative procedure based on the simplex method developed by Lagarias and colleagues [52]. In more detail, ${\lambda }_{\mathbf{x}}^{*}$ and ${\lambda }_{\mathbf{S}}^{*}$ have been obtained by applying such procedure to ${\epsilon }_{\mathbf{x}}\left(\lambda \right)$ and ${\epsilon }_{\mathbf{S}}\left(\lambda \right)$, respectively; in both cases the starting point of the simplex method was set equal to ${10}^{\left(-\frac{{\mathrm{SNR}}^{\mathbf{X}}}{10}\right)}$, which corresponds to the optimal value of ${\lambda }_{\mathbf{x}}$ in the case of white Gaussian signals [27]. The blue box in Figure 3 describes the inverse procedure to obtain an estimate of the cross-power spectrum and stresses the role of the regularization parameter in the two-step process.

4. Results

In this section we illustrate the results of our analysis. We begin with the description of the analytical dependence between ${\mathrm{SNR}}^{\mathbf{X}}$ and ${\mathrm{SNR}}^{\mathbf{S}}$, then we highlight how the optimal parameter for the reconstruction of the cross-power spectrum depends on ${\mathrm{SNR}}^{\mathbf{S}}$ and how this implies that the spectral complexity of the signal is behind such dependence. Finally we show how the reconstruction error ${\epsilon }_{\mathbf{S}}\left(\lambda \right)$ behaves for different values of the regularization parameters. As a byproduct, this analysis also confirms the results of Vallarino et al. [27] in the case of a more complex setting.

4.1. Analytical Relation between ${\mathrm{SNR}}^{\mathbf{X}}$ and ${\mathrm{SNR}}^{\mathbf{S}}$

From Equations (9) and (11) and reminding that $\mathbf{N}\left(t\right)\sim \mathcal{N}(0,{\alpha }^2\mathbf{I})$ it follows that

$${\mathrm{SNR}}^{\mathbf{X}}=10{log}_{10}\left(\frac{{\sum }_t{∥\mathbf{G}\mathbf{X}\left(t\right)∥}_2^2}{MT{\alpha }^2}\right);$$

(25)

and

$${\mathrm{SNR}}^{\mathbf{S}}=10{log}_{10}\left(\frac{{\sum }_f{∥\mathbf{G}{\mathbf{S}}^{\mathbf{X}}\left(f\right){\mathbf{G}}^{\top }∥}_F^2}{M{N}_f{\alpha }^2}\right),$$

(26)

where T is the number of time points and $N_f$ is the number of frequencies used to compute the cross-power spectrum. Observe that to derive Equation (26) we used the fact that the cross-power spectrum of a white noise Gaussian process of zero mean and covariance matrix ${\alpha }^2\mathbf{I}$ is ${\mathbf{S}}^{\mathbf{N}}\left(f\right)={\alpha }^2\mathbf{I}$.

By isolating ${\alpha }^2$ from Equation (25) and substituting in Equation (26) we obtain

$${\mathrm{SNR}}^{\mathbf{S}}=10{log}_{10}\left(\frac{T^{2}M{\sum }_f{∥\mathbf{G}{\mathbf{S}}^{\mathbf{X}}\left(f\right){\mathbf{G}}^{\top }∥}_F^2}{N_f{\sum }_t{∥\mathbf{G}\mathbf{X}\left(t\right)∥}_2^4}\right)+2{\mathrm{SNR}}^{\mathbf{X}}.$$

(27)

Equation (27) relates the signal-to-noise ratio of $\mathbf{X}\left(t\right)$ with that of ${\mathbf{S}}^{\mathbf{X}}\left(f\right)$. It shows that, for same levels of ${\mathrm{SNR}}^{\mathbf{X}}$, ${\mathrm{SNR}}^{\mathbf{S}}$ changes with the spectral complexity coefficient of the signals. In fact, the higher the spectral complexity coefficient, the higher the quantity ${∥\mathbf{G}{\mathbf{S}}^{\mathbf{X}}\left(f\right){\mathbf{G}}^{\top }∥}_F^2$. Intuitively, this happens because when the signal has a higher spectral complexity coefficient its cross-power spectrum is more peaked and thus it is stronger over the cross-power spectrum of the noise with respect to a signal with a lower spectral complexity coefficient.

4.2. Dependence of ${\lambda }_{\mathbf{S}}^{*}$ on ${\mathrm{SNR}}^{\mathbf{S}}$

As described in Section 3 we simulated several sensor level configurations, based on different combinations of spectral complexity coefficients, source locations, and ${\mathrm{SNR}}^{\mathbf{X}}$ levels. For each configuration, we collected the two optimal parameters ${\lambda }_{\mathbf{x}}^{*}$ and ${\lambda }_{\mathbf{S}}^{*}$ and we investigated their dependence on the signal-to-noise-ratio. In accordance with classical results from inverse theory [36], we found that ${\lambda }_{\mathbf{x}}^{*}$ depends on the signal-to-noise ratio. What is novel here is the relation between ${\lambda }_{\mathbf{S}}^{*}$ and both ${\mathrm{SNR}}^{\mathbf{X}}$ and ${\mathrm{SNR}}^{\mathbf{S}}$. Indeed, for increasing ${\mathrm{SNR}}^{\mathbf{X}}$, less regularization is needed, but such dependence varies with the MVAR models. On the other side, the dependence of ${\lambda }_{\mathbf{S}}^{*}$ on ${\mathrm{SNR}}^{\mathbf{S}}$ is neater and does not depend on the models. Figure 4 shows this result; on the left the regularization parameters for the cross-power spectrum reconstruction versus ${\mathrm{SNR}}^{\mathbf{X}}$ are shown, while on the right the same parameters are shown with respect to ${\mathrm{SNR}}^{\mathbf{S}}$. For the ease of presentation the figure shows the parameters related to one source location; while on the left lines corresponding to different MVAR models have different heights, on the right they overlap.

4.3. ${\lambda }_{\mathbf{S}}^{*}<{\lambda }_{\mathbf{x}}^{*}$ and Dependency from the Spectral Complexity

We also investigated the relation between the two optimal regularization parameters. Figure 5 shows the ratio $\frac{{\lambda }_{\mathbf{S}}^{*}}{{\lambda }_{\mathbf{X}}^{*}}$ versus ${\mathrm{SNR}}^{\mathbf{X}}$ for the simulated MVAR models. The ratio between the two parameters is always smaller than $\frac{1}{2}$, meaning that ${\lambda }_{\mathbf{S}}^{*}<\frac{1}{2}{\lambda }_{\mathbf{x}}^{*}$, as it was analytically proved in a simplified case in [27]. Further to this, the figure shows that for increasing spectral complexity coefficients this ratio gets smaller. This latter result is directly related to Equation (27). In fact, for same levels of ${\mathrm{SNR}}^{\mathbf{X}}$, signals with higher spectral complexity have higher ${\mathrm{SNR}}^{\mathbf{S}}$ and, thus, need less regularization.

4.4. The Reconstruction Errors

To show the benefit of using a value of the regularization parameter different from ${\lambda }_{\mathbf{x}}^{*}$ when estimating the cross-power spectrum, in Figure 6 we plotted the reconstruction errors ${\epsilon }_{\mathbf{S}}\left(\lambda \right)$ as a function of the regularization parameter (normalized by ${\lambda }_{\mathbf{x}}^{*}$) obtained when considering two illustrative realizations of the simulated sensor data. Specifically, we fixed the locations and time courses of the pair of interacting sources and we considered the corresponding simulated MEG data for two levels of ${\mathrm{SNR}}^{\mathbf{X}}$, namely ${\mathrm{SNR}}^{\mathbf{X}}=-20$ dB and ${\mathrm{SNR}}^{\mathbf{X}}=5$ dB. Similar results where obtained when considered the other source configurations.

As shown by Figure 6, for both the values of ${\mathrm{SNR}}^{\mathbf{X}}$ the value of the reconstruction error significantly decreases when ${\lambda }_{\mathbf{S}}^{*}$ is used instead of ${\lambda }_{\mathbf{x}}^{*}$. Specifically in this simulation, ${\epsilon }_{\mathbf{S}}\left(\lambda \right)$ drops from $0.99$ to $0.96$ when ${\mathrm{SNR}}^{\mathbf{X}}=-20$ dB, and from $0.92$ to $0.77$ when ${\mathrm{SNR}}^{\mathbf{X}}=5$ dB.

Notably, one may observe that the relative reconstruction errors shown in Figure 6 are rather large, being above $90\%$ in the low-SNR case and remaining above $75\%$ even in the high-SNR scenario. We point out that this fact is mainly due to the combined effect of two factors: first, Tikhonov regularization tends to produce reconstructions that are small but non-zero almost everywhere, as it reduces but does not cancel entirely backprojection of noise; second, in our simulations the true activity is zero everywhere but in two points. These two facts inevitably lead to large relative errors that, however, pleasantly decrease for increasing values of ${\mathrm{SNR}}^{\mathbf{X}}$.

5. Discussion and Conclusions

In the present work, we investigated the role of the spectral complexity of a time series, $\mathbf{x}\left(t\right)$, in the design of an optimal inverse technique for estimating its cross-power spectrum, ${\mathbf{S}}^{\mathbf{x}}\left(f\right)$, from indirect measurements of the time series itself. Motivated by an analysis pipeline widely used for estimating brain functional connectivity from MEG data, we reconstructed the cross-power spectrum in two steps: first, we estimated the unknown time series by using the Tikhonov method, then, we computed the cross-power spectrum of the reconstructed time series. In the present work, we used numerical simulations to study how the spectral complexity of $\mathbf{x}\left(t\right)$ impacts the value of the regularization parameter that provides the best reconstruction of the cross-power spectrum.

As a first analytical result, we related ${\mathrm{SNR}}^{\mathbf{X}}$ to ${\mathrm{SNR}}^{\mathbf{S}}$, i.e., the signal-to-noise ratio of the time series and the signal-to-noise ratio of the corresponding cross-power spectra. The obtained formula suggests that, for a fixed level of ${\mathrm{SNR}}^{\mathbf{X}}$, ${\mathrm{SNR}}^{\mathbf{S}}$ depends on the spectral complexity of $\mathbf{x}\left(t\right)$: the higher the spectral complexity coefficient the higher ${\mathrm{SNR}}^{\mathbf{S}}$. Intuitively this happens because a higher value of the spectral complexity coefficient corresponds to a more peaked cross-power spectrum that will emerge over the cross-power spectrum of the noise.

To test the effect of this result on the choice of the Tikohonov regularization parameter in a practical scenario, we simulated a large set of MEG data and applied the described two-step approach for estimating the cross-power spectrum of the underlying neural sources. In details, we simulated 1200 synthetic MEG data with varying ${\mathrm{SNR}}^{\mathbf{X}}$ generated by pairs of coupled point-like sources at varying locations and with different spectral complexities. For each simulated data, we computed the two parameters providing the best estimates of the time series (${\lambda }_{\mathbf{x}}^{*}$) and of the cross-power spectrum (${\lambda }_{\mathbf{S}}^{*}$), defined as the ones minimizing the relative ${\ell }_2$ norm of the difference between the true and the reconstructed time series/cross-power spectrum according to Definition 2. As shown by Figure 4, the results of our simulations highlighted a high correlation between the values of ${\lambda }_{\mathbf{S}}^{*}$ and of ${\mathrm{SNR}}^{\mathbf{S}}$.

Eventually, we focused on the relationship between the two parameters ${\lambda }_{\mathbf{x}}^{*}$ and ${\lambda }_{\mathbf{S}}^{*}$, whose ratio is shown in Figure 5. The figure points out that this ratio depends on the spectral complexity of the simulated time series. This fact may be understood in lights of the previous results, as ${\lambda }_{\mathbf{S}}^{*}$ depends on ${\mathrm{SNR}}^{\mathbf{S}}$ that in turns depends on the spectral complexity coefficient. Additionally, we found that, for all the simulated data, $\frac{{\lambda }_{\mathbf{S}}^{*}}{{\lambda }_{\mathbf{x}}^{*}}<\frac{1}{2}$, in line with the results shown in [27] for a simplified model where the neural time series were assumed to be white Gaussian processes. Moreover, when the spectral complexity coefficient increases ($c>5$ in our simulations) the ratio between the two parameters approaches 0.01. This agrees with the results shown in [26] where, by simulating sinusoidal signals, the authors suggested to use for connectivity estimation a parameter of two orders of magnitude lower. In fact, our numerical results indicate that the use of ${\lambda }_{\mathbf{S}}^{*}$ results in a substantially lower reconstruction error on the cross-power spectrum, particularly when the data have a high SNR.

The present work focuses on the cross-power spectrum as a connectivity metric. Even though the cross-power spectrum is the starting point for the computation of many connectivity metrics it would be interesting to directly investigate the behavior of the Thikonov regularization parameters when using such metrics. Future works will be devoted to this. It is also worth noticing that the definition of optimality when defining the regularization parameters is not univocal, since many metrics can be used. A common example is the area under the curve (AUC), which is the metric that was used in [26]. The use of different metrics would firstly strengthen our results and would also allow a more straightforward comparison with the results of [26]. Finally, the dependence of ${\lambda }_{\mathbf{S}}^{*}$ on ${\mathrm{SNR}}^{\mathbf{S}}$ suggests that an analysis of such dependence could be considered for the definition of a rule for choosing ${\lambda }_{\mathbf{S}}^{*}$ in practical scenarios.