Identifying oscillatory brain networks with hidden Gaussian graphical spectral models of MEEG

Identifying the functional networks underpinning indirectly observed processes poses an inverse problem for neurosciences or other fields. A solution of such inverse problems estimates as a first step the activity emerging within functional networks from EEG or MEG data. These EEG or MEG estimates are a direct reflection of functional brain network activity with a temporal resolution that no other in vivo neuroimage may provide. A second step estimating functional connectivity from such activity pseudodata unveil the oscillatory brain networks that strongly correlate with all cognition and behavior. Simulations of such MEG or EEG inverse problem also reveal estimation errors of the functional connectivity determined by any of the state-of-the-art inverse solutions. We disclose a significant cause of estimation errors originating from misspecification of the functional network model incorporated into either inverse solution steps. We introduce the Bayesian identification of a Hidden Gaussian Graphical Spectral (HIGGS) model specifying such oscillatory brain networks model. In human EEG alpha rhythm simulations, the estimation errors measured as ROC performance do not surpass 2% in our HIGGS inverse solution and reach 20% in state-of-the-art methods. Macaque simultaneous EEG/ECoG recordings provide experimental confirmation for our results with 1/3 times larger congruence according to Riemannian distances than state-of-the-art methods.


Introduction
In vivo identification of functional brain networks would greatly benefit from the inverse solutions for MEG/EEG electromagnetic observation modalities (Larson-Prior et al. 2013;Lopes da Silva 2013).A direct relationship to the activity of postsynaptic potentials (PSPs) is attributed to the latent brain variables generating the MEG/EEG (Paul L. Nunez 1974;Freeman 1975;Valdes et al. 1999).
Where the actual latent brain variables are local currents, vector process () in time-domain  with components representing locally summated PSP myriads within a spatial scale of millimeters (Paul L Nunez and Srinivasan 2006).See appendices I. Notation for variables and mathematical operations and II.Nomenclature for theoretical quantities.MEG/EEG observations () with millisecond temporal resolution are due to these latent brain variables (or activity) () converted according to a linear and stationary forward-operator   or Lead Field summarizing an Electromagnetic Forward Model (EFM) (Jörge Javier Riera and Fuentes 1998;Baillet, Mosher, and Leahy 2001;Grech et al. 2008).An inverse solution from data () provides estimates () for this latent brain activity within a natural timescale (Hämäläinen and Ilmoniemi 1994).Estimates () could therefore help identify functional brain networks in association with the actual PSP activity and within very small timescales, avoiding timedomain distortions (He et al. 2019;Reid et al. 2019).The target of all identification is estimating functional connectivity within a given functional network model interpreted as the set of parameters   governing coactivation in (), which is without loss of generality represented by a multivariate probability (()|  ) (K. Friston 2011;Valdes-Sosa et al. 2011).Therefore, inverse solutions estimating () are only a first step to subsequently estimate latent functional connectivity   (a multistep procedure).Estimation  ̂ commonly requires a second-step inverse solution from such brain activity pseudodata () given the functional network model (()|  ) (Razi and Friston 2016;Lee, Friston, and Horwitz 2006;Valdés-Sosa et al. 2005).Functional brain networks synchronized at specific frequencies or bands of frequency (Fig. 1a) are thought to underpin the periodic components or oscillations observed in the MEG/EEG (Niedermeyer and da Silva 2005).These oscillations are a spectral vector process (, ) at a frequency  obtained from the narrow band filtered, or alternatively Hilbert transform of the narrow band filtered vector process ().Oscillations (, ) then reflect analogous latent brain oscillations (, ) (Bruns 2004;Nolte et al. 2020;Hipp et al. 2012) according to the forward operator   of the EFM spectral equivalent.These so-called oscillatory brain networks are described by frequency-specific functional connectivity   () , which produces brain oscillations (, ) even within a millisecond timescale  (Moffett et al. 2017).Specific oscillatory network patterns exhibit a transient behavior which may be consider stationary in a timescale of approximately seconds (Engel, Fries, and Singer 2001;Varela et al. 2001;Vidaurre, Abeysuriya, et al. 2018;Vidaurre, Hunt, et al. 2018;Tewarie et al. 2019).

Bayesian inverse solution as a maximum a posteriori probability (MAP)
An inverse solution (equation 1-1) estimating the latent variables (categories  and ) from their data (cathegory ) is without loss of generality a Bayesian maximum a posteriori probability (MAP) (Bishop 1995).A MAP  ̂ is the optimum value or estimate for an a posteriori probability (|) computed for a given likelihood (|) explaining  upon  and positing an a priori probability () upon .The nominal Bayesian procedure (Mackay 1999) is a first-type MAP (MAP1) involving only a given first-type likelihood (|).The full Bayesian procedure (MacKay 2003) is a second-type MAP (MAP2) involving a second-type likelihood (|), which must be determined marginalizing a given first-type likelihood (|) explaining  upon (parameters) .An a priori probability (|) explains these parameters upon (hyperparameters) .
̂  The expected covariance matrix  ̂ () () (equations 1-5) is determined from an ensemble covariance matrix  ̂ () (), and its sampled estimator  ̂ () () in the expectation stage MAP1 (Dempster, Laird, and Rubin 1977) is also common for a first-step MAP1 (Mattout et al. 2006;K. Friston et al. 2008;Trujillo-Barreto, Aubert-Vázquez, and Valdés-Sosa 2004 A proof of concept for the HIGGS model in EEG simulations reveals that multistep inverse solutions (equations 1-2, 1-3) may produce a very distorted precision matrix  ̂ () in a second step MAP1 (equations 1-3).When even no ill-conditioning holds (ruling out leakage) in a first step MAP1 (equation 1-2) determining the brain oscillations (, ).In the same situation, the onestep inverse solution (equations 1-5, 1-6, 1-7) determines  ̂ () exactly.We employ the sparse hermitian graphical a priori model to determine  ̂ () via the unbiased hgLASSO algorithm for both multistep (equations 1-3) and onestep (equations 1-5) inverse solutions.Implementations of onestep inverse solutions that violate such unbiasedness conditions with other hermitian graphical a priori models also outperform multistep hgLASSO inverse solutions.The multistep procedure is then the systemic cause of distortions irremediable even with an exact second-step inverse solution for  ̂ ().This proof of concept also reveals that multistep inverse solutions employing a sparse real-valued graphical a priori model or Graphical LASSO (gLASSO) that approximates the HIGGS model further distorts  ̂ ().We confirm this difference in performance between multistep and onestep procedures by comparing EEG inverse solutions against their fine-grained ECoG reference inverse solutions in macaque simultaneous recordings.

GGS model and interpretations of functional connectivity
The PSP processes myriad fulfills some frequency-domain mixing conditions and consequently Gaussianity in the complex-valued follows for brain oscillations (, ) (Brillinger 1975;Rosenblatt 1956;Roweis and Ghahramani 1999;Freeman 1975).Then, without loss of generality, the Gaussian Graphical Spectral (GGS) model (equation 2-1) with complex-valued hermitian precision matrix   () is the functional network model for brain oscillations.This GGS model is valid at any frequency  within the human spectrum ( ∈ ) from 0.5 to 140 hertz (Niedermeyer and da Silva 2005;Le Van Quyen and Bragin 2007;Frauscher et al. 2018) and within a resting state or task transient (∀ ∈ ) from a block design of approximately 2 seconds (Engel, Fries, and Singer 2001;Varela et al. 2001;Vidaurre, Abeysuriya, et al. 2018;Vidaurre, Hunt, et al. 2018;Tewarie et al. 2019).
a This matrix may also be represent a time-frequency connectivity behaivor   ( 0 , ) for a transient  = ( 0 ,  0 + Δ) in a sliding in a slower timescale  0 .This timescale  0 is the domain for evolving network configurations as represented across task or resting state transients.The Hermitian Graphical Naïve (hgNaïve) estimator (a priori free P = 0 in equation 2-8) is the inverse of the sample cross-spectral matrix  ̂ () ←  ̂ −1 () and usually yields a quite dense matrix with many spurious connectivities, which suggests the use of different penalizations.Then, we introduce (equations 2-8) the L2 norm P = ‖•‖ 2 2 , the Hermitian Graphical Ridge (hgRidge) (Honorio and Jaakkola 2013; van Wieringen 2019), the L1 norm P = ‖•‖ 1 , the Hermitian Graphical LASSO (hgLASSO) (Tugnait 2019), known as Graphical LASSO (gLASSO), in the real variable (J.Friedman, Hastie, and Tibshirani 2008).In this paper, we shall emphasize that hgLASSO is the a priori model proposed to produce unbiased sparse estimates with the optimum of the target function (equations 2-5), where hgNaïve and hgRidge as well as gLASSO are hereinafter model violations.The critical issues are stable and scalable hgLASSO calculations in ultrahigh matrix dimensions preserving unbiasedness.We introduce a novel algorithm reformulating MAP1 (equations 2-7) with Bayesian hierarchical hgLASSO.This algorithm is plugged into multistep, and onestep methods are deferred to avoid interrupting the flow of exposition 2.3 hermitian graphical LASSO (hgLASSO).

HIGGS model underneath the Gaussian EFM spectral equivalent
As with latent brain activity (), the conversion of brain oscillations (, ) is given according to the same Electromagnetic Forward Model (EFM) (Jörge Javier Riera and Fuentes 1998;Baillet, Mosher, and Leahy 2001;Grech et al. 2008).The EFM forward operator   (equation 2-11) describes in the time domain (∀ ∈ ) and spectral domain (∀ ∈ ) a purely linear and stationary measurement process.These measurements possess no relation whatsoever with any other biological mechanism b and are only perturbed by a spectral noise process (, ) at the sensors.
The EM approximated MAP2 (maximization stage)  ̂ (+1) () and  ̂ (+1) () (equation 2-25) may be trapped at local optima but is explicitly Wishart and scalable to high dimensions combined with regularization priors (  ()) and  (  ()) that may facilitate the search closest to the global optimal (McLachlan and Krishnan 2007;Wills, Ninness, and Gibson 2009).As with the previous firsttype likelihood for   () (equations 2-6, 2-7, 2-8), any of the penalizations discussed in the previous section may be used here to regularize the EM procedure.See Lemma 1 in IV.HIGGS second-type maximum a posteriori and priors (Corollary to Lemma 1).
For the spectral sensor noise process, we assume that   () is known but a scalar factor   2 () to be estimated:   () =   2 ()  , for which we use an exponential prior with scale parameter   (equations 2-27), with p being the number of MEG/EEG sersors and  being the sample number.
(  2 ()|  ) = (  2 ()|  p) (II-27) In this paper, we define   as the identity matrix, but note that it might be used to encode spurious EEG sensor connectivity due to scalp leakage currents.  (equation 2-27) could be used to encode the instrumental noise inferior threshold that we define as 10% of the EEG signal (equation 2-28).

Unbiasedness of the MAP1 inverse solution with hgLASSO a priori probability
We leverage recent results on the distribution of high-dimensional estimators of precision matrices of Jankova and Van de Geer (JVDG) (Janková and van de Geer 2018) to provide statistical guarantees for the HIGGS inverse solution.First, we reduce the bias of  ̂ () at each iteration of EM maximization in the onestep hgLASSO inverse solution (equations 2-11) or second step in the multistep hgLASSO inverse solution (equations 2-6) by substituting  ̌ () () in  ̂ (), an unbiased estimator via "desparsification" (equation 2-14).This unbiased estimator is the complex-valued extension for the graphical LASSO (gLASSO) (J.Friedman, Hastie, and Tibshirani 2008) exposed in the JVDG theory (Janková and van de Geer 2018).With the fixed value of the regularization parameter   = √ (q) and  ≫ q, whose z-statistic (equation 2-16) possesses a Rayleigh distribution with variance 1 √2 ⁄ .We zero all values of  ̂ (; , ) with (; , ) lower than a threshold to ensure a familywise error of type I.It should be noted that this debiasing and thresholding yields every iteration of the EM maximization in the onestep hgLASSO inverse solution (equations 2-11) or the second step in the multistep hgLASSO inverse solution (equations 2-6) by substituting  ̂ () in  ̂ (), a statistically guaranteed thresholded connectivity matrix.

Scalable and stable MAP1 inverse solution with the hgLASSO algorithm
The solution for the Hermitian Graphical LASSO (hgLASSO) that we present here is a transformation of its prior (equation 2-3) by means of the extension to the complex variable of the scaled Gaussian mixture procedure (equation 2-17) (Andrews 1974) Unfortunately, there is no explicit solution for  ̂ ()

Simulations
We created simulations (Fig. 4) based on second-order stochastic stationary dynamics (equation 2-2) following from an underlying Gaussian Graphical Spectral (GGS) model.This hermitian precision matrix (Fig. 4 a) defining this GGS produces brain oscillations that resemble the so-called EEG Xi-Alpha model in the spectrum (Fig. 4 b) (Pascual-marqui, Valdes-sosa, and Alvarez-amador 1988).In other words, this autoregressive dynamic () is defined from (, ) (equations 2-4), produced by a hermitian precision tensor (connectivity)   () at all frequencies which also causes the oscillations in the alpha peak (Fig. 4 c).This hermitian precision tensor originates from spectral factors   () (Fig. 4 d) or the Hilbert transform of multiply lagged connectivity   () producing broad band dynamics (Fig. 4 e).This connectivity   () (Fig. 4 a) is employed as a baseline to compare the results and is more in correspondence with the actual target for state-of-the-art methods.The simulation pipeline that was repeated 100 times to produce the same amount of possible ground truth and observations is as follows.
1. Definition of random hermitian precision-matrix   ( 0 ) (Fig. 4a) (ground truth functional connectivity) given for a referential frequency component  0 meant to be the center of the alpha peak (Fig. 4b) ( 0 = 10Hz) in the cortical network subspace.The elements of this matrix must lie in the unitary circle of the complex plane for the reliability of binary classification measures used in this validation.2.
Construction of a synthetic spectral factorization for   ( 0 ) (Fig. 4c) or spectral directed model used to design a precision tensor at all frequencies ( ∈ ).This factorization was performed assuming that the reference spectral factors at ( 0 = 10Hz) are derived from a hermitian eigendecomposition   ( 0 ) =  † (equation 2-38), where   ( 0 ) is the directed transfer function at  0 extracted from the  † elements.These lags are associated with the spectral   ( 0 ) at  0 given the phase relation (equation 2-40), and the directed transfer function tensor   () at all frequencies ( ∈ ) (Fig. 4d) as in (equation 2-3) is then obtained from the definition for   () (equation 2-39) modulated by a Gaussian spectral form with center in the alpha peak The precision tensor   () (Fig. 4c) is then recomposed from the spectral factors (equation 2-41) assuming directed alpha transfer function tensor   () (equation 2-40) and Xi innovations precision tensor   () as in (equation 2-4).
The Xi factor was construed by following steps similar to alpha but modulated by a right-sided Gaussian spectral form with central frequency ( 0 = 0Hz).The precision matrix   was the inverse of the surface Laplacian in the whole cortical space, and the lags for the precision tensor construction were set to zero, corresponding to the properties of this process.
The precision tensor slicewise inverted   () =   −1 () (equation 2-41) (this is the ground truth covariance or crossspectrum   ()) is used to obtain complex-valued vectors of the Fourier transform (equation 2-42).Equivalent in this case to the Hilbert transform (, ) ( ∈ ) by means of a hermitian Gaussian random generator with sample number  = 600 at every frequency ( ∈ ).Brain oscillation time series () (Fig. 4e) are obtained by means of the inverse Fourier transform (ℱ −1 ) across frequencies ( ∈ ) of these samples (, ) to obtain the same number of time instances The Xi-Alpha process (, ) in the cortical network subspace (equation 2-42) was projected to obtain observations (Fig. 4g) at the sensor space () (equation 2-43), time-domain equivalent of (equation 2-9) with corresponding forwardoperator   (Fig. 4f), adding the Xi process and a white noise process in the whole cortical space.A white noise process () in the sensor space  was also added to the projected source space process.The composition of the confounding processes (Xi process defined by   and source and white sensor noises) with the alpha process defined by   ( 0 ) was adjusted to keep spectral energy in the alpha band (8-12 Hz) of 10% of the alpha process energy.() =   () + () (II-43) Two different simulations based on the Electromagnetic Forward Model (EFM) with EFM forward-operator   (equation 2-43) were defined: For an idealized "planar EFM" head (Fig. 4f1) and realistic "human EFM" head (Fig. 4f2) to produce the Xi-Alpha observations (Fig. 4g).The planar EFM (Fig. 4f1) was computed based on the bidimensional geometry of two concentric circles defining a planar cortex and scalp and a layout of equidistant planar scalp sensors (30 sensors).The human EFM (Fig. 4f2) was computed for the cortical surface extracted from a healthy subject T1 image, with the Boundary Element Method (BEM) implemented in SPM and a layout of the extended 10-20 system (30 sensors).A cortical network was defined on a subspace of 22 points for both the planar and human EFMs.For the planar cortex, these were equidistantly located, and for the human cortex, they were randomly located across different areas of the left (L) and right (R) hemispheres: Occipital (LO and RO), Temporal (LT and RT), Parietal (LP and RP) and Frontal (LF and RF).

Experimental confirmation
For the confirmation of HIGGS connectivity, we compared EEG against a more direct technique: electrocorticography (ECoG) (Fig. 5) leveraging the unique experimental setup that offers the advantage of large brain coverage ECoG on a healthy macaque (Nagasaka, Shimoda, and Fujii 2011).This comparison scenario is particularly interpretable in terms of the effect of volume conduction heterogeneities since the EEG is hidden under several tissue layers (Fig. 5 e2), and the ECoG is instead hidden under only one (Fig. 5 e1).See the macaque preparation for the surgical implantation of this ECoG in (Fig. 5 a) that, which is shown in the postsurgical X-ray (Fig. 5 b).The macaque sensor layouts for the simultaneous ECoG/EEG recordings consisted of 128 ECoG sensors placed upon the cortical surface in the left hemisphere and a low-density array of 20 EEG scalp sensors.The relative distribution of the macaque cortical surface segmented from the MRI is shown in (Fig. 5 c).The observations for EEG   () and ECoG   () were recorded simultaneously in the resting state.During the experimental session, the monkey was awake, blindfolded, and constrained to a sitting position.Both EEG and ECoG were synchronized to the trigger signal and downsampled to 1000 Hz, keeping 2 minutes of recordings.The artifact removal procedure included linear detrending with the L1TF package, average DC subtraction, and 50 Hz notch filtering.The spectral analysis (Fig. 5 d) of both ECoG (Fig. 5 d1) and EEG (Fig. 5 d2) signals reveals a larger power spectral density within the band (8-14 Hz) associated with the macaque alpha oscillations.The forward operators were obtained from a head conductivity model (Fig. 5 e) for the ECoG    (Fig. 5 e1) and EEG    (Fig. 5 e2) through FEM computations in SimBio using the macaque individual T1 MRI segmentation.The head model included five conductivity compartments: cortex, ECoG silicon layer, inner skull, outer skull, and scalp.We work at a subspace of cortical sources detected from ECoG, and in a similar fashion to the simulation study (Fig. 5 f), extracting the connectivity (Fig. 5 g) with all methods from both ECoG (Fig. 5 g1) and EEG (Fig. 5 g2) observations.We screen out the significant cortical network nodes (Fig. 5 f) from the ECoG by means of implementation for spectral analysis of the Structured Sparse Bayesian Learning (SSBL) (Paz-Linares et al. 2017).This implementation uses assumptions similar to those of HIGGS but with a diagonal covariance structure.We offer technical details in (Gonzalez-Moreira et al. 2020).The detection (Fig. 5 f) corresponding to the alpha band oscillations (Fig. 5 d1) yielded a large-scale network of distributed nodes, which strongly correlates with the findings for these ECoG cortical patterns in the left hemisphere (Q.Wang et al. 2019).The more extended Occipital (O) oscillations were accompanied by secondary cortical oscillations extended over the Temporal (T), Parietal (P), and Frontal (F) areas.This large-scale analysis differs from previously reported studies, which were limited to P <-> F interactions at ECoG sensor space as in (Papadopoulou, Friston, and Marinazzo 2019).We consider eight regions of interest (ROI) extracted manually  10Hz) by means of the eigen-decomposition.The process crossspectrum tensor is obtained by slice-inverting the precision tensor producing e. brain-oscillations by means of a hermitian random Gaussian generator.f.Two types of forward models, f1.planar and f2.human project these oscillations to the sensor space producing g. the Xi Alpha observations.from the left hemisphere cortical segmentation.These are Inferior Occipital (IO), Superior Occipital (SO), Posterior Temporal (PT), Anterior Temporal (AT), Inferior Parietal (IP), Superior Parietal (SP), Inferior Frontal (IF), and Superior Frontal (SF).

Unbiased functional connectivity via hgLASSO inverse solution of the GGS model
As an essential step (Fig. 6), we verify the properties of the MAP1 estimator for the precision matrix   () representing functional connectivity as described by a Gaussian Graphical Spectral (GGS) model of oscillatory brain networks at a generic frequency .See Materials and Methods section 2.1 for the interpretation of such a functional network model and MAP1 inverse solution  ̂ () for the functional connectivity.We remind the reader that the algorithm revindicated for this MAP1 inverse solution places the hermitian graphical LASSO (hgLASSO) a priori model upon the hermitian graph elements in   (), which explains brain-oscillations (, ) in the GGS model.In Materials and Methods section 2.3, we describe the implementation of the hgLASSO algorithm that we demonstrate here for the first time fulfilling stable and scalable unbiasedness conditions.Here, amplitudes |  ()| of the hermitian graph elements are employed as the proxy for functional connectivity and to illustrate the algorithmic properties.The hgLASSO a priori model is incorporated into both multistep MAP1 and onestep MAP2 inverse solutions, as described in Materials and Methods section 2.2 for the Hidden GGS (HIGGS) model explaining oscillations (, ).We illustrate the properties of our algorithm with an example at different dimensionalities, 30, 100, and 1000, yielding 900, 10,000, and 1000,000 functional connectivities to be estimated, respectively.Our validation strategy (Fig. 6) is to compare the empirical distribution of the empirical distribution of the unbiased statistic for the estimated precision matrix  ̂ () against the Rayleigh distribution for the theoretical   ().A Rayleigh distribution is derived from the complex-valued extension of (Janková and van de Geer 2018) for the GGS model described in Materials and Methods section 2.3.First, we employ a blockwise chained (blocks were mutually dependent) simulated precision-matrix structure   () as in (Tugnait 2019), illustrated here with hot colormaps of a typical trial (Fig. 6 a).These matrices were integrated into a hermitian Gaussian generator generating the sampled covariance matrix from sampled instances (a sample number of 100 times the matrix dimensions).Second,  ̂ (), a precision matrix, is estimated with our hgLASSO algorithm (Fig. 6 b), and third, from this estimate, the unbiased precision matrix  ( ̂ ()) (Fig. 6 c) is computed as proposed by the theory.Fourth, we obtain the histogram of a z-statistic () (Fig. 6 d) for this unbiased estimate at the null hypothesis matrix entries (zero precision-matrix entries   (; , ) = 0), and this histogram accurately resembles the theoretical Rayleigh distribution.Computing the z-statistic was performed by scaling the unbiased precision matrix using the theoretical variances  ̂(; , ).Fifth, we remove the values under the Rayleigh threshold (Rayleigh corrected precision-matrix) obtained from the 0.05 p value of the theoretical Rayleigh distribution.With the Rayleigh threshold, from the originally dense unbiased precision matrix (Fig. 6 c), we obtain an improved corrected result (Fig. 6 e) in comparison to the sparse biased hgLASSO precision matrix (Fig. 6 b).The hgLASSO (base of all estimated maps), as we illustrate in the ultrahigh dimensionality, a robust convergence pattern (Fig. 6 f) for multiple repetitions of the experiment (100 trials).The z-statistic histograms for these 100 trials reflect the high coincidence with the theoretical Rayleigh distribution (Fig. 6 g).

Erroneous multistep and exact onestep functional connectivity for the HIGGS model in simulations
The proof of concept for the HIGGS model (Fig. 7) employs random precision-matrices instances   () similar to the previous results in section 3.1 but yielding binary functional connectivity as defined with the amplitude of hermitan graph-elements |  ()|.Our binary "ground-truth" functional connectivity |  ()| (Fig. 7 left-center) is a much more plausible target of estimation errors via classification scores (Fig. 7 right).Simulations of the HIGGS model producing EEG observations (, ) are due to brain oscillations (, ) steaming from an alpha rhythm peaking at 10 Herts, which follows from an interpretation of the GGS model in terms of second-order stochastic stationary dynamics.See the description of these simulations in Materials and Methods section 2.4.The brain oscillations (, ) are projected through forward operators   that follow two types of EFM geometrical designs: 1) Planar EFM (Fig. 7 left-top), recreating an ideal circumstance where concentric circles recreate the cortical network, head and EEG sensor geometries, and EEG fields lead through homogeneous conductance.2) Human EFM (Fig. 7 left-bottom), recreating an average human EEG circumstance.Then, the GGS model precision matrix   () (Fig. 7 left-center), which describes brain oscillations (, ) at the alpha peak, is the target of the estimation  ̂ () with all types of inverse solutions exposed here (Fig. 7 a1-g1 and Fig. 7 a2-g2).See these inverse solutions in Materials and Methods section 2.2.First, three onestep MAP2 inverse solutions are implemented as successive approximations due to MAP1 inverse solutions within an Expectation-Maximization (EM) loop for 1) brain oscillations (, ) in the EM expectation stage and 2) functional connectivity  ̂ () in the EM maximization stage.The functional connectivity MAP1 (maximization stage) determines a GGS model precision matrix  ̂ () from the expected covariance matrix  ̂ () via the unbiased Hermitian Graphical LASSO (hgLASSO) algorithm (Fig. 7 a1 and Fig. 7 a2).See this algorithm in Materials and Methods sections 2.1 and 2.3.Here, we also introduce violations to the hgLASSO unbiasedness conditions in the functional connectivity MAP1 via the hermitian graphical Ridge (hgRidge) (Fig. 7 b1 and Fig. 7 b2) and hermitian graphical Naïve (hgNaïve) (Fig. 7 c1 and Fig. 7 c2).Second, the multistep MAP1 inverse solutions are implemented for 1) brain-oscillations (, ) via ELORETA (Fig. 7d1 and Fig. 7d2) (Pascual-Marqui et al. 2006) and LCMV (Fig. 7 e1 and Fig. 7 e2) ( Van Veen et al. 1997) in the first step and 2) functional connectivity  ̂ () in the second step and from the sampled covariance-matrix  ̂ () via the unbiased Hermitian Graphical LASSO (hgLASSO) algorithm.Here, we also introduce a violation to the hgLASSO hermitian graph elements in the functional connectivity MAP1 via the approximation of Graphical LASSO (gLASSO) (J.Friedman, Hastie, and Tibshirani 2008;G. L. Colclough et al. 2015).This violation is incorporated in the second step for either multistep MAP1 inverse solutions with the first step via ELORETA (Fig. 7f1 and Fig. 7f2) or LCMV (Fig. 7g1 and Fig. 7g2).The sampled covariance matrix  ̂ () is determined from the first-step brain oscillations (, ) defined as a filtered process (G.L. Colclough et al. 2015) and not its Hilbert transform (Wodeyar 2019).To achieve a fair comparison, we employed the optimal regularization parameter for the initial inverse solutions eLORETA and LCMV selected by generalized crossvalidation of the sampled covariance matrix for the data  ̂ ().The accuracy of the functional connectivity estimates was assessed by the classification scores of a binary ground true |  ()| (Fig. 7 left-center) by measures of the Receiver Operator Characteristic (ROC) (Fig. 7 right).These measures are given for inverse solutions | ̂ ()| due to the planar (Fig. 7right-top) and human (Fig. 7right-bottom) EFMs.The measures used in the radar graphs were the global AUC (area under the curve) and partial SENS (sensitivity), SPEC (specificity), PREC (precision), and RECALL (F1 or Fisher scores) measured at the optimal operating point of the ROC curve.

Confirmation of the performance of the functional connectivity for the HIGGS model in macaque EEG/ECoG experiments
We leverage an experimental setup (Fig. 8) offering for macaque 1) large-scale coverage resting-state ECoG observations   (, ) onto the left cortical hemisphere and 2) simultaneous EEG observations   (, ) (Nagasaka, Shimoda, and Fujii 2011).The ECoG observations   (, ) are from a high-density subdural sensor array that was placed surgically onto the macaque cortex and provides (Fig. 8 left-top).See the experimental preparation and inverse solutions from ECoG  ̂  () and EEG  ̂  () in Materials and Methods section 2.4, in this case targeting a functional network underlying the resting-state alpha rhythm.

Theory
The state-of-the-art connectivity inverse solution MEG-ROI-nets Graphical LASSO (gLASSO) (G.L. Colclough et al. 2015), which is employed here as a baseline here for comparison, was completely limited to only real-valued operations.The flaw is twofold due to instability in large scale and the lack of MEG-ROI-nets gLASSO algorithms capable of dealing with models that dwell in the complexvariable space of the Fourier or Hilbert transform under the GGS assumptions (K.J. Friston et al. 2012;Bien and Tibshirani 2011; C.-J. Hsieh et al. 2011;C. Hsieh 2014;Tugnait 2019;Schreier and Scharf 2010;Andersen et al. 1995;Bollhöfer et al. 2019;C.-J. Hsieh et al. 2013C.-J. Hsieh et al. , 2012)).It was therefore our task to provide with a solution to estimation errors produced by such algorithms as discussed in Introduction section 1.3.With same multistep methods the estimation errors of the state-of-the-art gLASSO approximating the hermitian graph-elements in this HIGGS (f1,g1;f2,g2) fall beyond even in ideal circumstances.In 100 trials of a binary connectivity ROC measures show these estimation errors and improvement with our direct solution with hgLASSO (a1,a2) which may not be possible with hgRidge or hgNaïve estimation (b1,c1;b2,c2).
A gLASSO inverse solution based upon a real-valued sampled covariance  ̂ () for (, ) defined as the filtered process for () or sometimes the Hilbert envelope −|(, )| is claimed to encode in  ̂ () all connectivity information in state-of-the-art neuroimaging practice but only does it partially.Here, we employed the MEG-ROI-nets gLASSO (G.L. Colclough et al. 2015) upon the real-valued sampled covariance matrix  ̂ () for the filtered process (, ) which does not uphold the Gaussian Graphical Spectral (GGS) model assumptions but is a less gross approximation than Hilbert envelope −|(, )| (Matthew J. Brookes, Hale, et al. 2011;Matthew J. Brookes, Woolrich, et al. 2011b;Giles L Colclough et al. 2015;Vidaurre, Abeysuriya, et al. 2018;Vidaurre, Hunt, et al. 2018;Tewarie et al. 2019;Hillebrand et al. 2012;Matthew J. Brookes et al. 2007;Larson-Prior et al. 2013;Hipp et al. 2012).For oscillatory network phenomena estimation of connectivity  ̂ () required both amplitude and phase information that was determined with the Hermitian Graphical LASSO (hgLASSO) (Nolte et al. 2020).We remind the reader a hgLASSO inverse solution introduced in Materials and Methods section 2.1 and section 2.3 is only resource to find estimates  ̂ () for the complex-valued precision matrix of such GGS model.At every frequency, the complex-valued sampled covariance matrix  ̂ () for (, ) defined as the Hilbert transform for () conforms to the so-called cross-spectrum.Only this quantity, or equivalently its matrix inverse or precision-matrix  ̂ (), fully encodes the dynamical properties of oscillatory components (, ) in activity ().Our hgLASSO algorithm showed consistency in simulations to the theoretical Rayleigh tendency of amplitudes for the hermitian graphelements in complex-valued Gaussian Graphical Spectral (GGS) models (Fig. 6).Such Rayleigh tendency in hermitian graph-elements (Andersen et al. 1995;Tugnait 2019; Schreier and Scharf 2010) was a natural extension of the Chi-Square tendency in the real-valued case (J.Friedman, Hastie, and Tibshirani 2007; Janková and van de Geer 2018).The complexity of the proposed algorithm was bounded by the hermitian matrix eigen decomposition and we shall note no gLASSO algorithm was available which fulfills statistical guarantees, performs stably in the large dimensionality, and deals with a complex-valued cross-spectrum.Our hgLASSO is therefore an optimal procedure amongst state-of the-art algorithms for this type of graphical model either in real or complex-variable (W.subspace of prescreened network nodes ruled out other sources of distortions and allowed us to evaluate in isolation the estimation errors exposed in Introduction section 1.2.These estimation errors are caused in theory by misspecification of the HIGGS model in inverse solution targeting connectivity estimates  ̂ ().Other sources of distortion were the localization errors and leakage in activity estimates (, ) that would be produced by multistep (in the first step) or onestep (successive approximations) upon the whole cortical space (G.L. Colclough et al. 2015;Palva et al. 2018; M. J. Brookes, Woolrich, and Barnes 2012;Wens et al. 2015).Inverse solutions for HIGGS connectivity  ̂ () were implemented via the hgLASSO, for both onestep, in successive approximations in an iterative fashion, and multistep, in a first step, estimating the cross-spectrum  ̂ () as described in Materials and Methods section 2.2.Our onestep connectivity inverse solution was an extension to the Expectation Maximization (EM) algorithm as with the identification of the covariance matrix  ̂ () in a Covariance Component Model (CCM) (C.Liu and Rubin 1994;Galka, Yamashita, and Ozaki 2004;K. Friston et al. 2006K. Friston et al. , 2007) ) or identification of the autoregressive-coefficients  ̂ () in a State Space Model (SSM) (Jorge J. Riera et al. 2004;Faes, Stramaglia, and Marinazzo, n.d.;Galka et al. 2004).Implementation of the hgLASSO algorithm defined the flavor of the HIGGS model onestep identification connectivity defined as the complex-valued precision matrix  ̂ ().HIGGS obtaining an unbiased, scalable and stable inverse solution for a locally optima GGS within EM iterations (McLachlan and Krishnan 2007) was not considered by previous CCM or SSM identification approaches.

Validation
We spoused our validation of estimation errors in HIGGS inverse solutions in electroencephalogram (EEG) in simulations (Fig. 7) and macaque real data (Fig. 8), since this is the "worst case scenario" compared to magnetoencephalogram (MEG) or electrocorticogram (ECoG).EEG is exquisitely sensitive to conductivity heterogeneity of head tissue layers, which produces quite large distortions of the electric lead fields   (Jörge Javier Riera and Fuentes 1998;Baillet, Mosher, and Leahy 2001;Grech et al. 2008).Additionally, the high conductivity of the scalp leads to blurring of the EEG potential, which aggravates the estimation errors of any activity estimates (, ) or its cross-spectrum  ̂ (), that are transferred to connectivity estimates  ̂ () (Van De Steen, Valdes-Sosa, and Marinazzo 2016).Thus, EEG results for connectivity  ̂ () presented here can be considered a lower bound for those of MEG or ECoG.The methods under evaluation fell into three different modalities: I) the onestep successive approximations, employing our hgLASSO unbiased estimator for  ̂ () (Fig. 7 a1, a1), and employing estimators hgRidge (Fig. 7 b1, b2) and hgNaïve (Fig. 7 c1, c2) which violated this unbiasedness condition in  ̂ ().II) the multistep methods with eLORETA (Pascual-Marqui et al. 2006) (Fig. 7 d1, d2) and LCMV (Van Veen et al. 1997) (Fig. 7 e1, e2), both employing the same hgLASSO unbiased estimator for  ̂ ().It has been argued (Wodeyar 2019) that these types of multistep estimators reduce connectivity estimation errors.III) the same multistep methods with eLORETA (Fig. 7 f1, f2) and LCMV (Fig. 7 g1, g2) but employing MEG-ROI-nets gLASSO estimator (J.Friedman, Hastie, and Tibshirani 2008) for the real-valued  ̂ () implemented in the MEG-ROI-nets package (G.L. Colclough et al. 2015) which does not uphold the GGS assumptions.Even with the relatively simple ground truth some methods produced large estimation errors in these simulations (Fig. 7) indistinctly, employing I) an idealized or planar (Fig. 7 left top) and II) a human (Fig. 7 left bottom) Electromagnetic Forward Model (EFM).In simulations (Fig. 7 right) partial classification measures RECALL, PREC, and SENS quantify and confirm the essential flaw with multistep methods employing II) hgLASSO or III) gLASSO.The perfect multistep estimation in the planar EFM (Fig. 7 a1), which benefited from the renowned statistical goodness of eLORETA and our unbiased hgLASSO algorithm, was undone in the human EFM (Fig. 7 a2).Another multistep solution based on LCMV (Fig. 7 e1, e2), with excellent performance anywhere else, does not reach expectations even in the planar EFM ideal circumstances (Fig. 7 e1).Only the onestep hgLASSO sparse inverse solution showed no major estimation errors (Fig. 7 a1, a2) and outscores any other non-sparse onestep (Fig. 7 b1,b2,c1,c2) or sparse multistep (Fig. 7 d1, d2, e1, e2) solutions in both planar (right-top) and human (right-bottom) models.Connectivity leakage measured as 20% in PREC or RECALL (shown as colormap false positives) is a systematic property of these multistep methods appearing even in a relatively simple human network model.The situation is even worse when using the gLASSO approximation for these methods, which corrects leakage in either simulation model (Fig. 7 f1, f2, g1, g2) as claimed before (Wodeyar 2019;G. L. Colclough et al. 2015) but produces extreme estimation errors, yielding random classification according to all classification measures.gLASSO fails under the specific assumption of a GGS model, which bases connectivity on hermitian graph elements.Previous simulation studies validating this correction method did not corroborate them under the GGS assumption and were limited to even simpler network models (up to 5 nodes) (G.L. Colclough et al. 2015).The qualitatively sparse pattern (Fig. 7) obtained in simulations via the onestep (Fig. 7 a1, a2) and multistep eLORETA (Fig. 7 d1, d2), with both methods employing the unbiased hgLASSO inverse solution, resembled the sparsity simulated in the ground truth (Fig. 7 left center) and minimized connectivity estimation errors.Sparsity renders factual biological circumstances such as efficiency in largescale network phenomena (Bullmore and Sporns 2009), as we illustrated with the macaque ECoG (Fig. 8) with the sparse alpha network (Fig. 8 left center) and sparse connectivities (Fig. 8 a1, b1) obtained from any (onestep or multistep) method employing the unbiased hgLASSO estimator.Remarkably, the dimensionality of the macaque ECoG alpha network (Fig. 8 left-center) ruled out the possibility of using current implementations of the gLASSO algorithm implemented in (J.Friedman, Hastie, and Tibshirani 2008;C. Hsieh 2014;C.-J. Hsieh et al. 2013), or by MEG-ROI-nets (G.L. Colclough et al. 2015).Confirmation in macaque shows the connectivity determined for this type of network from the ECoG (Fig. 8 left-top) and EEG (Fig. 8 left-bottom) by means of the onestep hgLASSO (Fig. 8 a1, a2) and multistep eLORETA hgLASSO (Fig. 8 b1, b2).The macaque experiment (Fig. 8) confirms the benefit in our onestep hgLASSO inverse solution, targeting a severely ill-posed inverse problem in recovering connectivities for a large-scale alpha network (Fig. 8 left-center) based only on EEG 19-sensor observations (Fig. 8 left-bottom).Indeed, for both the onestep and multistep methods, there are striking similarities in the ECoG (Fig. 8 a1, b2) and EEG (Fig. 8 a2, b2) solutions.Our onestep hgLASSO solutions (Fig. 8 a1, a2) maximize the sparsity.With the Kullback-Leibler divergence (KLD) kernel (Fig. 8 a3, b3), we illustrate contraposition between these sparsity and distortion levels, measured as multivariate relative entropy for GGS models.This measure represents divergence (incongruence) from the identity matrix for the product between the EEG precision matrix and ECoG covariance that would express a multivariate effect of distortions locally (for each graph element).The norm for this matrix, the KLD metric and others (Reimanian and Log-Euclid) shows a 1/3 congruence improvement.For a vector  or matrix  (as defined in I-1) the -th, or , -th element.
(I-10) (), or (), or () Types of probability density functions for (as defined in I-1) a random scalar , or random vector , or random matrix .
(I-13) (|  ,   ), or (|  ,   −1 ) Real-valued Gaussian distribution for a vector  conditioned to (as defined in I-12) the mean  and real valued symmetric ensemble covariance-matrix   or precision-matrix   (as defined in I-12).
(I-15) (|) Exponential distribution for the scalar  conditioned to (as defined in I-12) the scalar parameter of shape .
(I-16) (|, ) Gamma distribution for the scalar  conditioned (as defined in I-12) to the scalar parameters of shape  and rate .
Sum operator along index , valid for , , ,  in matrices or vectors (as defined in 4.4) that can be defined within a countable set  (as defined in 4-1) or from 1 up to a maximum .

𝑡=1
Product operator along index , valid for , , ,  in matrices or vectors (as defined in 4.4) that can be defined within a countable set  (as defined in 4-1) or from 1 up to a maximum .(II-3) , with T = || and  ∈  Domain of time for processes, local neural currents or MEG/EEG/ECoG observations and their perturbations,  with size T, instances in this space are defined as  ∈ .
(II-4) , with F = || and  ∈  Domain of frequency for spectral processes, spectral local neural currents or spectral MEG/EEG/ECoG observations and their perturbations,  with size F, instances in this space are defined as  ∈ .
(II-7) (), and (, ) Local neural currents or cortical activity (), in time  (as defined in II-3), and its Hilbert transform, cortical spectral activity or cortical oscillations (, ), in time  (as defined in II-3) and frequency  (as defined in II-4).Either vectors are defined in the space  (as defined in II-1).
(II-14)  ̅  (), and  ̅  () First type maximum a posteriori estimator (bar  ̅ ) of a covariance-matrix  (as defined in 5-12), and which if dependent of the second type maximum a posteriori estimators.
(II-15)  ̅  (), and  ̅  () First type maximum a posteriori estimator (bar  ̅ ) of a covariance-matrix  (as defined in 5-12), and which if dependent of the second type maximum a posteriori estimators.

Proof of Lemma 1
The general strategy the identification of the HIGGS model (equation III-1) is based on the EM algorithm: a second-type maximum a posteriori estimation, performed iteratively for () (parameters) and  (hyperparameters) (Dempster, Laird, and Rubin 1977;C. Liu and Rubin 1994;McLachlan and Krishnan 2007;Wills, Ninness, and Gibson 2009) Analyzing the exponent (equation III-7) from the Bayesian expression (equation III-6) is enough to find the structure of the posterior distribution (()|(),  () ).
̂(+1) =   ( () (| ̂ ))  ()  The final expression of the target functions above for the sources ℒ () (  ) (equation IV-5), and residuals ℒ () (  ), (equation IV-6), defining the solutions of iterated Gaussian Graphical Spectral (GGS) models, is due to the priors upon the precision-matrices   (equation IV-7) and   (equation IV-8).(  ) = (Π(  ⊙   )|  ) Where the prior structure of the GGS precision-matrices is represented wuith the ad-hoc   for the sources (equation IV-6), and   for the residuals (equation IV-7).The scale (regularization) parameters control the level of penalization exerted by these priors, given by   that impose Π the L1 or L2 norm upon the multivariate   , and   that may be interpreted as the residuals inferior limit represented by the scalar precison   2 due to a reparameterization of the residuals precision-matrix   =   2   .The estimators for the second-type maximum a posteriori are then expressed as (equation IV-9), for  ̂ (+1)   ̂ The hierarchical representation of the Gibbs prior with hermitian graphical LASSO exponent Π( ⊙ ) = ‖ ⊙ ‖ 1 , can be built on corollaries of the Andrews and Mallows Lemma (Andrews 1974), for the extension of real-valued Laplace probability density function to the real-valued or complex-valued matrix case, by considering improper density functions or simply more general measurable spaces.By the Andrews and Mallows Lemma in the real-valued Laplace, also for the real-valued or complex-valued case the integral representation holds (equation V-3).
(, ) = ∏  1  (|(, )||0,  2 (, )  ⁄ )  ( 2 (, )|1,   2  2 (, ) 2 ⁄ ) q , ∎ (V-8) Remark: In Lemma 2 we stablish a statistical equivalence between the Gibbs probability density function with argument in the hermitian graphical LASSO prior and a hierarchical improper representation through a Gaussian probability density function of the conditional expectation |, and variances (weights)  with Gamma probability density function.Remarkably, this representation we are not using a probability density function in the strict mathematical sense but a measure density of a product of measurable spaces for | and .

VI. Concavity of the first-type maximum a posteriori with hgLASSO LQA prior (Lemma 3) Lemma 3 (concavity of the Local Quadratic Approximation for the hermitian graphical LASSO prior):
The target function ℒ(, ) (− of the improper probability function for the measurable product space | × ) may be represented as (equation VI-1) an it is strictly concave on the intercept of the region of positive definiteness of its arguments (equation VI-2) and the region comprehended by the set of inequalities (equation VI-3).
Then, ℒ(, ) has a minimum within these regions defined in (equation VI-2) and (equation VI-3) given by the intercept of the system (equation VI-4) and (equation VI-5).

Proof of Lemma 3
With the hierarchical representation of the hermitian graphical LASSO precision-matrix  prior, extension of Andrews and Mallows Lemma (Lemma 2), we attain a modified target function ℒ(, ) of the precision-matrix that builds in the combination of the sampled covariance-matrix (|) Wishart likelihood dependent on the local quadratic approximation (improper probability density function of | × ) (equation VI-1).Other terms are related to the normalization constant of the precision-matrix | Gaussian prior and the precision-matrix variances (weights)  Gamma prior.To build the target function ℒ(, ) we may employ the Bayes theorem for the improper posterior probability density function (, |) based on the Wishart likelihood (|) and improper probability density function for (| × ) = (|)() (equation VI-6).
1 2 (VI-30) Given the positive definiteness of  and , it can also be deduced that the solution  (equation VI-30) is positive definite by following the steps in proposition 3 of Lemma 4. ∎ According to this local quadratic approximation, we can define a new random matrix through the scaling transformation  ̃=  ⊘  (standard precision-matrix), so that its prior is a Gibbs probability density function of the squared L2 norm (equation VI-31).
( ̃) ∝  −  2 ‖ ̃‖2 2 (VI-31) For the solution based of this standard precision-matrix we also standardize the sampled covariance-matrix due to the transformation  ̃= ( −1 ⊘ ) −1 that keeps the stochastic properties of the Wishart distribution when conditioned to  ̃.
Positive definiteness and hermiticity of its solution  ̃, as a natural property of this estimator, is required.Also, it can be checked that the solution  ̃ (equation VIII-12) must share the same eigenspace with  ̃, since from (equation VIII-11) two auxiliary second order matrix equations hold, i.e.  ̃2 +  ̃ ̃−  q =  and  ̃2 +  ̃ ̃−  q =  q , given the left and right multiplication by the Standard Precision matrix.This imply that any solution  ̃ commute with  ̃, and thus they share the eigenspace.∎ Let's show first that the proposed solution is positive definite.Given the positive definiteness and hermiticity of the complex-valued matrix  ̃ it admits a singular value decomposition  ̃=  ̃ ̃ ̃ † with real and positive singular values  ̃.Thus, the argument in the matrix square root of formula (equation VIII-12) admits a singular value decomposition of the kind  ̃2 + 4 q =  ̃( ̃ + 4 q ) ̃ †, where its singular values, given by  ̃ + 4 q , are also real and positive.In consequence the square root term (equation VIII-12) has also real and positive singular values given by the following singular value decomposition: √ ̃2 + 4 q =  ̃√ ̃ + 4 q  ̃ † (VIII-13) Finally, the singular value decomposition of the standard precision matrix estimator  ̃ (equation VIII-12) can be expressed as (equation VIII-14).The product of diagonal matrices always commutes, so we can directly check that (equation VIII-15) is also the singular value decomposition of  ̃ ̃.To check that the proposed estimator (equation VIII-12) satisfies the equation (equation VIII-11) it is enough to check that it also satisfies the pair of equations  ̃ ̃=  q −  ̃2 and  ̃ ̃=  q −  ̃2.The left side in both equations are equal, i.e.

Fig
Fig. 1| Illustration of the ontological levels involved and basic multistep identification of connectivity.a. Generation by the electromagnetic forward model of the periodic components or oscillations in MEG/EEG observations (, ) from brain oscillations (, ) emerging from functional graph-elements or functional connectivity   () defining a cortical oscillatory network at a given frequency .b.Functional connectivity distortions due to multistep inversesolutions of the optimal cortical oscillations (, ) explaining the observations and functional connectivity  ̂ (), from pseudodata (, ).Even with with perfect spatial localization estimates (, ) by any first step inverse-solution produce however false positives and false negatives functional graphelements in the second-step inverse-solution  ̂ ().

Fig. 2|
Fig. 2| Connectivity distortions in the real-valued approximation of a hermitian GGS model.a. Undirected network as defined by the amplitude of hermitian graph-elements or entries of the hermitian GGS precision-matrix.b.Binary GGS-precision-matrix amplitudes   () are perfectly retrievd by estimation  ̂ () based on a hermitian GGS model with hermitian graphical LASSO prior and distorted by estimation based on a real-valued GGS model with graphical LASSO prior.

Fig
Fig. 3| Theoretically correct procedure for the estimation of connectivity and to eliminate distortions.a.The ideal but impractical solution is to bypass intermediate estimators (, ) onestep search of optimal connectivity  ̂ () explaining the observations.b.Our implementation of efficient sequential approximations to the onestep procedure by means of Expectation Maximization iterations of ̂( ) (, ) and  ̂ () () produces perfect functional graph- elements due to statistical guaranties.

Fig
Fig.4|Specification of a linear autoregressive dynamics in the frequency domain based on the Xi and Alpha model a. Alpha connectivity is the ground true to be mixed later with that of the Xi process.The starting point are these Alpha and Xi precision-matrices at 10Hz for an oscillatory network process represented in a binary map that corresponds to edges in such network graph-elements.b.This is extended for all frequencies via an amplitude-phase transformation to recreate c. the precision-tensor employing d. directed transfer function tensor for the Alpha process by means of its composition in spectral factors with the Xi innovations The factors are derived from the central slice (10Hz) by means of the eigen-decomposition.The process crossspectrum tensor is obtained by slice-inverting the precision tensor producing e. brain-oscillations by means of a hermitian random Gaussian generator.f.Two types of forward models, f1.planar and f2.human project these oscillations to the sensor space producing g. the Xi Alpha observations.

Fig
Fig. 5| Confirmation of HIGGS connectivity based in EEG that is recorded simultaneously with ECoG implanted in the macaque.a. Surgical preparation for implantation onto the macaque cortex macaque of a high-density ECoG array embedded in a silicon layer.b.A post-surgical X-ray image showing this implantation.c.Digital preparation based on the macaque MRI segmentation of the cortex, inner skull, outer skull and scalp, and conductance model for this segmentation based on registration with the MRI of the high-density ECoG implantation (128 blue sensors) and low-density EEG layout (20 green sensors).d.Power spectral density of the simultaneous recordings for ECoG (d1) and EEG (d2) highlighting the alpha band within 8-14Hz employed later as data for the identification of network connectivity.e. Electromagnetic Forward Model (EFM) employed to compute the connectivity inverse-solutions from the ECoG (e1), with Lead Field dependent on two conductance layers (cortex and silicon), and EEG (e2), with Lead Field dependent on four conductance layers (cortex, inner skull, outer skull and scalp).ECoG recordings and their Lead Fields provide a more fine-grained reference for confirming connectivity estimators and measures of distortions for the EEG f.Significant cortical subspace based on the ECoG alpha band data to facilitate the comparison in terms of connectivity due to low-density of the EEG.The detected subspace is revealing a large-scale alpha cortical network that extended over the Inferior Occipital (IO), Superior Occipital (SO), Posterior Temporal (PT), Anterior Temporal (AT), Inferior Parietal (IP), Superior Parietal (SP), Inferior Frontal (IF) and Superior Frontal (SF) areas.g) For the subspace covered by this network the HIGGS connectivity inverse-solutions were computed from the ECoG (g1) and EEG (g2).

Fig
Fig. 6| Experiment to evaluate unbiasedness conditions of the Hermitian Graphical LASSO (hgLASSO) algorithm in different scales at the top: low (first row), high (second row) and ultrahigh (third row).The hgLASSO algorithm performs the MAP1 for the presicion-matrix as described in Materials and Methods sections 2.1 and 2.3.From left to right in hot colormap a. the typical simulated precision-matrix, b. the hgLASSO estimation, c. the debiased estimator, with d. the empirical distribution, compared of the theoretical Rayleigh distribution, and e. the corrections due to thresholds of this distribution.At the bottom for the ultra-high dimensionality the robust convergence pattern in 30 iterations of the hgLASSO likelihood in 100 trials (f) and the seemingly perfect Rayleigh statistics in these 100 trials (g).
On the one hand, a fine-grained EFM forward-operator    for the ECoG observations provides accurate functional network detection (Fig.8left-center), as well as functional connectivity by any onestep (Fig. 8a1) or multistep (Fig. 8 b1) inverse solutions  ̂  ().On the other hand, sensor low density and several tissue layers yield a coarse-grained EFM forward operator    for the EEG observations (Fig. 8 left-bottom) and therefore as expected incongruent inverse solutions  ̂  () (Fig. 8 a2 and Fig. 8 b2) regarding an ECoG reference  ̂  () (Fig. 8 a1 and Fig. 8 b1) as proposed in (Q.Wang et al. 2019).We employ the Kullback-Leibler Divergence (KLD) (Kullback 1968) as the proxy to measure this incongruence.The KLD distance measures is the relative entropy between pairs of multivariate probability distributions.In this case, the distributions are Gaussian Graphical Spectral (GGS) models  ((, )| ̂  ()) and  ((, )| ̂  ()) expressed as a function of the estimated precision matrices  ̂  () and  ̂  () summarizing the GGS properties.A multivariate effect of distortions locally (for each graph element) for EEG  ̂  () and regarding ECoG  ̂  () may be measured by means of the KLD kernel for onestep (Fig. 8 a3) and multistep (Fig. 8 b3) inverse solutions.For a quantitative analysis of the performance (Fig. 8 right), we report distances based on the KLD metric and others (Riemannian and Log-Euclid).

Fig
Fig. 7| Simulated experiment of a relatively simple oscillatory network HIGGS model (left-center) producing observations through 1. an idealized "Planar" (left-top) or 2. average "Human" (left-bottom) EFM.Estimation errors regarding a ground-true hermitian precision-matrix with binary amplitudes (leftcenter) may be judged in the colormaps of the estimation under the Planar (a1-g1) and Human (a2-g2) EFM.These estimation are characteristic of multistep methods (based on first-step ELORETA and LCMV) in an average EEG situation (d2,e2) even employing the hgLASSO unbiased estimator.With same multistep methods the estimation errors of the state-of-the-art gLASSO approximating the hermitian graph-elements in this HIGGS (f1,g1;f2,g2) fall beyond even in ideal circumstances.In 100 trials of a binary connectivity ROC measures show these estimation errors and improvement with our direct solution with hgLASSO (a1,a2) which may not be possible with hgRidge or hgNaïve estimation (b1,c1;b2,c2).

Fig
Fig. 8| Experimental confirmation of performance in macaque simultaneous EEG/ECoG recordings (left).A large-scale alpha oscillatory network (leftcenter) screnned from ECoG with certain method is the target of HIGGS inverse-solutions.Incongrunce base on comparing these solutions under the ECOG (left-top) and EEG (left-bottom) EFM.Employing as reference the fine-grained ECOG solutions (a1,b1), incongrunce of the EEG solutions are measured with Kulback-Leibler, LogEuclidean and Remaniean metrics (right) confirming the performance in Fig. 2.This incongruence due to leakage and localization errors is represented with the Kulback-Leibler Divergence (KLD) kernel "EEG/ECoG + EEG\ECoG -2I" (a3,b3), for EEG and ECoG precision-matrices.
The scaled Gaussian mixture (equations 2-17) is a type of representation of the hgLASSO prior (equations 2-3) by a Gaussian distribution of   () conditioned to other hyperparameters (variances ()) with a gamma distribution.As a result, we obtain a Local Quadratic Approximation (LQA) of the target function used in (equations 2-6 or 2-11) ℒ((), ()) in terms of the weighted hermitian graphical Ridge (hgRidge) due to the modified LASSO prior (equations 2-18).The estimation derived from the LQA poses a unique solution, given by Lemma 3 in VI.Concavity of the first-type maximum a posteriori with hgLASSO LQA prior (Lemma 3).

Supplementary Information I. Notation for variables and mathematical operations
This work was supported in part by the National Nature and Science Foundation of China (NSFC) under founding No. 61871105 and the CNS program of UESTC under founding No. Y0301902610100201.
(I-1) , or , or , or  Respectively denote a scalar  (italic lowercase), or vector  (bold italic lowercase), or matrix  (bold capital), or set  (double struck capital).(I-2) (), or (, ) (I-21) ‖‖  ,  = 1,2, or ‖‖  ,  = 1,2 L1 or L2 norms of the vector  or matrix .Extreme values of the scalar function , correspondingly minimum or maximum, in the matrix argument  , valid for a scalar  or a vector .Zeros of the scalar function  in the matrix argument , valid for a scalar  or a vector .
(II-1) , with E = || and  ∈  Physical space of MEG/EEG/ECoG observations , with size E, sensors in this space are defined as  ∈ .(II-2) , with G = || and  ∈  Physical space of local neural currents or cortical generators of the MEG/EEG/ECoG , with size G, cortical generators in this space are defined as  ∈ .

valued Andrews and Mallows Lemma: Local Quadratic Approximation (LQA) of the hermitian graphical LASSO (hgLASSO) prior (Lemma 2) Lemma 2 (corollary of Andrews and Mallows for the hermitian graphical LASSO):
The second-type maximum a posteriori estimates for the scalar precision  ̂ 2 (+1) , a reparameterization   =   2   for the hermitian graphical model of the residuals, can be obtained in terms of the zero (equation IV-14) due to the direct differentiation of the target function (equation IV-10).The measurable space of the random matrix  with Then employing the expressions (equation VI-8) and (equation VI-10) and substituting the L2 norm (equation VI-9) we may express the target function ℒ(, ) (equation VI-7) as a compact form of the precision-matrix  precision-matrix,  and its prior (equation.The concavity of the target function ℒ(, ) (equation VI-11) may be analyzed in terms of positive definiteness of the Hessian operator ℋ computed by a block array of the second order matrix derivatives over the arguments  and  (equation VI-12).The block hessian operator (equation VI-12) is positive definite if and only if its diagonal blocks are positive definite.Given the expression of ℒ(, ) in (equation VI-11) we can deduce the structure of the hessian diagonal blocks from the following expressions (equation VI-13) and (equation VI-14). for (equation VI-13) above we obtain the expression in (equation VI-15) and  for (equation VI-14) above we obtain the expression in (equation VI-16).