The role of spectral complexity in connectivity estimation

The study of functional connectivity from magnetoecenphalographic (MEG) data consists in 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.


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 to 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 Figure 1. Schematic representation of the differences between the regularization parameter providing the best time series estimate (λ x ) and the one providing the best cross-power spectrum estimate (λ S ). The first one provides an optimal reconstruction of the neural activity, but it may not lead to an optimal estimate of the cross-power spectrum; vice versa λ S provides an optimal reconstruction of the cross-power spectrum at the expense of a sub-optimal estimate of the time series. 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 optimize 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 a paradigmatic inversion technique and we study the interplay between the regularization parameter providing the reconstructed neural time series minimizing the relative error in 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 [23] 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. [24] 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 modelled as a multivariate autoregressive process.
The plan of the paper is as follows. Section 2 introduces the problem in a formal way. Section 3 describes how the synthetic data are simulated and analysed. Section 4 presents the results of the analysis. Our conclusions are offered in Section 5.

Forward model
Let X(t) = (X 1 (t), . . . , X N (t)) ∈ R N be a multivariate stationary stochastic process whose realisations x(t) can not be observed and let Y(t) = (Y 1 (t), . . . , Y M (t)) ∈ R M be the process whose realizations y(t) are used to infer information on x(t). Let Y(t) and X(t) be related by the following equation where G ∈ R M×N is the forward matrix and N(t) = (N 1 (t), . . . , N M (t)) ∈ R M is the measurement noise, that is here assumed to be a white Gaussian process with zero mean and covariance matrix α 2 I, i.e. N(t) ∼ N (0, α 2 I), independent from X(t).

Cross-power spectrum
We are interested in reconstructing the cross-power spectrum of X(t), which describes the statistical dependencies between each pair of time series (X j (t), X k (t)) j,k∈{1,...,N} . The cross-power spectrum is a one parameter family of N × N matrices S X ( f ), whose (j, k) − th element is defined as whereX j ( f , T) is the Fourier transform of X j (t) over the interval [0, T], defined aŝ and X H is the Hermitian transpose of X [25] Given a realization x(t) of the process X(t), the cross-power spectrum S X ( f ) can be estimated via the Welch method [26], which consists in partitioning the data in P overlapping segments multiplied by a window function, {w(t)x p (t)} P p=1 , computing their discrete Fourier transformx p ( f ) = 1 L ∑ L−1 t=0 x p (t)w(t)e −2πit f L and averaging: where L is the length of each segment and W = 1 L ∑ L−1 t=0 w(t) 2 . It is often the case that the data reaches 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 x(t) of the process X(t), and the corresponding cross-power spectrum S x ( f ), 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 2 −norm over the frequencies of S x j,k ( f ), j, k = 1, . . . , N, that is 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 S x ( f ) are considered because S x ( f ) is Hermitian.

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 x(t), an estimate of the cross-power spectrum can be obtained from y(t). Such estimate can be achieved through a two-step process [22]: is obtained by solving the inverse problem associated to equation (1).
Here we consider the Tikhonov regularized solution [27] of the problem which is defined as where λ is a proper regularization parameter and · 2 is the 2 -norm. ii. Then, the corresponding estimate of the cross-power spectrum S x λ ( f ) is computed from the reconstructed time series using the Welch method, as described in the previous section.
When applying this two-step process the regularization parameter λ in equation (6) has to be set for the computation of x λ (t). Thus, the problem naturally arises of the choice of such parameter, which can be set in order to optimally reconstruct either x λ (t) or S x λ ( f ). We define optimality through the minimization of the normalised 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 x(t) as and the optimal parameter for the reconstruction of S x ( f ) as where · F is the Frobenius norm; ε x (λ) and ε S (λ) 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 λ. In fact, they assume their maximum value when either λ is very high and thus x λ (t) is negligible with respect to x(t), or when λ is too small and thus, vice versa, x(t) is negligible with respect to x λ (t). This definition may appear overly complex compared to, e.g., a mere 2 -norm of the difference; however, in the presence of sparse data where only few time series are non-zero, the simple 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 2 -norm. A similar definition has been introduced in [28].
In experimental contexts, where x(t) is not known, the choice of the optimal regularization parameter is crucial. This matter is widely discussed in literature [29][30][31][32], and many criteria have been proposed. Such criteria apply to equation (1) and can be used to set the regularization parameter λ 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 X(t) related to such model as 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 will be presented in Section 4; however we first need to relate the cross-power spectrum of the unknown S X ( f ) with that of the data S Y ( f ).
By computing the cross-power spectrum of both sides of equation (1) and from the linearity of the Fourier transform it follows that where the mixed terms S XN ( f ) and S NX ( f ) are negligible thanks to the independence between X(t) and N(t). 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 S X ( f ) related to such model as This definition is in line with the definition of SNR X for the signal, the main difference being in the use of the Frobenius norm rather than the 2 -norm, motivated by the fact that we are working with matrices rather than vectors.

Generation and analysis pipeline of the MEG simulated data.
In this section we will 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.

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,33,34]. The two equations that are of interest here read as where E(r, t) and B(r, t) are the electric and magnetic fields at location r and time t, µ 0 is the magnetic permeability in vacuum and J(r, t) is the total electric current that flows inside the brain. The latter is the sum of two contributions J(r, t) = J p (r, t) J p (r, t) being the primary current directly related to the brain activity, while J v (r, t) = −σ(r)∇V(r, t) is the induced volume current due the non-null conductivity σ(r) of the brain, V(r, t) being the electric scalar potential. The manipulation of Maxwell's Equations leads to the Biot-Savart Equation where Ω is the volume occupied by the brain, the first term B 0 (r, t) = µ 0 4π Ω J(r , t) r−r |r−r | 3 dv is the magnetic field induced by the primary current while the second term is related to the volume current.
Solving the forward problem requires the computation of these two contributions knowing the primary current. While 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 {Ω j } j=1,...,J and to replace volume integration with surface integration. In this way Biot-Savart equation becomes where ∂Ω i,j is the contact surface between regions Ω i and Ω j and n i,j (r ) is the unit vector normal to the surface ∂Ω i,j at r 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 (16) after having computed V(r, t) by solving the equation which follows from equation (13) by applying the divergence. For further details on the MEG forward problem we refer the reader to [1].

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 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 y(t) = (y 1 (t), . . . , y M (t)). Now, observing that the magnetic field B depends linearly on the primary current J p , the magnetic field in correspondence of the sensors of the instrument is where r k , k = 1, . . . , N, is the location of the k-th brain parcel, G(r k ) ∈ R M×3 is the corresponding leadfield matrix, and {q k (t)} k=1,...,N are the electric current intensities along the three orthogonal direction of the N dipoles within the brain at time t and n(t) is the measurement noise. The l-th column of G(r k ) contains the measurement at a sensor level when a unit current dipole is placed at location 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 [35]. In this case the electric current intensities are scalar quantities (we refer to them as {q k } k=1,...,N ) and the leadfield matrices are column vectors (we refer to them as {G k } k=1,...,N ).
Let us define x(t) := (q 1 (t), . . . , q N (t)) (19) and reassembling equations (19) and (20) in to equation (18), we get which can be interpreted as a realization of equation (1). From now on we will refer to 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 [36]. 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.

Data generation
We simulated N mod = 10 pairs of active sources, (z 1 (t), z 2 (t)) , with unidirectional coupling from the first to the second; their time series follow a multivariate autoregressive (MVAR) model of order P = 5 [37,38] The non-zero elements a i,j (k) of the coefficient matrices were drawn from a normal distribution of zero mean and standard deviation γ, and T=10000. We retained only coefficient matrices providing (i) a stable MVAR model [37] and (ii) pairs of signals (z 1 (t), z 2 (t)) such that the 2 -norm of the strongest one was less than 3 times the 2 -norm of the weakest one. In order to obtain time series with different spectral complexity coefficients we set γ 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. Finally, the resulting time series (z 1 (t), z 2 (t)) 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 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 SNR X evenly spaced in the interval [−20dB, 5dB].
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.

Inverse model
Source space time series were reconstructed using the Tikhonov method, also known as minimum norm estimate (MNE) [2] within the MEG community. For each combination of source time series, source locations and SNR X level we computed the optimal regularization parameters λ * x and λ * S by minimizing the reconstruction errors ε x (λ) and ε S (λ), defined in Definition 2. The minimization procedure was achieved by using the Matlab built in function fminsearch. 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.

Results
In this section we illustrate the results of our analysis. We will begin with the description of the analytical dependence between SNR X and SNR S , then we will highlight how the optimal parameter for the reconstruction of the cross-power spectrum depends on SNR S and how this implies that the spectral complexity of the signal is behind such dependence. As a byproduct, this analysis also confirms the results of Vallarino et al. [24] in the case of a more complex setting.

MVAR time-courses simulation
Assign time-courses to sources in the brain

Projection to sensor level
Add sensor noise MEG data

Neural activity estimation
Cross-power spectrum estimation Inverse problem: λ has to be chosen Step 1 Compute cross-power spectrum: keep the same λ Step 2 Forward problem Inverse problem and SNR S = 10 log 10 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 (24) we used the fact that the cross-power spectrum of a white noise Gaussian process of zero mean and covariance matrix α 2 I is S N ( f ) = α 2 I. By isolating α 2 from equation (23) and substituting in equation (24) we obtain SNR S = 10 log 10 Equation (25) relates the the signal-to-noise ratio of X(t) with that of S X ( f ). It shows that, for same levels of SNR X , SNR S changes with the spectral complexity coefficient of the signals. In fact, the higher the spectral complexity coefficient, the higher the quantity GS X ( f )G 2 F . 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.

Dependence of λ *
S on SNR S As described in Section 3 we simulated several sensor level configurations, based on different combinations of spectral complexity coefficients, source locations and SNR X levels. For each configuration we collected the two optimal parameters λ * x and λ * S and we investigated their dependence on the signal-to-noise-ratio. In accordance with classical results from inverse theory [30], we found that λ * x depends on the signal-to-noise ratio. What is novel here is the relation between λ * S and both SNR X and SNR S . Indeed for increasing SNR X less regularization is needed, but such dependence varies with the MVAR models. On the other side, the dependence of λ * S on SNR 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 SNR X are shown, while on the right the same parameters are shown with respect to SNR 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.

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

Discussion and conclusions
In the present work, we investigated the role of the spectral complexity of a time series, x(t), in the design of an optimal inverse technique for estimating its cross-power spectrum, S x ( f ), 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 Figure 4. Optimal regularization parameters for the reconstruction of the cross-power spectrum (λ * S ) as a function of SNR X (left) and SNR S (right). Different colors correspond to different MVAR models. On the left, the lines have different heights, while on the right they overlap, meaning that the dependence of λ * S on SNR S is neater with respect to SNR X . present work, we used numerical simulations to study how the spectral complexity of x(t) impacts the value of the regularization parameter that provides the best reconstruction of the cross-power spectrum.
As a first analytical result, we related SNR X to SNR 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 SNR X , SNR S depends on the spectral complexity of x(t): the higher the spectral complexity coefficient the higher SNR 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 we 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 SNR 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 estimate of the time series (λ * x ) and the best estimate of the cross-power spectrum (λ * S ). As shown by Figure 4, the results of our simulations highlighted a high correlation between the values of λ * S and of SNR S . Eventually, we focused on the relationship between the two parameters λ * x and λ * 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 λ * S depends on SNR S that in turns depends on the spectral complexity coefficient. Additionally, we found that, for all the simulated data, λ * S λ * x < 1 2 , in line with the results shown in [24] 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 [23] where, by simulating sinusoidal signals, the authors suggested to use for connectivity estimation a parameter of two orders of magnitude lower.
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 behaviour 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 [23]. The use of different metrics would firstly strengthen our results and would also allow a more straightforward comparison with the results of [23]. Finally, the dependence of λ * S on SNR S suggests that an analysis of such dependence could be considered for the definition of a rule for choosing λ * S in practical scenarios.
Funding: EV, AS, SS and MP have been partially supported by Gruppo Nazionale per il Calcolo Scientifico.

Conflicts of Interest:
The authors declare no conflict of interest.