Lag-invariant detection of interactions in spatially-extended systems using linear inverse modeling.

Measurements on physical systems result from the systems' activity being converted into sensor measurements by a forward model. In a number of cases, inversion of the forward model is extremely sensitive to perturbations such as sensor noise or numerical errors in the forward model. Regularization is then required, which introduces bias in the reconstruction of the systems' activity. One domain in which this is particularly problematic is the reconstruction of interactions in spatially-extended complex systems such as the human brain. Brain interactions can be reconstructed from non-invasive measurements such as electroencephalography (EEG) or magnetoencephalography (MEG), whose forward models are linear and instantaneous, but have large null-spaces and high condition numbers. This leads to incomplete unmixing of the forward models and hence to spurious interactions. This motivated the development of interaction measures that are exclusively sensitive to lagged, i.e. delayed interactions. The drawback of such measures is that they only detect interactions that have sufficiently large lags and this introduces bias in reconstructed brain networks. We introduce three estimators for linear interactions in spatially-extended systems that are uniformly sensitive to all lags. We derive some basic properties of and relationships between the estimators and evaluate their performance using numerical simulations from a simple benchmark model.

Introduction

Global patterns in complex systems emerge through interactions between large numbers of units. One way of characterizing the behavior of such systems is by applying network theory to the matrix of pair-wise interactions [1]. If inversion of the forward model is ill-posed, this approach leads to bias in the estimated interaction matrix and hence to distortions of reconstructed network topology [2]. Bias in estimated interaction matrices is well-known in the field of MEG and EEG and prevents a clear view on the functional organization of the human brain. Another example are local field potentials (LFP’s), which are used in invasive neurophysiological studies [3]. Although less ill-posed than the EEG/MEG inverse problems, the LFP inverse problem still requires regularization and this leads to bias in the reconstruction of current generators in neural tissue [4]. Although linear, these operators have large null-spaces and are very sensitive to small perturbations, as evidenced by the high condition numbers of their discretizations. In the current study, we propose a method to reduce the bias in estimated interaction matrices that can be applied to observations obtained from arbitrary linear forward models.

Estimation of system interactions from indirect sensor measurements is particularly well-developed in the field of EEG and MEG, which are non-invasive techniques for detecting magnetic fields outside the head (MEG) and electric potentials on the scalp (EEG) set-up by current generators in activated brain tissue [5]. In contrast to blood-oxygen level-dependent (BOLD) functional magnetic resonance imaging (fMRI), MEG and EEG provide direct measures of neural activation with high temporal resolution and, as such, are indispensable tools in fundamental and clinical human neuroscience. Although direct analysis of MEG and EEG sensor data can be useful in disentangling cognitive processes, seizure detection, and classification of pathologies, it does not provide insight into the spatiotemporal patterns of neural activation underlying the sensor data. This applies in particular to functional interaction analysis [6, 7], which aims to characterize the temporal relationship between active sources [8], and is therefore commonly carried out in source-space [6, 9]. Source-space interaction analysis, however, generally yields biased estimators due to residual mixing of the reconstructed source signals. This undesirable effect is referred to as field spread in the case of MEG and volume-conduction in the case of EEG and is jointly referred to as signal leakage. Signal leakage is the main obstacle in the non-invasive study of human brain interactions and has attracted considerable attention from the research community.

Methods of dealing with signal leakage in source-space interaction analysis can roughly be divided into two categories, corresponding to the type of interaction that is considered: signal or amplitude interactions. Signal interactions refer to interactions between reconstructed source-signals proper, whereas amplitude interactions refer to interactions between the signals’ amplitude fluctuations. Signal interactions characterize interactions between fast and usually oscillatory signals on short time-scales, typically a fraction of the signals’ oscillation period. Examples are the Pearson correlation coefficient, coherence, and phase-synchronization [10]. Amplitude interactions characterize interactions between slowly varying amplitude envelopes on time-scales in the BOLD-fMRI range (0.01-0.1 Hz). Methods of reducing leakage in signal interactions are based on the observation that leakage is instantaneous and thus cannot explain the observation of interactions with non-zero lag. Methods of this type therefore quantify lagged interactions (i.e. interactions with non-zero lag) and discard all zero-lag interactions [11-15]. Methods of reducing leakage in amplitude interactions remove all linear instantaneous interactions between the reconstructed source signals, before computing their amplitude envelopes. This is done by applying a whitening transformation to the signals and is referred to as signal orthogonalization. The amplitude envelopes are subsequently low-pass filtered and functional interactions are quantified by the Pearson correlation coefficient [16-19]. This approach has enabled the discovery of the electrophysiological basis of human hemodynamic resting-state networks [20, 21].

This study is concerned with signal interactions. Although the above mentioned methods have proven valuable in various applications, their sensitivity decreases with decreasing (absolute) lag, which leads to false negatives and distortions in the topology of reconstructed functional networks [2]. This is a serious limitation, especially given that intra-cortical electrophysiological signals are known to interact with small lags under many experimental conditions and have been hypothesized to be critical for effective neuronal communication and routing of information [22, 23]. The main reason why measures of signal interactions are unable to separate true from spurious zero-lag interactions is that the formal relationship between signal leakage, and the forward and inverse operators is not exploited. In [24, 25] this relationship is exploited for linear inverse operators by linking signal leakage to the algebraic structure of the resolution operator. This method, however, is designed for seed-based interactions and is asymmetric, without having a causal interpretation, which complicates its interpretation. Furthermore, straightforward symmetrization and generalization to the multivariate case leads to a variation on an earlier proposed bias-correction scheme that has been shown to yield only marginal improvements [26]. In [27], a subspace projection method is proposed for suppressing the contribution of source-power to the sensor-space spectral matrix. This yields a bias-corrected spectral matrix, which is subsequently projected to source-space using a non-linear scanning method. The method is based on the observation that contributions to the sensor-space spectral matrix that are due to source-power (spurious interactions) and source coherence (true interactions) lie in different linear subspaces of the vectorspace of n × n complex-valued matrices, where n denotes the number of sensors. It is in this last step that the relationship between signal leakage and the forward matrix is exploited. The source-space spectral matrix is subsequently estimated by carrying out a non-linear scanning procedure in which the corrected sensor-space spectral matrix is compared with the (normalized) cross-product of leadfields of every pair of source locations. The main limitation of this method is the use of a non-linear scanning method, as such methods are known to be effective only for a small number of active sources [28].

The bias-reduction method proposed in [27] can in principle be combined with any inverse operator. This was noted in [27] but not worked out. In the current study we combine the bias-correction proposed in [27] with an arbitrary linear inverse operator to obtain a lag-independent estimator of system interactions that can be applied to reconstruct interactions in distributed source-activity. This requires to invert a linear forward model defined on vectorized matrix spaces. We also propose two novel bias-corrected estimators. The first estimator is obtained by interchanging the order of inversion and bias-correction. Thus, instead of correcting the sensor-space spectral matrix and subsequently projecting the corrected matrix to source-space, the sensor-data is first projected to source-space to yield reconstructed system activity, from which the spectral matrix can be estimated and corrected using a subspace projection in source-space. The second estimator is obtained by first solving the same vectorized inverse problem as above and subsequently correcting the reconstructed spectral matrix using the same source-space projection as above. This last estimator hence combines the ingredients of the other two estimators.

In Background we describe the linear forward model and linear inverse operator used in this study (), define system interactions in the time-frequency domain using spectral matrices and formulate null-hypotheses and the construction of confidence regions (), introduce the benchmark model that will be used in evaluating the proposed methodology (), and motivate the use of complex-valued test-statistics for lag-independent detection of system interactions (). In Estimation of system interactions we describe uncorrected source- and sensor-space based estimators of system interactions and derive their filtering properties. In The geometry of signal leakage we provide a geometric characterization signal leakage in linear estimators of system interactions and provide some analytical examples using simple source configurations. In Suppression of signal leakage we use this characterization to construct three corrected interaction estimators and describe how their efficacy can be quantified. In Comparative performance we compare the performance of the corrected estimators with respect to suppression efficacy (), bias-reduction (), and detection power ().

Background

Linear inverse modeling

In the time-frequency domain, the electromagnetic forward model can be represented by a linear where x(t, ω) is a p-dimensional vector of brain activity, y(t, ω) is an n-dimensional vector of sensor measurements, L is the forward matrix, e(t, ω) is an n-dimensional vector of measurement noise, and t and ω denote time and angular frequency, respectively. A model of this form describes the relation between brain activity in the form of current source density (CSD) at p locations within the brain and induced measurements at n sensors, which can be electric potentials on the scalp as in electroencephalography (EEG), on the cortical surface as in electrocorticography (ECoG), or inside the brain as in local field potentials (LFP) recordings, as well as magnetic flux outside the head as in magnetoencephalography (MEG). These different types of measurement correspond to different forward matrices, which contain the relevant geometric and conduction properties of the volume conductor under study (i.e. the head) and are obtained by discretizing and solving the quasi-static approximation to Maxwells’ equations [5]. Electromagnetic inverse modeling refers to the reconstruction of the CSD or any derived quantities from observed sensor measurements. What makes inverse modeling challenging is the fact that p is generally much larger than n.

A popular way of obtaining a reconstruction x⋆(t, ω) of x(t, ω) is by solving the following penalized least-squares problem: where λ ≥ 0 is the noise regularization parameter [5, 28]. This problem has the unique solution where the inverse operator is given by

In the neuroimaging community, x⋆(t, ω) is referred to as the minimum norm reconstruction of x(t, ω) and several generalizations of it have been developed [28]. The properties of a linear inverse operator are summarized by its resolution matrix. Given a forward matrix L and an inverse operator L♯, the resolution operator associated with L and is defined as

Its columns are known as point-spread functions and its rows is cross-talk functions [29, 30]. Note that in case of a single active source of unit strength at the j-th location, the reconstruction is equal to the j-th point-spread function. We will need these notions later on when analyzing the functional connectivity estimators.

Detection of system interactions

We focus on the detection of interactions from ongoing recordings of the system’s activity i.e. without perturbing the system. In neuroscience, this corresponds to recording brain activity in the absence of sensory stimulation and cognitive tasks. The proposed methodology, however, can easily be adapted to detect interactions from recorded system responses to perturbations. We assume that ongoing system activity can be reasonably well modeled by a Gaussian stochastic process. Without loss of generality, we assume the process to have expectation zero. Such a process is completely determined by its spectral matrix S(ω), which is a function of (angular) frequency ω and is defined as where the brackets denote taking the expectation over time and the superscript † denotes taking the conjugate-transpose. The spectral matrix of a process is positive semi-define and conjugate symmetric, i.e. , where the bar denotes entry-wise complex-conjugation and can hence be thought of as a covariance matrix in the time-frequency domain. Measurement noise is also modeled by a Gaussian stochastic process and for simplicity we will assume that S(ω) = σ2 I, where I denotes the identity matrix i.e. the measurement noise is white, has the variance at each sensor, and is uncorrelated across sensors. The observed sensor data y(t, ω) therefore is also a stochastic process and has a spectral matrix S(ω) given by

In neuroscience, system interactions are referred to as functional connectivity, which is a broad term that encompasses notions based on different frameworks such as stochastic processes, weakly-coupled oscillators, chaos theory, and information theory [10]. We will adopt the framework of stochastic processes and hence characterize system interactions by the (non-diagonal entries of) the spectral matrix S(ω). Thus, the interactions between system locations k and l at frequency ω is characterized by the (k, l)-th entry of S(ω), which we will denote by γ(ω) and is a complex number whose magnitude and phase measure, respectively, the strength and lag (i.e. latency) of the interaction at frequency ω.

This study is concerned with the construction of test-statistics for detecting the presence of system interactions based on observed sensor data. This problem is formalized by the following null-hypothesis: against the two-sided alternative H1: γ(ω)≠0, for all pairs (k, l). The test-statistics we consider are complex-valued estimators of γ, together with their real and imaginary parts and are described in detail in and . The performance of the test-statistics will be evaluated by computing the distance between H0 and 100(1 − α)% confidence regions. Throughout the text we will refer to this distance as the sensitivity of a test-statistic and to the matrix carrying the sensitivities of all pairs (k, l) as the sensitivity matrix. For complex test-statistics, 100(1 − α)% ellipsoidal confidence regions in the complex-plane are constructed and for its real and imaginary parts, 100(1 − α)% confidence intervals are constructed by simulating the statistics’ sample distribution and calculating the 100(1 − α/2)% and 100 α/2% percentiles. To compare the sensitivity of a complex test-statistic to that attained by the combined use of its real and imaginary parts, we set α = 0.05 for the complex test-statistic and α = 0.025 for its real and imaginary parts. This is done to correct for the fact that in the latter case, two tests need to be performed instead of one. The sample distributions of the different test-statistics are approximated by generating 103 realizations of the simulated sensor data. In practice, these distributions can be approximated using appropriate bootstrap schemes.

Benchmark model

To illustrate the methodology, we use a simple forward model in which electric potentials are measured that are induced by a current source density (CSD). The electric potential V generated by a current source density ρ in an infinite, homogeneous, and isotropic volume conductor is described by Poisson’s equation: where Δ denotes the Laplace operator and where, without loss of generality, we have set the electrical conductivity to 1. In a neurophysiological context, V is the local field potential (LFP) inside nervous tissue that is generated by transmembrane currents which are described macroscopically by the current source density ρ [31, 32]. We consider a one-dimensional source-space segment of length 4 mm and a one-dimensional array of 11 electrodes with an inter-electrode spacing of 0.4 mm and that is located parallel to the source-space segment at a height of h = 0.5 mm (see Fig 1A).

Fig 1

Set-up of the test model.

A. Source-space segment of length 4 mm (horizontal axis at height 0 mm), recording electrodes (black dots at height 0.5 mm), and the locations of two active sources of neuronal activity (red dots). Also shown are the sensitivities of the electrodes to activity of the two sources (two black curves). B. Observed spectral matrix of the recorded electric potentials that are induced by non-interacting sources of unit strength and in the presence of measurement noise (σ = 0.1).

The electric potential at location y that is generated by a current monopole ρ(x′) = δ(x′ − x) of unit strength at location x is given by the free-space fundamental solution to Poisson’s equation:

The CSD is modeled by two monopoles located at locations -1 and 1 mm and correspond to the red dots in Fig 1A. The source-space is discretized using a spacing of 0.1 mm and hence gives rise to a 10 × 81 forward matrix L that describes the relation between the (unobserved) CSD and the (observed) electric potentials. The (k, l)-th entry of L is V(x, y), where x and y denote the locations of the k-th electrode and the l-th source-point, respectively. The forward model has the form of Eq (1).

The two active sources are modeled by stationary oscillatory stochastic processes with spectral matrix where 0 ≤ γ ≤ 1 is the interaction strength (i.e. coherence) and 0 ≤ ϕ < 2π is the lag. Thus, the sources have the same amplitude and interact with strength γ and a latency that is a fraction ϕ/2π of their average oscillation period. Measurement noise with spectral matrix S = σ2 I is added, where σ2 is measured in units of the maximal eigenvalue of LL. Table 1 lists the model parameters and their nominal values. Fig 1B shows the observed spectral matrix S that is generated by source activity with ϕ = 0, γ = 0, and σ = 0.1. Note that although the two active sources can be discerned, due to the mixing of source activity by the forward model (two black curves), it is not immediately clear if, and to what extent, the sources are interacting.

Table 1

Model parameters, their symbols, and nominal values/ranges.

Parameter	Unit	Value
Array length	mm	4.0
Segment length	mm	4.0
Intra-sensor distance	mm	0.4
Intra-source distance	mm	0.2
Height of array (h)	mm	0.5
Interaction strength (γ)	-	0.3
Interaction lag (ϕ)	deg	0-90
Noise level (σ)	-	0-0.1

Test-statistics

A commonly used statistic in neuroscience for testing the null hypothesis H0: γ = 0 of the absence of interactions between brain locations k and l is the imaginary part of some estimator of γ. Fig 2A shows Timag as a function of lag . It shows that Timag is highest for imaginary interactions i.e. when system activity at locations k and l is coherent with lag one-fourth of the average oscillation cycle, and decreases to zero for purely real (i.e. instantaneous) interactions. Similarly, the sensitivity of the real part is highest for instantaneous connectivity and decreases to zero for purely imaginary interactions (see Fig 2A). As a consequence, interactions with small/large lags might remain undetected. In neuroscience this leads to bias in measures that are derived from estimated interaction matrices such as network-theoretical measures that are often used to characterize the brain’s functional organization [1, 33, 34].

Fig 2

Real, imaginary, and complex test-statistics.

A. Real-and imaginary part of the imaginary test-statistic as a function of lag. The strength of interaction was set to γ = 1 so that the curves range between zero and one. B. Sensitivity of the real, imaginary, and complex test-statistics at the true interaction-pair (k0, l0), as a function of lag C. Sensitivity of the imaginary test-statistic for all interaction-pairs and for a lag of 0 degrees. D. Sensitivity of the real test-statistic for all interaction-pairs and for a lag of 0 degrees. E. Sensitivity of the imaginary test-statistic for all interaction-pairs and for a lag of 90 degrees. F. Sensitivity of the real test-statistic for all interaction-pairs and for a lag of 90 degrees.

To illustrate the dependence of the real and imaginary test-statistics on interaction lag, we computed their sensitivity as a function of lag, for being the (k0, l0)-th entry of the sample spectral matrix of the reconstructed source activity and k0 and l0 denote the true source locations. The number of samples was set to N = 50, the correlation level was set to γ = 0.3, and no measurement noise was added. Source activity was reconstructed using the inverse operator (see Eq (2)) where λ was set to a small value (λ = 10−20). All other parameters were chosen as in Table 1. The sensitivities are shown in Fig 2B. Observe that the sensitivity of Timag is zero for lags smaller than about 45 degrees, which means that interactions at lags less than 45 degrees remain undetected. For lags larger than 45 degrees, its sensitivity starts to increase and is maximal at a lag of 90 degrees.Treal essentially behaves in the opposite way: it is highest at lag zero and decreases towards zero with increasing lag. In particular, interactions with lags larger than about 80 degrees remain undetected. Important to note is that the non-zero values of Treal for lags between 45 and 80 degrees are in fact spurious and when these are corrected for, Treal essentially behaves the same as Timag but reflected in the 45 degrees axis.

The properties of Timag and Treal considered above are local in the sense that we only considered their behavior at a single pair of system locations. Although locally they behave the same, up to reflection, globally they behave very differently and this is in fact the reason why imaginary test-statistics have become popular in neuroscience [11-15]. We first consider the global behavior of the imaginary test-statistic. We consider the extreme cases of ϕ = 0 degrees and ϕ = 90 degrees. The sensitivities for all source-pairs are shown in Fig 2C (0 degrees) and Fig 2E (90 degrees). In the case of ϕ = 0 degrees, none of the source-pairs show significant interaction. Thus, on the one hand, the true interaction at (k0, l0) remains undetected, but on the other hand, no spurious interactions are present at any other source-pair. In the case of ϕ = 90 degrees, sensitivity is highest at the location (k0, l0) of true interaction and is otherwise low. Thus, also in this case, spurious interactions are (almost) absent. The global behavior of Timag therefore allows for the detection of interactions as long as the lag is not too small. Fig 2D shows that the global behavior of Treal is very different. In particular, for both ϕ = 0 and ϕ = 90 degrees, their are many source pairs that exhibit spurious interactions, especially around the true source locations. Moreover, although for ϕ = 0 degrees, the test-statistic is weakly sensitive to the true interaction, in practice it will not be detected because of the presence of spurious interactions. Thus, Treal does not allow for the detection of interactions, irrespective of the lag. Fig 2B also shows that, locally, the complex-valued test-statistic is roughly uniformly sensitive to interactions i.e. it is non-zero and roughly equal for all lags. Unfortunately, the global behavior of Tcompl does not allow for the detection of true interactions because of the high level of spurious connectivity in its real part.

Subspace projections for suppressing spurious interactions

The goal of this study is to reduce the spurious connectivity in Tcompl so that it allows for the detection of system interactions at all lags. Spurious connectivity will be reduced by correcting the test-statistics using appropriately defined subspace projections. We consider three types of corrections and the associated corrected source-space spectral matrices. The matrices differ in the order in which the inverse operator and the projection operators are applied and in which space the subspace projection is applied (source- or sensor-space). The first type is obtained by first projecting the time-frequency coefficients of the sensor data to source-space by using an arbitrary linear inverse operator. Subsequently, the reconstructed time-frequency coefficients of the source-activity are used to estimated the activity’s spectral matrix. This is done using the sample spectral matrix since this is an unbiased estimator. Lastly, the estimated spectral matrix is corrected by applying a suitably defined subspace projection. Note that this projection is defined in the vectorspace of source-space spectral matrices. Fig 3A provides an illustration. The second type is obtained by switching the order of the inverse and subspace projection operators as illustrated in Fig 3B. The second type is also obtained by switching the order, but applying the subspace projection in sensor-space instead of source-space as illustrated in Fig 3C. The third of the three types of corrections described above has been proposed in [27]. Note that whereas the source-based estimator makes use of the reconstructed source-activity, the two sensor-based estimators avoid this step by applying an inverse operator in the space of spectral matrices instead of in signal space. The difference between the two sensor-based estimators is in which space the subspace projection is applied (source- or sensor-space). We will provide a comparative performance analysis of these three types of corrected estimators.

Fig 3

Subspace projections for suppression of spurious interactions.

A. Source-based estimator with correction in source-space. First the time-frequency coefficients of the sensor data are projected to source-space using an arbitrary linear inverse operator. Next the spectral matrix of the reconstructed time-frequency coefficients are estimated and subsequently corrected by applying a subspace projection. B. Sensor-based estimator with correction in source-space. First the spectral matrix of the sensor-space time-frequency coefficients is estimated. Next, the estimated sensor-space spectral matrix is projected to source-space using the tensor product of an arbitrary linear inverse operator with itself and subsequently corrected by applying a subspace projection. C. Sensor-based estimator with correction in sensor-space. First the spectral matrix of the sensor-space time-frequency coefficients is estimated. Next, the estimated sensor-space spectral matrix is corrected by applying a subspace projection and subsequently projected to source-space using the tensor product of an arbitrary linear inverse operator with itself.

Estimation of system interactions

Source-space based estimation

Measurements in the time-frequency domain at a particular frequency ω yield a n × k complex-valued data matrix Y, where n and k denote the number of sensors and samples, respectively. The observed data is related to source activity X and sensor noise E through where X and E are complex-valued p × k and n × k matrices, respectively. If X could be observed directly, the best estimator for the spectral matrix S is the sample spectral matrix When only Y can be observed directly, a natural way to estimate S is to first reconstruct the source activity X be applying an inverse operator to the observed sensor data and subsequently to compute the sample spectral matrix of the reconstructed source activity. Let X⋆ = L♯ Y be the reconstructed source activity, where L♯ is some linear inverse operator, and define the estimator where denotes the sample spectral matrix of the observed sensor activity. Note that is Hermitian as required for a proper estimator of S. We will refer to this estimator as source-based to emphasize that it is defined in terms of the (reconstructed) source activity. It will be convenient to write it in vectorized form as where ⊗ denotes the (Kronecker) tensor product. The estimator can be expressed in terms of X and E by first writing where denotes the sample spectral matrix of the measurement noise and denotes the sample cross-spectral matrix between source activity and measurement noise. Consequently, where R1 = (L♯ L ⊗ L♯ L), and F1 = L♯ ⊗ L♯ are the associated resolution and inverse operators, respectively. Since source activity and measurement noise are uncorrelated and and are unbiased estimators of S and S, respectively, the expectation of is

We now compute the filtering properties of , which will be used in the next section to establish a relationship between the source- and sensor-based estimators. Let r ≤ n be the rank of L and let L = UDV be the singular value decomposition of L, where U is an n × r orthonormal matrix, V is an p × r orthonormal matrix, and D is the r × r diagonal matrix that holds the non-zero singular values δ1 ≥ δ2 ≥ ⋯ ≥ δ of L in decreasing order. Note that U and V are bases for sensor and source space, respectively. Using the singular basis of the forward matrix, the resolution operator R1 can be represented as where the filter coefficients are given by

The filter coefficients range between zero and one and measure to what extent a contribution of the form of the true source spectral matrix is retained in the reconstructed spectral matrix. In particular, the diagonal entries of the filter matrix (i.e. the gains) are

Sensor-space based estimation

The source-space based estimator is obtained by first projecting the observed sensor data to source-space and subsequently computing the sample spectral matrix of the reconstructed source-activity. An alternative way of estimating the source-space spectral matrix is to first estimate the sensor-space spectral matrix and subsequently to project it to source-space [27]. We will refer to this estimator as sensor-based and denote it by . It is obtained by writing the sample sensor-space spectral matrix in vectorized form and subsequently projecting it to source-space by using the inverse operator F2 = (L ⊗ L)♯: where

Alternatively, this estimator can also be derived by solving the following convex optimization problem: where the minimization is done over S. The gradient of this objective function has a unique minimum in .

The sensor-based estimator can be written in terms of X and E as where R2 = F2(L ⊗ L) is the associated resolution operator. Furthermore, its expectation it given by

Comparing the expressions for and , we see that both are obtained by applying a linear inverse operator to the vectorized sensor-data : For the source-based estimator, the inverse operator is L♯ ⊗ L♯ and for the sensor-based estimator, the inverse operator is (L ⊗ L)♯. It is easily verified that in the absence of regularization (i.e. λ = 0), both estimators reduce to the Moore-Penrose inverse of L⊗L. In the limit of strong regularization, the inverse operators corresponding and reduce to λ−1(L⊗L) and λ−2(L⊗L), respectively.

The general relation between and is obtained by comparing their filtering coefficients. The resolution operator R2 corresponding to can be represented as with filter coefficients

In the absence of regularization, the filter coefficients of both estimators are equal to one, so that all contributions of the form are completely retained. The gains are

Because of the extra term in the numerator of , the source-based estimator provides stronger damping than the sensor-based estimator. The expressions for the filter coefficients also show that there is a natural correspondence between λ1 and , which allows the estimators to be directly compared by setting .

The geometry of signal leakage

Signal leakage refers to bias in reconstructions of brain activity that is due to the incomplete unmixing of sensor signals by application of an inverse operator. When measuring electric potentials as in EEG, ECoG, and LFP, signal leakage is called volume conduction and when measuring magnetic fluxes such as in MEG, signal leakage is called field spread [6]. Here we provide a geometric characterization of signal leakage that pertains to the functional connectivity estimators defined in textbfEstimation of system interactions and that applies to any linear inverse operator L♯. This characterization will form the basis for defining functional connectivity estimators that actively suppress signal leakage.

The expectations of the estimators defined in are of the form for j = 1, 2. The second term on the right-hand-side of Eq (4) corresponds to signal leakage due to measurement noise that is projected to source-space by the inverse operator. The first term on the right-hand-side contains true source-space functional connectivity in vectorized form, but in a distorted way due to multiplication by R⊗R. For a square matrix M, let M+ and M− denote the matrices obtained from M by setting its off- and on-diagonal entries, respectively, to zero. Thus, M can be decomposed as M = M+ + M− and the same is true for its vectorization. We can hence write the expectations as

We now interpret each of the three terms on the right-hand-side of Eq (5).

Note that the first term on the right-hand-side of Eq (5) only depends on the diagonal entries of S, i.e. on the power of (active) sources and not on their interactions. In this sense, all non-zero off-diagonal entries correspond to spurious connectivity. We refer to this term as the leakage term. Eq (5) shows that the leakage term lies in the subspace of that is spanned by the 1 + (i − 1)(p + 1)-th columns of R for 1 ≤ i ≤ p. We write B for the matrix that contains these columns of R and refer to the column space of B as the leakage subspace. The leakage term can be written explicitly as where denotes the power of the i-th source and denotes the i-th point-spread function of R.

The second term on the right-hand-side of Eq (5) does not depend on source power but only on the interactions between (active) sources. Its off-diagonal entries are non-spurious in the sense that they vanish in the absence of interactions. Its entries are generally biased, however, except when r = e for 1 ≤ i ≤ p, which is impossible since n < p. We refer to this term as the interaction term. Eq (5) shows that it lies in the subspace of that is spanned by the columns of R that are not in B. We refer to this space as the connectivity subspace. The interaction term can be written explicitly as: where γ denotes the (i, i′)-th entry of S i.e. the interaction between source i and i′, and where the sum runs over all ordered pairs (i, i′) with i < i′.

Lastly, the third term on the right-hand-side of Eq (5) is independent of source activity and only depends on the measurement noise and as such is entirely spurious. We refer to this term as the (projected) noise term. Eq (5) shows that it is contained in the column space of F, which we will refer to as the (projected) noise subspace. It can be written explicitly as where μ denotes the (i, i′)-th entry of S and denotes the i-column of F, which can be thought of as a noise point-spread function of the i-th sensor. The sum runs over all ordered pairs (i, i′) with i ≤ i′.

To illustrate the effects of signal leakage on the estimated spectral matrix , we consider two simple source configurations in the absence of measurement noise. We first consider a pair of point sources at k and l, with variances and , respectively, that interact with strength γ. Thus, the leakage and interaction terms are given by and respectively. The squared norm of the estimated leakage term therefore is where the brackets 〈, 〉 denote the standard inner product, and the squared norm of the estimated connectivity term is given by

By comparing these expressions and using the fact that , which follows from the positive semi-definiteness of S, we see that the inequality is implied by which is trivially true. Therefore, for two interacting point sources, signal leakage is always stronger than (estimated) interaction strength.

As a second example we consider a spatially extended source covering m locations. If the source is not too large, then r ≈ r for all locations k, l within the source. Furthermore, if the source is homogeneous, say with variance σ2, then the leakage term is approximately equal to mσ2 r⊗r, where r denotes the shared point-spread function within the source (r ≈ r ≈ r). Moreover, if the interactions within the source are homogeneous, say with shared strength γ, then the connectivity term is approximately equal to γm(m − 1)r⊗r. Thus, the strength of signal leakage, relative to that of the estimated connectivity, is approximately equal to This shows that for spatially extended sources (i.e. m large), the estimated intra-source connectivity is dominated by interaction effects and signal leakage can hence be neglected. In particular, suppressing signal leakage is not necessary in this situation.

Suppression of signal leakage

Suppression in source-space

In we showed that can be decomposed into three terms that lie in different linear subspaces of and correspond to leakage, interaction, and projected noise, respectively. We now define a sequence of linear operators acting on that increasingly suppress the leakage term, while retaining the interaction term as much as possible. To construct this sequence, we exploit the fact that the leakage subspace is spanned by the columns of the matrix B. Thus, if is a singular value decomposition of B and the dimension of the leakage subspace is d, the first d columns of U form an orthonormal basis for the leakage subspace. For 0 ≤ k ≤ d, we define the linear operator on by where is the matrix carrying the first k columns of U. Note that is the projection onto the orthogonal complement of the span of the first k left-singular basis vectors of B. In the case k = 0, the projection reduces to the identity operator and in the case k = d, the projection nullifies any vector that lies in the leakage subspace.

We now define a sequence of corrected interaction estimators by for 0 ≤ k ≤ d. Note that and that larger projection ranks k provide increasingly stronger suppression of the leakage term. In particular, for k = d, the leakage term is completely suppressed. The expected value of the corrected estimator is given by where is the corrected resolution operator and is the corrected inverse operator.

An interesting property of the operator is that it does not alter the imaginary part of the interaction term . This follows from the fact that the subspaces of symmetric and anti-symmetric matrices are orthogonal to each other. Indeed, the column space of is contained in the space of symmetric matrices, so the space onto which projects contains the space of anti-symmetric matrices. Since the imaginary part of the interaction term is anti-symmetric, it is invariant under the action of . Thus, does not suppress the imaginary, i.e. lagged part of the system interactions. In the rest of the text, we will refer to this type of suppression as suppression in source-space. Also note that a global scaling of the interaction-strength scales the projected interaction term by the same factor and leaves the projected leakage term unchanged. Specifically, let μ ≥ 0 be a scaling factor and consider the scaled source-space spectral matrix , then

Suppression in sensor-space

As proposed in [27], a similar correction can be applied in sensor-space. Instead of first projecting the sensor-data to source-space and subsequently computing and correcting the sample source-space spectral matrix, first the sample spectral matrix of the observed sensor-data is computed and corrected, and subsequently projected to source-space. From Eq (3) it follows that the expected value of the vectorized sample sensor-space spectral matrix is given by which can be decomposed as similar to the decomposition of in Eq (5). This allows to define sensor-level equivalents of the leakage, interaction, and projected noise subspaces defined in as well as a sensor-level equivalent of the leakage correction . Concretely, let B be the matrix carrying the 1 + (i − 1)(p + 1)-th columns of L⊗L for 1 ≤ i ≤ p, let B = UΛV be a singular value decomposition of B, and let d3 denote the dimension of the (sensor-level) leakage subspace. Then the first d3 columns of U form an orthonormal basis for the (sensor-level) leakage subspace and we hence can define the projection on by where U( is the matrix carrying the first k columns of U. Completely analogous to the definition in , we define a sequence of corrected sensor-space spectral estimators by for 1 ≤ k ≤ d3, and subsequently project these to source-space to obtain the following corrected estimator for the source-space spectral matrix: where the inverse operator is given by

Note that so that and correspond to different ways to correct . The expected value of is given by where is the associated resolution operator. Thus, the expected value of the sensor-space corrected estimator has the same form as the source-space corrected estimators and therefore has the same leakage structure in source-space. We do observe that, in contrast to the source-space corrected estimators, the norm of does not need to be a decreasing function of k. In particular, although for k = d3 the sensor-space leakage term is completely suppressed, this need not be true for the source-space leakage term . Like the sensor-space corrected estimators, the imaginary part of the system interactions is not suppressed. In the rest of the text, we will refer to this type of suppression as suppression in sensor-space.

Effectiveness of suppression

The effectiveness of the projections defined in the previous two sections in suppressing signal leakage in source-space interactions does not only depend on the extent to which the leakage term is suppressed, but also on the extent to which the interaction term is suppressed. Indeed, because the leakage, interaction, and noise terms are not separately observable, the projections need to be applied to their sum, i.e. the estimated spectral matrix. Consequently, they might not only suppress leakage, but also partly suppress interactions and projected sensor noise. The strengths of the corrected leakage and interaction terms, relative to those of the respective uncorrected terms can be quantified using the generalized Rayleigh quotient. Let S be any source-space spectral matrix. We consider the generalized Rayleigh quotient with positive semi-definite matrices and G where is the Gram matrix of the corrected resolution operator . Thus, for any source-space spectral matrix S, the generalized Rayleigh quotient is given by

One minus the Rayleigh quotient quantifies the extent to which the expected value of the uncorrected estimator is reduced by the correction . Below we will refer to this as the suppression level. The suppression levels of the leakage and interaction terms are defined as and , respectively. The fact that does not suppress the imaginary part of the system interactions can be expressed as for all k. For the source-space corrections (i.e. j = 1, 2) the range of is confined to the unit interval because then reduces to . Furthermore, complete suppression of S by a rank-k correction is equivalent to being contained in the column space of and no suppression is equivalent to being contained in the orthogonal complement of . In the ideal case, the leakage term is completely suppressed and the interaction term is entirely retained, i.e. and for some k. This can also be interpreted geometrically by noting that the angle α(S) between R Vec(S) and the column space of U equals which shows that the suppression level and α(S) are related through a monotonically decreasing function.

Comparative performance

Effectiveness

We calculated the suppression levels of the interaction term as a function of the regularization-level λ for five increasing interaction lags of 0, 20, 40, 50, and 70 degrees. Suppression of interactions with a lag of 90 degrees is zero because the imaginary part of interactions is invariant under the action of the correction projection. All other parameter values were set as in Table 1. The projection ranks were set to their maximal values so that the leakage terms are completely suppressed. Fig 4A shows the results for the source-based estimator. Observe the presence of three regimes: An under-regularization regime (λ < −2) in which the interaction is not suppressed, a regime in which the level of regularization is appropriate (−2 < λ < 1) and in which the level of suppression sensitively depends on the value of λ, and an over-regularization regime (λ > 1) in which the interaction is completely suppressed. This behavior can be observed for all lags and for the sensor-based estimators as well (see Fig 4B and 4C). The absence of suppression in the under-regularization regime implies that the interaction term is contained in the orthogonal complement of the leakage subspace. As λ increases, the angle between the interaction term and the leakage subspace decreases and this is reflected in stronger suppression. For large values of λ, the angle approaches zero so that the interaction term will be contained in the leakage subspace and hence will be completely suppressed. These results make clear that the effectiveness of the correction decreases with increased regularization-levels. In practice this means that the correction is less effectiveness for noisy data. It is therefore not advised to use regularization-levels that have been optimized for the estimation of source power instead of source interactions. For [35] have demonstrated that the optimal regularization-level for interactions is lower than that for source-power, at least for the ridge inverse operator (see also [36]).

Fig 4

Suppression levels.

A. Suppression levels of the source-based estimator as a function of regularization-level and for three interaction lags (0, 20, 40, 50, and 70 degrees). B. Same as in A. but for the sensor-based estimator with suppression in source-space. C. Same as in A. but for the sensor-based estimator with suppression in sensor-space. D. Suppression of projected sensor noise for all three estimators and as a function of regularization-level. In panels A, B, and C, the interaction lags increase in the direction of the arrows. In all panels, the parameters were set as in Table 1.

Fig 4 also shows that for all three estimators and for all regularization-levels, suppression levels decrease with increasing interaction lag ϕ. This is a consequence of the imaginary part of the interactions being unaffected by the corrections and can be understood in the following way. Let s = γcosϕ + iγsinϕ be the interaction with strength γ and lag ϕ. The effect of the correction on s is that its real part is suppressed, and its imaginary part is retained. Hence, applying the correction to s gives rγcosϕ + iγsinϕ, for certain 0 ≤ r ≤ 1, where r = 1 and r = 0 correspond to no and complete suppression, respectively. The suppression level is now given by which is a decreasing function of the interaction lag ϕ. Since the suppression levels are (roughly) increasing functions of the regularization-level, not exceeding a given upper-limit on suppression requires weaker regularization-levels for smaller lags. For example, for the source-based estimator (see Fig 4A) not exceeding a suppression level of 0.3 requires a regularization-level of at most λ = 0, if the lag is 70 degrees. However, for instantaneous interactions, it requires a regularization-level of at most λ = −1.5. In practice this means that the detection of interactions with small lags requires data that is less noisy.

Fig 4D shows the suppression levels of the projected sensor noise for all three estimators. The projected noise term is independent of system activity and its suppression level is therefore independent of the lag. Observe that for weak and strong regularization-levels the noise is completely suppressed, implying that for these levels, the projected noise term is contained in the leakage subspace. For intermediate regularization-levels, the noise term is somewhat less suppressed, particularly in the case of the sensor-based estimator with correction applied in sensor-space. For the other two estimators, the noise is always suppressed stronger than the interaction term is, hence correction in source-space increases the signal-to-noise ratio of the estimated source-space spectral matrix.

Bias reduction

Here we assess to what extent the biases of the estimators are reduced by the leakage corrections. We here understand “bias” to be any function of the to-be-estimated parameter and the expected value of the estimator. Different functions capture different aspects of the estimators’ expected value, some of which might be more relevant more a particular application than others. We quantified the bias by the ratio of the estimators’ absolute expected value, averaged over all but the true interaction-pair, to itself plus the absolute expected value at the true interaction-pair (k0, l0). In calculating the numerator, pairs of the form (j, j) were excluded since these pertain to source power and not to source interaction. Furthermore, to take into account the limited resolution of the inverse operators, we relaxed the definition of “true interaction-pair” to include the eight surrounding pairs (k0±1, l0±1). The bias takes on values between zero and one, where zero corresponds to the estimators’ expected value being zero at all but the true interaction-pair and one corresponds to the estimators’ expected value being zero at the true interaction-pair. We computed the bias as a function of noise-level. Below, we refer to “performance” as one minus the bias. As in the previous section, the projection ranks of the corrected estimators were set to their maximal values so that signal leakage is completely suppressed. To isolate the effect of the correction from the issue of selecting an appropriate regularization-level, we let the regularization-level range between −8 and 0 in steps of 0.5 and selected the value that maximized performance. This was done for both the uncorrected and corrected estimators. The interaction lag was set to ϕ = 0 degrees and all other parameter values were chosen as in Table 1.

Fig 5A (blue and red bars, respectively) shows the performance of the uncorrected source- and sensor-based estimators as a function of noise-level. Performance decreases with increasing noise-level, which is to be expected since the norm of the projected sensor noise is proportional to the noise-level. Also notice that the performance of the sensor-based estimator decreases slower with increasing noise-level than that of the source-based estimator. The reason for this lies in the different filtering properties of the respective inverse operators that are amplified for increasing regularization levels. These differences cause the sensor-based estimator to have a somewhat higher spatial resolution than the source-based estimator. This is illustrated in Fig 5B and 5C, which show the expected values of the source- and sensor-based estimators, respectively, for σ = 0.05 and λ = 10−2.2. The true interaction-pair is denoted by the white circles. Note that the interaction is a bit better visible in the sensor-based estimator, which has a somewhat higher resolution. In particular, the reconstructed sources and their interactions are a bit more localized than in case of the source-based estimator. On the other hand, the sensor-based reconstruction seems to suffer a bit more from artifacts than the source-based estimator, which could affect statistical significance testing. Lastly, we observe that the biases of both estimators are relatively large. For example, a bias of 0.5 means that the average amplitude of the reconstruction at the true source-pair and that at other source-pairs are equal. In practice this means that the interaction will not be detected by pair-wise tests. The main reason for the high biases is the presence of signal leakage around the seed locations. The performance of the corrected estimators are shown in Fig 5D. Observe that the performance of all three corrected estimators is several times higher than that of both uncorrected estimators. The differences are largest in the absence of noise, but can present across noise-levels. To illustrate the differences, we calculated the expectation values of the corrected estimators for σ = 0.05 and λ = 10−3. They are shown in Fig 5 (panels E, F, and G). The interaction stands out clearly, due to the complete suppression of signal leakage.

Fig 5

Bias-reduction through leakage correction.

A. Performance on the test-model of the uncorrected estimators as a function of noise-level. Performance is defined as one minus the bias and ranges between zero and one. B. Expected value of the source-based estimator for σ = 0.05 and λ = −2.2. C. Same as B but for the sensor-based estimator. D. Performance on the test-model of the corrected estimators as a function of noise-level. E. Expected value of the corrected source-based estimator for σ = 0.05 and λ = −3. F. Same as E but for the sensor-based estimator with correction in source-space. G. Same as E but for the sensor-based estimator with correction in sensor-space. In panels B, C, E, F, and G, the true interaction-pair is designated by the white-circles.

Detection power

To assess the ability of the corrected source- and sensor-based estimators to detect system interactions, we used them as statistics for testing the null-hypothesis of no interaction (see ). We thus computed the sensitivity matrices of the test-statistics , , and corresponding to the source- and sensor-based estimators, where k denotes the projection rank, which was set either to k = 0, corresponding to the uncorrected estimators, or to its maximal value k = 21. Detection power was quantified as where M denotes the number of ordered interaction-pairs whose sensitivity is lower than the sensitivity at the true interaction-pair (k0, l0) and p is the number of sensors. Detection power takes on values between 1/(p(p − 1)/2 ≈ 0.05 and 1. If the sensitivity at the true interaction-pair is the highest of all p(p − 1)/2 ordered source-pairs, M = p(p − 1)/2 − 1 and detection power is 1 and if the sensitivity at the true source-pair is the lowest, M = 0 and detection power is 1/(p(p − 1)/2. The significance-level for computing the sensitivity matrix was set to α = 0.025 for the real and imaginary test-statistics and to α = 0.05 for the complex test-statistic. Thus, when carrying out independent and uncorrected pairs-wise tests for significant interactions across all ordered pairs, the significant pairs (at the chosen significance level) correspond to the pairs with non-zero detection-power. Furthermore, the number of samples was set to N = 100, the noise-level to σ = 0.01, and the regularization-level to λ = 10−2. The interaction lag was varied between 0 and 90 degrees. All other parameters were chosen as in Table 1.

Fig 6A (red curve) shows the detection power of the real part of the uncorrected source-based test-statistic as a function of lag. It is 0 for lags larger than about 65 degrees, which means that such interactions remain undetected. Between 0 and 60 degrees, its sensitivity gradually decreases. Furthermore, the detection power is at most about 0.8, which is rather low, for it means that 20% of the false positive interaction pairs have lower p-values than the true interaction pair. Fig 6B and 6C show that the real part of the uncorrected sensor-based test-statistic essentially behaves the same, although the detection power is higher (up to about 0.9). Fig 6D, 6E and 6F (red curves) show the detection power of the imaginary part of the uncorrected test-statistics. The red curves are hidden behind the green curves because the imaginary part of the test-statistics is invariant under the correction projection. Note that the imaginary parts only detect interactions with lags larger than about 40 degrees. However, for larger lags, the detection power is 1, hence the true interaction-pair has the smallest p-value of all significant interactions. Thus, for lags larger than 40 degrees, the imaginary part of both the source- and sensor-based test-statistics will detect the true interaction with probability 100 × (1 − α/2). Fig 6G, 6H and 6I (red curves) show the detection power of the complex-valued uncorrected test-statistics. Its sensitivity is high for all lags (at least 0.8), particularly for the sensor-based test-statistic. That is to say, the complex-valued test-statistic is lag invariant.

Fig 6

Detection power.

A. Detection power of the real uncorrected (red) and corrected (green) source-based test-statistic. D. Detection power of the imaginary uncorrected (red) and corrected (green) source-based test-statistic. G. Detection power of the complex uncorrected (red) and corrected (green) source-based test-statistic. Panels B, E, and H: Same format as panels A, D, and G, respectively, but for the sensor-based test-statistic with correction in source-space. Panels C, F, and I: Same format as panels A, D, and G, respectively, but for the sensor-based test-statistic with correction in sensor-space.

The detection power of the respective corrected test-statistics is also shown in Fig 6 (green curves). The real part has detection power 1 for lags up to about 50 degrees and is 0 for larger lags. The detection power of the imaginary part roughly behaves in the opposite way; it is 0 for lags up to about 30 degrees and is 1 for larger lags. The real and imaginary parts are hence complementary. The complex-valued test-statistics, on the other hand, have detection power 1 for all lags. We conclude that the corrected complex-valued test-statistics detect the true system interaction with probability 100 × (1 − α) irrespective of the lag.

Reconstruction of functional networks

In the previous sections we have considered the case of two active and interacting sources. In practical applications, however, more than two sources might be active and interacting, thus forming a functional network. In this section we consider the detection power of the uncorrected and corrected estimators to detect interactions in functional networks and how detection power depends on the strength of the interactions. Since the difference in performance between the three correction methods is much smaller than that between the uncorrected and corrected estimators and because the sensor-based estimators performed slightly better than the source-based estimators, we restrict the analysis to the sensor-based estimator with correction in source-space and concentrate on the effect of the correction. We considered a network of K evenly-spaced sources that are ordered from left to right in the one-dimensional source-space (see Section Benchmark model). The network structure is characterized by the spectral matrix vv†, where is defined as with ϕ = 2π(n − 1)/K, for n = 1, ⋯, K. The diagonal entries of the spectral matrix are set to 1 to that γ is the shared coherence between the sources. Note that the phases ϕ1, ⋯, ϕ are evenly-spaced on the unit circle and that ϕ1 = 0 and ϕ = 2π − ϕ2 and hence model a traveling wave of activity that propagates from left to right and has wavelength equal to the distance between the first and the last source. This model thus comprises a fully-connected network with K different delays that are spatially organized as a propagating wave. The number of sources was set to K = 4. The network structure is shown in Fig 7A, 7B and 7C, which display, respectively, the real and imaginary parts and the absolute value of the spectral matrix. For this illustration γ was set to 0.5. The absolute value of the spectral matrix (right panel) shows a fully-connected network of K = 4 sources. The K = 4, the delays between neighboring sources equals π/4. Thus, the real part of the spectral matrix (left panel) only shows the interaction between the first and the third source and between the second and the fourth source and the imaginary part of the spectral matrix (middle panel) only shows the interaction between all neighboring sources and, in addition, between the first and the third source. Thus, the real and imaginary parts “see” different interactions and only the absolute value shows the full network topology.

Fig 7

Reconstruction of functional networks.

A. Real part of the true spectral matrix. B. Imaginary part of the true spectral matrix. C. Absolute value of the true spectral matrix. D. Sensitivity matrices of the real (top row), imaginary (middle row), and complex (bottom row) uncorrected sensor-based test-statistic as a function of interaction strength γ ranging from 0 to 1 in steps in 0.1. E. Same as D. but for the corrected test-statistics.

The number of samples was set to N = 100, the noise-level to σ = 0.01 and the regularization-level to λ = 10−2. All other parameters, except the interaction parameters, were chosen as in Table 1. The interaction strength γ ranged from 0 to 1 in steps of 0.1. A complete power analysis as carried out in the previous section involves conducting a separate analysis for each pair of sources, which is cumbersome and not very instructive. It is more illustrative to display the actual sensitivity matrices. The sensitivity matrices for the uncorrected test-statistics are shown in Fig 7D. We make two observations. First, in the absence of interaction (i.e. γ = 0), the real and complex test-statistics have high sensitivity at and around the true source locations. When statistical tests are carried out this leads to false positives around the source locations. These spurious interactions have been referred to as first-order [37] and are due to the mere presence of active sources. No spurious interactions are present in the imaginary statistic because first-order spurious interactions are instantaneous. Second, when the interaction strength in increased, the sensitivity of all three test-statistics increases at and around the true interactions. Also the complex test-statistic, however, is sensitive to both real and imaginary interactions. The spurious interactions surrounding the true interactions have been referred to as second-order [37] or ghost interactions [38, 39]. They are independent of source-power and are due to the existence of true interactions. Note that all three test-statistics are equally affected by second-order spurious interactions. Incidentally, the local nature of second-order spurious interactions has been exploited in [39] to mitigate their effects. Fig 6E shows the sensitivity matrices of the corrected test-statistics. As in the previous sections, the projection rank was set to its maximal value of k = 21 so that the leakage term is completely suppressed. The figure shows that, in the absence of interaction, the sensitivity of the real and complex test-statistics is now zero hence no spurious interactions will be detected. Furthermore, when the interaction strength is increased, the sensitivity of all three test-statistic increases at and around the true interactions. Again, only the complex test-statistic is sensitive to both real and imaginary interactions. We conclude that, in this simple scenario, the corrected complex test-statistic is sensitive for all interactions, i.e. irrespective of their lag, and does not suffer from first-order spurious interactions. A price that is paid for the correction, however, are the artifacts in between the true source interactions, which appear as yellow dots (see right most panel in the top row of Fig 6E) and which are a consequence of the interaction term not being orthogonal to the leakage term (see ).

Discussion and conclusions

Summary and relevance

Detection of functional interactions in linearly-mixed systems is complicated by signal leakage, which refers to the incomplete unmixing of source signals due to the ill-posedness of the inverse problem under study [6]. If the forward model is (nearly) instantaneous, such as the forward models used for EEG/MEG [5] and local field potential (LFP) recordings [3], signal leakage is usually dealt with by discarding all instantaneous interactions, whether true or spurious. In practice, this is done by discarding the real part of complex-valued test-statistics and only retaining their imaginary part [11-15]. The drawback of this approach is that the sensitivity of the resulting test-statistics decreases with decreasing interaction lag, leading to false negatives and distortions in the topology of reconstructed functional networks [2]. More recent approaches exploit the relationship between signal leakage and the algebraic structure of linear inverse operators [24, 25, 27]. In [27] signal leakage is suppressed by correcting the sensor-space spectral matrix and subsequently projecting it to source-space using a non-linear dipole scanning algorithm [28]. To be applicable to more general source configurations, however, the correction method needs to be combined with more general inverse methods. In the current study we showed how to combine the correction with an arbitrary linear inverse operator and proposed two alternative methods in which the order of correction and projection is swapped, leading to novel sensor- and source-based estimators of system interactions. We demonstrated that all three estimators enable lag-independent detection of system interactions. As such, our study contributes to the development of reconstruction methods for system interactions from linearly-mixed observations and we expect applications in neuroscience and the physical sciences.

The geometry of signal leakage

We have provided a general characterization of signal leakage in source-space as characterized by spectral matrices that applies to all linear inverse operators. In particular, (the expected value of) the sample estimator of the source-space spectral matrix can be decomposed into three parts that are contained in different linear subspaces of the vector space of complex-valued (Hermitian) matrices. The first and second parts only depend on source-power and interactions, respectively, and we have referred to the corresponding subspaces as the leakage and interaction subspaces, respectively. The third part only depends on projected measurement noise and we have referred to the corresponding subspace as the noise subspace. In [27] a similar characterization is described for the sample estimators of the sensor-space spectral matrix. The difference between the source- and sensor-based characterizations of signal leakage is that, in the first case, the leakage and interaction subspaces are defined in terms of the resolution operator associated with a given inverse operator, whereas in the second case, the resolution operator is replaced by the forward matrix. These source- and sensor-space characterizations of signal leakage formalize the notions used in [37-39] in the special case in which connectivity is defined in terms of spectral matrices. In particular, the leakage term formalizes the equivalent notions of “first-order interactions” [37] and “spurious interactions” [38, 39]. The interaction term can be further decomposed into a term corresponding to true connectivity and a term corresponding to a type of spurious interaction that it independent of source power and which has been referred to as “second-order” [37] or “ghost” [38, 39] interactions. The current study focused on suppressing first-order interactions, whereas [38, 39] proposed a method to “bundle” second-order interactions and thereby reducing the number of false positives. However, our characterization of second-order interactions in terms of tensor-products of point-spread functions might be helpful in bundling these interactions.

Source- versus sensor-based estimation of system interactions

We have described how leakage can be suppressed by applying projection operators to estimated source- or sensor-space spectral matrices. The source-space spectral matrices themselves can be estimated in two ways. The most obvious way is to project the observed sensor data to source-space and subsequently estimate the source-space spectral matrix from the reconstructed source-activity. An alternative way is to first estimate the sensor-space spectral matrix and subsequently projecting it to source-space by an inverse operator obtained by vectorizing a bilinear forward model in matrix space. This is essentially what has been proposed in [27]. We provided some insight into this estimator by giving an alternative characterization as the minimizer of a regularized cost function in matrix space. We have derived the filtering coefficients of both estimators and established a natural correspondence, which we exploited in the evaluation of their performance. In particular, the first estimator regularized with level λ can be associated with the second estimator regularized with level λ2. Corresponding filtering coefficients remain somewhat different, however, and this is reflected in the second estimator having a slightly higher spatial resolution.

Leakage suppression in source- and sensor-space

We have formulated and evaluated three methods for suppressing signal leakage in reconstructed system interactions, each of which involves the application of a suitable projection operator. The projection operator can be applied in source-space, to either the source- or sensor-based estimator of the source-space spectral matrix, or it can be applied in sensor-space, the latter corresponds to the correction method proposed in [27]. The sensor-based corrected estimators were more efficient in suppressing spurious interactions than the source-based corrected estimator and correction in source-space was more efficient than correction in sensor-space (Fig 4). The reductions in bias due to the corrections were comparable for the three methods (Fig 5). We subsequently used the corrected estimators as statistics for testing for significant interactions. Application to the model system demonstrated that the corrected test-statistics were able to detect system interactions, irrespective of the lag, with higher power than the uncorrected test-statistics. The detection power was the same for all three corrected test-statistics (Fig 6). Lastly, we showed that the corrected, but not the uncorrected, test-statistics are able to reconstruct functional networks comprising more than two sources (Fig 7). Taken together, we conclude that, out of the three test-statistics, the sensor-based test-statistic with correction in sensor-space is to be preferred, although the differences are small. We thus advice to use this test-statistic in practical applications to EEG/MEG and LFP recordings.

Scope and limitations

This study focussed on theoretical aspects of reconstructing system interactions in linearly-mixed systems and did not address practical aspects. Two of these aspects are the selection of the regularization-level of the inverse operator and the rank of the correction operators. With regard to selecting appropriate values for the regularization parameters, several techniques can be explored including cross-validation [40], L-curve methods [41], residual analysis [42], and empirical Bayes [43]. Which method is most suited for which type of system and forward model is largely an empirical question that can be addressed through realistic numerical simulations. An interesting observation in this context is that the appropriate level of regularization depends on which aspect of the system under study is of interest. Thus, in [35] it was observed that estimation of brain interactions from MEG data requires less strong regularization than the reconstruction of the activity time-courses. In [36] this observation is formally investigated. Selection of the optimal rank of the projection operators remains an open question that can potentially yield interaction estimators with less bias and more powerful statistical hypothesis tests. In [27], the optimal projection rank was chosen on the basis of simulated data with a specific source configuration and parameter settings. For example, two sources were assumed to be active with equal strengths. However, no evidence was given to show that the obtained projection ranks are (close to) optimal for other source configurations and a data-driven selection method would therefore be preferable. The difficulty lies in the fact that the leakage, interaction, and noise terms are not separately observable, but only their superposition is. A data-driven method might be based on the sequence of estimates of the interaction matrices.

To determine to what extent the proposed estimators can successfully be applied to experimental EEG/MEG and LFP data requires characterizing and quantifying the effects of different factors that are present in such recordings such as cortical background activity, measurement noise, and inaccuracies in forward modeling, such as sensor positions, dipole orientations, and conduction properties of the volume conductor. Furthermore, although in this study we have considered functional networks comprising several sources, simulations of larger functional networks are required to obtain a better understanding of which types of network topologies can be successfully reconstructed and which ones are more challenging. Realistic simulations that take into account all of the above effects, however, cannot currently be carried out by standard implementations of the estimators due to its high memory demands. For example, implementation of the source-based estimator for a source-space of p = 104 voxels requires computing the singular value decomposition of the matrix B, whose dimensions are p2 × p, which is not infeasible with standard algorithms. Similar considerations apply to the sensor-based estimator. In [27], these computational problems were circumvented by down-sampling the human cortex by a factor ten. Although this approach is suitable when using adaptive spatial filters for the inverse modeling, they are not appropriate when working with non-adaptive spatial filters as used in the current study, because such filters are unable to actively suppress source-activity from undesired locations. Thus, a natural follow-up study would be to use iterative reconstruction methods in combination with realistic EEG/MEG volume-conductor models and simulated functional networks to characterize the performance of the proposed methods.

(M)

Click here for additional data file.

(M)

Click here for additional data file.

(M)

Click here for additional data file.

(M)

Click here for additional data file.

(M)

Click here for additional data file.

(M)

Click here for additional data file.

(M)

Click here for additional data file.

(M)

Click here for additional data file.

(M)

Click here for additional data file.

(M)

Click here for additional data file.

(M)

Click here for additional data file.

(M)

Click here for additional data file.

(M)

Click here for additional data file.

19 Jun 2020

PONE-D-20-10601

Lag-invariant detection of interactions in spatially-extended systems using linear inverse modeling

PLOS ONE

Dear Dr. Hindriks,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

The reviewer raised important issues to the presentation and evaluation of the whole methodology.

I recommend to answer to all his comments one by one and re-submit the revised manuscript.

Important issues are the description of the formulas in a clear fashion to everyone and also to

submit in a public database your code for further validation from other researchers.

Please submit your revised manuscript by Aug 03 2020 11:59PM. If you will need more time than this to complete your revisions, please reply to this message or contact the journal office at plosone@plos.org. When you're ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

We look forward to receiving your revised manuscript.

Kind regards,

Stavros I. Dimitriadis

Academic Editor

PLOS ONE

Additional Editor Comments:

The reviewer raised important issues to the presentation and evaluation of the whole methodology.

I recommend to answer to all his comments one by one and re-submit the revised manuscript.

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. In your Data Availability statement, you have not specified where the minimal data set underlying the results described in your manuscript can be found. PLOS defines a study's minimal data set as the underlying data used to reach the conclusions drawn in the manuscript and any additional data required to replicate the reported study findings in their entirety. All PLOS journals require that the minimal data set be made fully available. For more information about our data policy, please see http://journals.plos.org/plosone/s/data-availability.

Upon re-submitting your revised manuscript, please upload your study’s minimal underlying data set as either Supporting Information files or to a stable, public repository and include the relevant URLs, DOIs, or accession numbers within your revised cover letter. For a list of acceptable repositories, please see http://journals.plos.org/plosone/s/data-availability#loc-recommended-repositories. Any potentially identifying patient information must be fully anonymized.

Important: If there are ethical or legal restrictions to sharing your data publicly, please explain these restrictions in detail. Please see our guidelines for more information on what we consider unacceptable restrictions to publicly sharing data: http://journals.plos.org/plosone/s/data-availability#loc-unacceptable-data-access-restrictions. Note that it is not acceptable for the authors to be the sole named individuals responsible for ensuring data access.

We will update your Data Availability statement to reflect the information you provide in your cover letter.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Partly

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: In this paper, the author proposes a novel method to address the important problem of estimating phase-synchronization while accounting for spurious interactions due to MEG field spread or EEG volume conduction. The method is based on a subspace projection to suppress the contribution of source power to the sensor-space spectral matrix. In contrast to the original paper proposing this approach (Ossadtchi et al. (2018)), this paper demonstrates how it can be combined with any linear inverse operator to obtain a lag-independent estimator of system interactions. The method is novel and promising and the paper is well-written. However, a number of concerns should be resolved before the paper could be considered for publication. Please see attached for details.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Nitin Williams

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

Submitted filename: Hindriks_review.pdf

Click here for additional data file.

21 Sep 2020

Below, the reviewer’s comments are displayed in black, authors responses in blue, and citations from the revised manuscript in green.

In this paper, the author proposes a novel method to address the important problem of estimating phase-synchronization while accounting for spurious interactions due to MEG field spread or EEG volume conduction. The method is based on a subspace projection to suppress the contribution of source power to the sensor-space spectral matrix. In contrast to the original paper proposing this approach (Ossadtchi et al. (2018)), this paper demonstrates how it can be combined with any linear inverse operator to obtain a lag-independent estimator of system interactions. The method is novel and promising and the paper is well-written. However, a number of concerns should be resolved before the paper could be considered for publication:

Major concerns:

1. The method separates signal leakage-related and interaction-related contributions to sensor-space spectral matrix. However, it has been shown that the effect of signal leakage depends on the underlying true interaction strength (Fig. 3 in Palva et al. (2018)). In that sense, there is a portion of variance in the sensor-space spectral matrix that cannot be attributed to mere signal leakage or mere interaction but is due to the combination of the signal-leakage and underlying interaction strength. Does the method account for this?

For example, consider interactions between two different pairs of signals. The source-power for both pairs of signals are the same, while source-coherence for the two pairs is different. Due to the differences in source-coherence, the effect of ‘signal leakage’ will also be different. Now, would this difference be reflected in the ‘signal leakage’ term? Alternatively, does this constitute a case where the ‘source power’ and ‘source interactions’ do not lie in different linear subspaces, as the method assumes?

When the interaction-strength is varied while keeping source-variance fixed (i.e. when source-coherence is varied), the corrected leakage term remains unchanged and the interaction term is scaled by the same factor. That is so say, the effect of source-coherence is a different relative weighting of the corrected leakage and interaction terms. I have added a formal description of this effect to the manuscript (see Section Suppression in source-space).

On a related note, I did not quite understand how the first term in equation 5 can have non-zero non-diagonal terms when in the previous line, it is mentioned the M can be split into M+ and M-, where M+ has all off-diagonal terms set to zero and M- has all diagonal terms set to zero. Given M+ has all off-diagonal terms set to zero, would we not expect the first term of equation 5 to also have all off-diagonal terms set to zero?

In matrixform (i.e. non-vectorized form), the first term of Equation (5) is given by R*S^+*R^T, with R the resolution matrix. Although S^+ is indeed a diagonal matrix, R*S^+*R^T not necessarily is (the product of a diagonal matrix with two non-diagonal matrices is not necessarily a diagonal matrix). In fact, the non-diagonal entries of R*S^+*R^T correspond to first-order leakage.

2. Related to point 1 above, how does the method’s performance compare to obtaining corrected interaction strength by subtracting estimated interaction strength from interaction strength obtained from a surrogate dataset? The surrogate dataset is such that it has also been subject to MEG field spread, its ‘source power’ is the same as the original dataset but ‘source interaction’ is set to zero. Would the correction by the proposed method be essentially equivalent to this? It would be important to compare the method’s performance with correction based on such a surrogate dataset.

The reviewer raises an interesting point here. If the source power were known it could in principle be used to remove first-order leakage without suppressing second-order leakage. I tried this approach and I found that the estimated source power is much lower than the true power, which is because the power is smeared out over many locations due to the squared-norm penalty in the penalized least squares problem that is used to obtained the inverse operators. As a consequence, the correction has practically no effect. Because of this negative result I decided to not include it in the manuscript.

3. The simulations based on two nodes are instructive and provide a basic demonstration of the method. However, since the assessment of the method is exclusively based on simulations, additional simulations should be done with more realistic simulated data. In particular, the correction should be applied to infer functional networks with multiple nodes. Between 80 and 120 nodes would make the network similar to those often estimated in MEG. The ‘true’ network should have interaction strengths ranging from 0 to 1 and the ability of the method to recover these strengths should be assessed. The amplitude and spectral properties of time series should also be similar to MEG data. For this, models of region dynamics such as Neural Mass Models (NMMs) (Jansen & Rit (1995)) could be used. Since properties of data and dimensionality of such a model would be similar to that encountered with MEG, good performance of the method on such a dataset would recommend its application to MEG.

The reviewer raises an important point here. In the revised manuscript I have added a section named Reconstruction of functional networks in which the performance of the methods is assessed in reconstructing functional networks comprising several sources. Network dynamics is designed so that the interaction latencies cover the entire circle [0,360] degrees. Performance was assessed for interaction strengths ranging between 0 and 1.

I would like to clarify that the goal of the simulations is not to see how well the methods will perform on experimental EEG/MEG data, but to introduce the methods and to provide a proof-of-principle. I agree with the reviewer that realistic simulations need to be performed in order to assess the performance on experimental EEG/MEG data. Besides using a different forward model, this requires the inclusion of many more aspects such as cortical background activity, sensor noise, inaccuracies in the forward modeling, etc. Furthermore, the methods are computationally too expensive to be implemented on a high-resolution source-space. I discuss these issues in the discussion section of the revised manuscript (see the last paragraph):

To determine to what extent the proposed estimators can successfully be applied to experimental EEG/MEG and LFP data requires characterizing and quantifying the effects of different factors that are present in such recordings such as cortical background activity, measurement noise, and inaccuracies in forward modeling, such as sensor positions, dipole orientations, and conduction properties of the volume conductor. Furthermore, although in this study we have considered functional networks comprising several sources, simulations of larger functional networks are required to obtain a better understanding of which types of network topologies can be successfully reconstructed and which ones are more challenging. Realistic simulations that take into account all of the above effects, however, cannot currently be carried out by standard implementations of the estimators due to its high memory demands. For example, implementation of the source-based estimator for a source-space of p = 10^4 voxels requires computing the singular value decomposition of the matrix B_j, whose dimensions are p^2\\times p, which is not infeasible with standard algorithms. Similar considerations apply to the sensor-based estimator. In [Ossadtchi, 2018], these computational problems were circumvented by down-sampling the human cortex by a factor ten. Although this approach is suitable when using adaptive spatial filters for the inverse modeling, they are not appropriate when working with non-adaptive spatial filters as used in the current study, because such filters are unable to actively suppress source-activity from undesired locations. Thus, a natural follow-up study would be to use iterative reconstruction methods in combination with realistic EEG/MEG volume-conductor models and simulated functional networks to characterize the performance of the proposed methods.

The use of neural mass models makes it more difficult to assess the performance of the methods proper and to disentangle it from the effects of pre-processing (e.g. bandpass filtering). This is because what effectively matters for performance is the distribution of the time-frequency coefficients of the sensor signals. Thus, by directly using time-frequency coefficients, performance can be disentangled from indirect issues. In fact, the time-frequency coefficients of linearized neural mass models with Gaussian white-noise input are independent and normally-distributed, just as in the manuscript, so that, as far as performance is concerned, the results are the same. I do acknowledge, however, that such models can be appropriate for testing entire processing pipelines as they are used in practice.

4. The paper describes three variants of the proposed method. For a user of the method, are there any technical reasons for choosing one variant over another (apart from simulation results)? Does one method perform better than others in all cases? If not, in what situations should one use one variant over another? Clear description of technical differences between variants and guidance on when to use one variant over another would be useful. Alternatively, the paper could focus on a single variant if one variant performs best in all cases or if there are only small differences in performance between the variants.

For completeness, I prefer to keep all three methods. I added a paragraph in the discussion section I which I compare the performance of the three methods and draw a conclusion about which one to use in practice:

The sensor-based corrected estimators were more efficient in suppressing spurious interactions than the source-based corrected estimator and correction in source-space was more efficient than correction in sensor-space (Figure 4). The reductions in bias due to the corrections were comparable for the three methods (Figure 5). We subsequently used the corrected estimators as statistics for testing for significant interactions. Application to the model system demonstrated that the corrected test-statistics were able to detect system interactions, irrespective of the lag, with higher power than the uncorrected test-statistics. The detection power was the same for all three corrected test-statistics (Figure 6). Lastly, we showed that the corrected, but not the uncorrected, test-statistics are able to reconstruct functional networks comprising more than two sources (Figure 7). Taken together, we conclude that, out of the three test-statistics, the sensor-based test-statistic with correction in sensor-space is to be preferred, although the differences are small. We thus advice to use this test-statistic in practical applications to EEG/MEG and LFP recordings.

5. All simulations were done with single examples. Where applicable, it would be good to perform simulations with multiple examples so that variability of method’s performance can be assessed. In this regard, particularly for Figs 3-5, it would be good to see confidence intervals marked for all line plots.

Sample variability is assessed in the sections Detection Power and Reconstruction of functional networks. In these sections, the statistical power of the different test-statistics is approximated using Monte Carlo simulations. Specifically, the calculation of each point (on each curve) in Figure 6 of the revised manuscript (which corresponds to Figure 5 of the unrevised manuscript) is based on a large number of realizations of the generative model (10^3) that are used to approximate the statistics’ sample distributions and confidence regions (see Section Detection of system interactions). Thus, confidence regions are already incorporated in calculating these curves. Sections Effectiveness and Bias reduction cover performance aspects that have no inherent variability but relate to asymptotic properties of the estimators (effectiveness and bias).

6. Early in the paper, it would be good to have a simple diagram illustrating the method’s operation, to complement the mathematical description.

Thanks for the suggestion: This will indeed be helpful to the reader. I added a section named Subspace projections for suppression of spurious interactions that contains a figure illustrating the three different correction methods (see Figure 3 of the revised manuscript).

7. In the Discussion, it would be good to describe any limitations of the method and particularly, in what practical situations might assumptions of the method be violated. Under what situations might errors be introduced into the estimation of corrected interaction strength (for e.g. low SNR? low lags?) How sensitive is the method to errors in the forward model? Further, how might the method work differently for EEG and for MEG data, and why? Some mention of these points is made while describing the method, but a fuller treatment in Discussion would be good. Any intuition on strengths, assumptions and limitations of the method would be very useful to future users.

The reviewer raises important questions relating to the application of the proposed methods to experimental EEG/MEG data. The manuscript focuses on the formulation and basic properties of the methods and I acknowledge that in order to be applicable to experimental EEG/MEG data, several questions needs to be addressed. These questions are ideally addressed in a separate study that entirely focuses on the application of the method to EEG/MEG data. In the discussion section I added a paragraph in which I discuss these issues (see cited text in response to comment 3.

8. A MATLAB and/or Python implementation of the proposed method could be made available to the community, in keeping with principles of Open Science.

I fully agree and I will upload the documented Matlab code to a public repository.

Minor concerns:

1.Some typographical errors should be corrected: Page 3, final paragraph: “varies applications” should be “various applications” Page 10, final paragraph: “die to the mixing” should be “due to the mixing”

Corrected.

2.I might be mistaken, but in Page 17, would it not be that: T_imag=Imag(gamma(k,l))=abs(gamma(k,l))sin(theta(k,l)) , i.e. rather cos and T_real=Real(gamma(k,l))=abs(gamma(k,l))cos(theta(k,l)), i.e. rather than sin

Yes indeed. Thanks for spotting these errors.

3.Fig. 1B - z-axis labels should be specified

Done. I have referred to the values on the color-axis as “arbitrary units” (a.u.) because I have set the tissue conductivity in Poisson’s equation to the arbitrary value of 1.

4.Fig. 2A - colour of lines for ‘real’ and ‘imag’ too similar, make colour of lines different and include legend

Fig 2B - line colours and marker types too similar, change line colours and include legend

Fig 2C - white dots showing true interaction could be change to bigger black dots for visibility, z-axis labels should be specified Fig. 2 caption - “strength of interaction wat set” should be “strength of interaction was set”

Done.

5.Fig. 3 - should include legend for black lines in A, B and C!

Done.

6.Fig. 4 - white dots showing true interactions could be bigger black dots for visibility, z-axis labels of images should be specified, colorbar scales could be equal to aid comparisons

Done. I left scaling of the color bars as it was because the overall scale is unimportant for performance: It’s the relative strength, i.e. across points in source-space, that determines performance and this is better visible now.

7.Fig. 5 - x-axis labels could be 10 to 90, since 90 degrees is often mentioned in text. Panels A,C,E should be labelled ‘Source-based’ and changed to A,B,C. Panels B, D, F should be labelled ‘Sensor-based’ and changed to D,E,F.

Done.

Submitted filename: Response to Reviewers.pdf

Click here for additional data file.

9 Nov 2020

Lag-invariant detection of interactions in spatially-extended systems using linear inverse modeling

PONE-D-20-10601R1

Dear Dr. Hindriks,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Stavros I. Dimitriadis

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Dear Author

After carefully reading the revised manuscript and your response to the reviewer's comments,

I am glad to announce you that your paper has been accepted for publication in its

current form.

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: No

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: I thank the author for addressing my comments and commend the author on methodological development which is a valuable contribution to the field. A couple of typographical errors I spotted:

1. Figure 5, caption: In panels B, C, 507E, F, and G, the true interaction-pair is designated by the "white-circles." Should be "black circles"

2. Discussion, final paragraph: "not infeasible with standard algorithms". Should this be "infeasible with standard algorithms"?

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Nitin Williams

24 Nov 2020

PONE-D-20-10601R1

Lag-invariant detection of interactions in spatially-extended systems using linear inverse modeling

Dear Dr. Hindriks:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Stavros I. Dimitriadis

Academic Editor

PLOS ONE