A methodological framework for inverse-modeling of propagating cortical activity using MEG/EEG

The prevailing view on the dynamics of large-scale electrical activity in the human cortex is that it constitutes a functional network of discrete and localized circuits. Within this view, a natural way to analyse magnetoencephalographic (MEG) and electroencephalographic (EEG) data is by adopting methods from network theory. Invasive recordings, however, demonstrate that cortical activity is spatially continuous, rather than discrete, and exhibits propagation behavior. Furthermore, human cortical activity is known to propagate under a variety of conditions such as non-REM sleep, general anesthesia, and coma. Although several MEG/EEG studies have investigated propagating cortical activity, not much is known about the conditions under which such activity can be successfully reconstructed from MEG/EEG sensor-data. This study provides a methodological framework for inverse-modeling of propagating cortical activity. Within this framework, cortical activity is represented in the spatial frequency domain, which is more natural than the dipole domain when dealing with spatially continuous activity. We define angular power spectra, which show how the power of cortical activity is distributed across spatial frequencies, angular gain/phase spectra, which characterize the spatial filtering properties of linear inverse operators, and angular resolution matrices, which summarize how linear inverse operators leak signal within and across spatial frequencies. We adopt the framework to provide insight into the performance of several linear inverse operators in reconstructing propagating cortical activity from MEG/EEG sensor-data. We also describe how prior spatial frequency information can be incorporated into the inverse-modeling to obtain better reconstructions.


Introduction
The dominant view on the dynamical organization of large-scale cortical dynamics is that it is organized into functional networks of interconnected nodes. This view is widely accepted among researchers working with functional MRI, magnetoencephalography (MEG), and electroencephalography (EEG). It motivates the use of graph theory in analyzing fMRI and MEG/EEG data and is currently the most used approach Bassett and Bullmore (2006) . Cortical activity, however, is also known to exhibit properties that are typically ascribed to spatially-extended systems such as traveling and spiraling waves and activity flow from sources to sinks. We will refer to these sorts of dynamics as propagation . Regarding electrophysiological activity, studies using local field potential recordings have demonstrated that virtually all known oscillations in the mammalian cortex display propagation behavior. This includes induced gamma oscillations in rabbit primary sensory cortices Freeman and Barrie (2000) , beta oscillations in macaque motor cortex Rubino et al. (2006) ; Takahashi et al. (2015) , stimulus-locked responses in cat primary visual cortex Benucci et al. (2007) , traveling waves in macaque visual area 4 triggered by saccadic eye movements Zanos et al. (2015) , and theta oscillations in rat and human hippocampus Lubenov and Siapas (2009) ;Zhang and Jacobs (2015) . Furthermore, several of these studies have shown that the observed wave patterns modulate spiking activity, so that these patterns cannot be explained away as epiphenomena. The studies cited above recorded cortical ac- Ede et al. (2015) ; Patten et al. (2012) . The complicated relationship between propagating waves of cortical currents and the induced topographic patterns recorded by EEG/MEG sensor arrays, however, troubles the interpretation of such studies. Although the EEG/MEG forward model itself is linear, the relationship between the variables of interest in studying propagation dynamics is highly non-linear. This applies for example to phase-dynamics, propagation speeds, and wavelengths. The relationship is further complicated by the geometry of the human cortex. As an example, consider propagation of alpha oscillations. EEG sensor-data give the impression that they propagate between frontal and occipital regions with high speeds and with wavelengths in the order of tens of centimeters. This has given rise to the cortico-cortical propagation hypothesis, according to which propagation of alpha oscillations is mediated by white-matter fibers, which have conduction speeds of about 10 m/s. However, the EEG sensor-data are also consistent with slow propagation of alpha oscillations over small distances, mediated by intra-cortical axons Hindriks et al. (2014) , which have conduction speeds of about 0.3 m/s. According to this hypothesis, the high speeds measured on the scalp are a consequence of the local curvature of the cortex, giving rise to high angular velocities of the electric field at large distances from the cortex Hindriks et al. (2014) .
Another way of investigating propagation of cortical activity in humans is by using electrocorticographic (ECoG) recordings, which are possible in a small group of subjects, for example for presurgical functional localization in patients with intractable epilepsy. Such recordings have been used to investigate propagation of several types of cortical activity, including alpha and gamma oscillations Bahramisharif et al. (2013) ; Halgren et al. (2017) ; Muller et al. (2016) ; Zhang et al. (2018) . Although ECoG sensors are much closer to the cortical surface than EEG sensors and the volume currents do not have to penetrate the highly resistive skull, the convoluted nature of the human cortex makes it difficult to interpret the observed spatiotemporal patterns. Indeed, cortical sulci can be several centimeters deep and are highly curved, hence give rise to the same effects as described above for the scalp EEG alpha oscillations. Even potentials recorded from within cortical tissue have to be interpreted with care as they can be contaminated by volume-conducted currents from other cortical layers and even other regions Schroeder (2011, 2015) . This is especially problematic for investigating propagation dynamics because phases are more susceptible to volume-conduction than amplitudes are. In particular, phases exhibit phase-contraction , which refers to the phases of volume-conducted currents being closer to each other than those of the underlying primary currents, which leads to overestimation of propagation speeds Hindriks et al. (2016) .
Yet another way of studying propagation of cortical activity in humans is to apply source-modeling to EEG/MEG sensor-data and subsequently to characterize the propagation behavior of the reconstructed activity. Several studies have applied source-modeling to EEG/MEG data of propagating cortical activity, including spontaneous and TMSinduced sleep slow waves Bernardi et al. (2019) ; Bersagliere et al. (2018) ; Castelnovo et al. (2016) ; Massimini et al. (2007) ; Murphy et al. ( , 2009) ; Siclari et al. (2014) , evoked K -complexes Riedner et al. (2011) , sleep spindles Siclari et al. (2014) , and high-voltage infraslow waves Dennison (2019) . Although the results of such studies are compelling, it is difficult to assess the quality of the reconstructed waves because the true generator locations are usually unknown. In Wennberg and Cheyne (2013) , source-modeling of K -complexes was carried out using high-density EEG and MEG recordings in epilepsy patients for which the true generator locations were known from intracranial recordings. Several different forward models and inverse methods were compared and it was found that all combinations of forward models and inverse methods incorrectly localized K -complexes to deep midline structures. The authors of Wennberg and Cheyne (2013) note that these findings are relevant not only to K -complexes, but to extended superficial cortical sources in general, and that source-modeling of such activity should be conducted with circumspection. Methodolog-ical studies might help in resolving these conflicting experimental findings and, more generally, inform us about the quality with which propagating cortical activity can be reconstructed from EEG/MEG sensordata. Methodological studies can provide information about the effects of relevant experimental factors such as the number of electrodes and the type of montage in EEG, the coil configuration in MEG, features of the source-space and the inverse method, as well as the effects of the dynamics of cortical activity itself such as location, propagation speed, spatial extent, and wavelength. We are, however, not aware of any methodological study that specifically deals with propagating cortical activity.
In this study we describe a methodological framework for sourcemodeling of propagating cortical activity using linear inverse operators. We restrict to linear operators since their resolution properties can be assessed directly using closed-form expressions. The framework is not limited to linear operators, however, and can be extended to iterative inverse operators, although numerical simulations might be required to study their resolution properties. The framework exploits the spherical topology of the cortex to express cortical activity and linear inverse operators in the Fourier basis on the sphere, i.e. spherical harmonics . This basis is more natural than the dipole basis when studying propagation dynamics and spatially-extended activity in general, because the basis vectors correspond to spatial frequencies, or equivalently, wavelengths, instead of discrete cortical locations. Furthermore, the basis vectors are complex-valued and hence elegantly deal with phase-relationships. We define several tools that can be used to analyze cortical activity and their reconstructions, and to assess the resolution properties of linear inverse operators. The angular power spectrum of cortical activity describes how the activity's energy is distributed across spatial frequencies and is extensively used in cosmology to represent cosmic background radiation Miller et al. (1999) . The angular gain and phase profiles of a linear inverse operator describe the operator's spatial filter characteristics as a function of spatial frequency. The gain profile characterizes the operator's sensitivity to cortical activity of different spatial frequencies and the phase profile characterizes the distortion of phase-relationships of cortical activity at different spatial frequencies and hence relates to timing aspects. The angular resolution matrix characterizes signal leakage within and across spatial frequencies.
We will illustrate the use of the framework by evaluating the performance of several linear inverse methods in reconstructing propagating cortical activity from MEG and EEG sensor-data. The inverse operators include the minimum norm estimate (MNE) Hämäläinen et al. (1993) , weighted minimum norm estimate (wMNE) Lin et al. (2006b) , standardized low resolution electromagnetic tomography (sLORETA) Pascual-Marqui (2002) , andHarmony Petrov (2012) . We describe a simple way to simulate complex propagation patterns, which is to specify the current flow from sinks to sources by means of latency maps and evaluate performance by simulating oscillations that originate form a single source and propagate isotropically outwards with a fixed speed of 0.3 m/s, which is the typical conduction speed of intra-cortical axons. We considered wavelengths in the range 40-250 mm, corresponding to oscillations in the delta, theta, and alpha frequency bands. Performance was quantified by the spatial coherence between the true and the reconstructed waves. Coherence amplitude is a measure for the fidelity of the reconstructed waves, without taking into account a possible discrepancy in timing. The latter is measured by the phase of coherence. Our main findings are the following. First, MNE, wMNE, and sLORETA reconstructions have low coherence amplitude over the entire wavelength range. In particular, due to the noisy appearance of their reconstructions, it is difficult to infer the presence of a wave by inspecting the reconstructions. For all four methods, coherence phase is small over the entire range of wavelengths, indicating that wave timing is not distorted. Thus, although the reconstructions of MNE, wMNE, and sLORETA appear noisy, it might still be possible to infer wave features such as the locations of the current source and its propagation speed. Second, and in contrast to the other three methods, Harmony reconstructions have high coherence amplitude over the entire wavelength range and do not have a noisy appearance.
To assess how the reconstructed waves differed from the true waves, we compared their angular power spectra. This showed that MNE, wMNE, and sLORETA overestimate the power of short-wavelength activity. Harmony does not suffer from such a bias and therefore yields smooth reconstructions. To understand why MNE, wMNE, and sLORETA yield biased reconstructions, we inspected their angular gain and phase profiles and angular resolution matrices. The profiles showed that MNE, wMNE, and sLORETA are sensitive to cortical activity with wavelengths in a wide range, thus behaving as broadband spatial filters. Furthermore, their angular resolution matrices showed that long-wavelength activity tends to leak into short-wavelength activity. These two features are responsible for the noisy-appearance and low fidelity of MNE, wMNE, and sLORETA reconstructions. In contrast, Harmony does not overestimate short-wavelength power and yields reconstructions that appear smooth. Its angular gain profile showed that it is sensitive to long-wavelength activity but strongly suppresses short-wavelength activity, and thus behaves as a lowpass spatial filter. We also show that if the angular spectrum of the true waves is used as a prior in Harmony, performance can be further increased, at least in theory.

EEG/MEG Forward modeling
The magnetoencephalographic (MEG) and electroencephalographic (EEG) forward models relate the current-source density within brain tissue to the recorded magnetic field and scalp potential, respectively Hämäläinen et al. (1993) . For the purpose of studying propagating activity, it will be convenient to model the current-source density in the time-frequency domain. Thus, let s ( t, f ) be a p -dimensional vector with complex-valued entries, where p denotes the number of source locations. The j -th entry of s ( t, f ) represents the current-source density at the j -th location at epoch t and frequency f . Furthermore, let the n -dimensional complex-valued vector b ( t, f ) be the magnetic field (in case of MEG) or scalp potential (in case of EEG) at the n sensors at epoch t and frequency f induced by s ( t, f ). The MEG/EEG forward model in the time-frequency domain is then given by where e ( t, f ) denotes measurement noise and L is the n × p MEG/EEG leadfield matrix, which is constructed on the assumption that current dipoles are perpendicular to the local cortical surface Hämäläinen et al. (1993) . A reasonable model for the measurement noise is that it follows a complex-Gaussian distribution with expectation zero and covariance matrix C e . We model the noise by uncorrelated complex-Gaussian random variables with shared standard-deviation sensor ⟨| ( , ) |⟩, where < | Ls ( t, f )| > denotes the absolute value of Ls ( t, f ) averaged over sensors and sensor is a free parameter. Since in our simulations, the absolute value of s ( t, f ) is independent of t (see Section 2.3 ), the standarddeviation as defined above is in fact a constant. Note that L is realvalued, reflecting the instantaneous nature of the MEG/EEG forward model Hämäläinen et al. (1993) . In this study we use the source-space and MEG/EEG leadfield matrices from the sample data-set provided by MNE Software Gramfort et al. (2015) . The source-space consists of cortical meshes of a single subject with a total of 7498 vertices. The MEG scanner was a 306-channel Neuromag vectorview system and the EEG comprised 60 electrodes. Leadfield errors are modeled by perturbing each column of the leadfield matrix (i.e. each leadfield) according to some probability distribution. In Hosseini et al. (2018) perturbations are modelled by uniform distributions on ellipses in n -dimensional Euclidean space. We follow a similar approach using Gaussian distributions. We distinguish between perturbations along the direction of the leadfield, i.e. radial perturbations , and perturbations perpendicular to the direction of the leadfield, i.e. tangential perturbations , as these two types of perturbations have different effects on the quality of the source reconstructions. The radial and tangential perturbations are denoted by the stochastic matrices rad and tan , respectively. Incorporating the leadfield errors into the forward model Eq. (1) gives the following extended forward model: (2) Radial perturbations are modeled by univariate Gaussian distributions along the directions of the leadfields and are assumed to be proportional to the Euclidean norms of the leadfields and independent across leadfields. Thus, the k -th column of rad has the form rad || || , where l k denotes the k -th column of L and || l k || its Euclidean norm, follows a univariate standard Gaussian distribution, and rad is a parameter that controls the strength of the radial errors. A tangential perturbation to the k -th leadfield is obtained by projecting a zeromean n -dimensional isotropic Gaussian perturbation onto the ( − 1)dimensional sphere with radius || l k ||. The strength of the perturbation is controlled by the standard-deviation tan of the Gaussian distribution. Thus, the k -th column of tan has the form where follows an n -dimensional standard Gaussian distribution and the perturbed leadfield is hence given by Note that the norm of the perturbed leadfield equals the norm of the unperturbed leadfield so that the perturbation is indeed tangential. A less ad hoc approach to modeling leadfield errors is to directly perturb the MEG/EEG sensor locations, the source locations and orientations, and the conductivity parameters of the forward model (in the case of EEG). The advantage of the latter approach is that the errors have clear physical interpretations. However, it requires computation of all leadfields for every realization of the errors, which is computationally very expensive, especially when using a large number of realizations as done in this study. Furthermore, performance is effectively determined by the statistics of radial and tangential errors. Because the focus of our study is on the mathematical framework, we will follow the effective approach described above.

EEG/MEG Inverse modeling
We focus on inverse operators that can be represented by linear transformations from sensor-space to source-space. Popular linear inverse operators in the field of EEG/MEG are the minimum norm estimate (MNE) Hämäläinen et al. (1993) , weighted MNE (wMNE) Lin et al. (2006b) , and standardized low-resolution electromagnetic tomography (sLORETA) Pascual-Marqui (2002) . These inverse operators yield reconstructed current-source densities ̂ ( , ) of the form where minimization is done w.r.t. s ( t, f ), ≥ 0 is the regularization parameter, W is the regularization matrix, which is used to incorporate prior assumptions about the source activity, and || || denotes the Euclidean norm. Setting the gradient of the objective function Eq. (5) to zero and solving for ̂ ( , ) gives where the inverse operator L ♯ is given by in which I is the identity matrix and = ( ) −1 is the prior source covariance matrix and the superscript T denotes matrix transpose.
The inverse operators listed above are all of the form specified by Eq. (7) . Specifically, MNE is obtained by taking = , wMNE is obtained by setting = and normalizing the leadfields to have Euclidean norm 1, and sLORETA is obtained by setting = and adjusting the inverse operator by left-multiplication with a diagonal weighting matrix. In addition to these inverse operators, we consider a more recently proposed linear inverse operator called Harmony Petrov (2012) , which expresses C in the basis of real-valued spherical harmonics, which are the basis functions of the spatial Fourier transformation on the sphere. To apply Harmony, a surface-based source-space with spherical topology is required so that cortical activity patterns correspond to scalar fields on the sphere (in fact two fields; one for each hemisphere). Every cortical activity pattern can be represented by a unique superposition of spherical harmonics. As in the case of the Fourier basis on the line, the spherical harmonics are ordered to increasing spatial frequency, i.e. decreasing wavelength. The advantage of this representation is that source reconstructions can be penalized on the basis of their spatial frequency content, which is natural when dealing with spatially continuous cortical activity. To demonstrate that, apart from spherical harmonics, Harmony can also be formulated using graph harmonics, we implemented Harmony using the eigenvectors of the graph-Laplace operator defined on the triangulated cortical surface Shuman et al. (2013) . An advantage of graph harmonics is that they are exactly orthogonal, irrespective of how the cortical surface is sampled and, moreover, they do not require the source-space to have a spherical topology and can therefore be used for more general source-spaces. The prior source-space covariance matrix C , when expressed in the Laplacian eigenbasis, was taken to be the diagonal matrix with j -th diagonal entry equal to exp (− ∕100) , where j indexes the graph harmonics. We will come back to Harmony in Section 3.4 where we consider more general prior covariance matrices for Harmony.

Generative model of propagating cortical activity
The inverse operators listed in the previous section will be evaluated by using the following generative model of cortical activity: where the index j denotes location, s 0 ( t ) is a driving signal, a j ≥ 0 is the amplitude and j the delay at location j . Activity at location j is hence an amplified ( a j > 1) or attenuated ( a j < 1) delayed (with j ) copy of the driving signal s 0 and the propagation of the driving signal is determined by the vectors = ( 1 , ⋯ , ) and = ( 1 , ⋯ , ) , which we will refer to as amplitude and latency maps , respectively. Negative and positive delays correspond to "early " and "late " with respect to the driving signal and sources and sinks correspond, respectively, to local minima and maxima of when regarded as a function of position on the cortical surface. Cortical activity flows or propagates from early to late regions, thereby creating complex spatiotemporal activity patterns. A more general model can be constructed by using multiple driving signals, each with its own amplitude and latency map, and using the superposition principle. See Mitra et al. (2014Mitra et al. ( , 2015 for a similar generative model of propagating cortical activity in BOLD-fMRI. The driving signal depends on the experimental paradigm. For example, resting-state activity might be modeled by using a stationary stochastic process, whereas for induced and evoked responses, transient signals might be used. Although the driving signal could be interpreted as being physiological, for example as subcortical input, the above model does not require such an interpretation and is merely a convenient way to model propagating cortical activity. Fig. 1 A shows an example latency vector that induces complex spatiotemporal dynamics. Green and red regions are "early " and "late ", respectively, to the driving signal and the minima in early regions and maxima in late regions correspond to sources and sinks, respectively. Note that in the Fourier domain, s j ( t ) is given by where s 0 ( f ) denotes the Fourier transform of s 0 at frequency f and ( ) = − 2 . Cortical background activity can be incorporated into the generative model by adding a p -dimensional random vector that follows a complex-Gaussian distribution with expectation zero and specified covariance matrix. In Sections 3.4.2 and 3.4.3 we consider the effect of background activity and we model it by independent complex-Gaussian random variables with standard-deviation source ⟨| ( ) |⟩. The propagation patterns defined above can model arbitrarily complex spatiotemporal activity patterns. To evaluate the inverse operators we focus on focal oscillations of frequency f that propagate isotropically through the cortex with speed v and dissipate linearly with distance from the ignition location (but see Section 3.4.3 ). Thus, if the oscillations originate from location q on the cortical surface, their amplitude at location q ′ is ( ′ ) = exp (− ( , ′ )∕ ) , where > 0 is the characteristic length and d ( q, q ′ ) is the geodesic distance between q and q ′ . The latency at q ′ is modeled as ′ = ( , ′ )∕ and the wavelength is = ∕ . We restrict to critically-damped waves, i.e. waves for which = ∕2 . For such waves, about one wavelength is observable, in contrast with under-and overdamped waves, for which > /2 and < /2, respectively. Fig. 1 B provides an illustration. Note that for very weak damping (i.e. large ), the wave's amplitude is approximately constant throughout the cortex, giving rise to global oscillations that are synchronous throughout the entire cortex. For very strong damping (i.e. small ), the wave looses its internal structure and can hence be approximated by a point source. The average distance between neighboring vertices on the used cortical mesh is about 5 mm with some outliers up to 22 mm. According to Nyquists' theorem, the minimal wavelength that can be simulated is about 40 mm. Furthermore, the largest geodesic distance between two points on the cortex is about 250 mm. We will consider wavelengths in the range 40-250 mm and set = 0 . 3 m/s, which is the typical conduction speed of action potentials in intra-cortical axons. Fig. 1 C shows a snapshot of an isotropically propagating focal oscillation with a wavelength of 100 mm, originating from the middle section of the right insular cortex. Notice that about one wavelength is visible.

Measure for reconstruction quality
We measure the quality of a reconstruction ̂ ( ) of cortical activity s at frequency f by its spatial coherence with the Fourier transform s ( f ) of the true wave at frequency f : where || || denotes the Euclidean norm and the superscript † denotes the complex-conjugate transpose. Note that for linear inverse operators, and in the absence of measurement noise, ̂ ( ) = 0 ( ) ̂ ( ) , where ̂ ( ) = ♯ ( ) is the reconstruction of u ( f ) and that the coherence simplifies to The coherence is a complex number and can hence be written as Coherence amplitude a ( f ) measures the extent to which the reconstructed wave resembles the true wave, while ignoring a possible difference in timing. The difference in timing is measured by the coherence phase ( f ). In the case of oscillatory waves we will divide the coherence phase by 2 . For example, ( ) = means that the true and reconstructed waves are anti-correlated at frequency f . Depending on the research question, other measures might be appropriate. For example, when cortical phase or amplitude dynamics are of interest, spatial phase-coherence and amplitude-envelope correlations, respectively, might be more appropriate.

Angular power spectra
Let z be a cortical activity pattern in the frequency domain, i.e. a complex-valued function defined on the spatially continuous cortex. When using a surface-based source-space with spherical topology, z can A. Example latency map. Green and red correspond to early and late, respectively, relative to a global driving signal. The local extrema within early and late regions correspond to sources and sinks of activity. When viewed as a function of time, cortical activity flows out of sources and into sinks, thereby creating complex spatiotemporal dynamics. B. Amplitudes of isotropically propagating oscillations as a function of geodesic distance from the ignition location for three different dynamical regimes: overdamped (black), underdamped (blue), and critically-damped (red). C. Example of an isotropically propagating focal oscillation originating from the middle section of the right insular cortex. Propagation speed of the oscillations is 0.3 m/s and the wave has a wavelength of = 100 mm. D. Angular power-spectral density of the wave shown in C. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) be split into a left-and right-hemispheric part, say z L and z R , respectively. Using the mapping between the convoluted cortical hemispheres and a standard spherical cortex, z L and z R can be expressed in spherical coordinates ( , ) and represented as a superposition of complex-valued spherical harmonics: where h denotes either left or right, i.e. h ∈ { L, R }, is the spherical harmonic of degree l and order m , and ℎ is the (complex-valued) coefficient corresponding to . The spherical harmonics form a unitary basis for complex-valued scalar fields on the sphere and constitute a generalization of the Fourier basis in flat space. They are given by for = 0 , 1 , 2 , ⋯ and − ≤ ≤ and where denotes the associated Legendre polynomial of degree l and order m . The constant ensures has unit norm. Note that there are 2 + 1 harmonics of degree l and ( + 1) 2 harmonics of degree at most l . Harmonics of degree l have an associated wavelength of = 4 ∕(2 + 1) , where r denotes the radius of the spherical cortex. For the convoluted cortex, this relation only holds approximately, since the mapping between the convoluted and spherical cortices only approximately preserves distances. In practice, the cortex is sampled and hence the number of sampled harmonics that are orthogonal cannot exceed the number of samples, i.e. the dimension of the source-space. We set the maximum degree to 60, thus giving (60 + 1) 2 = 3721 harmonics for each hemisphere, which is just a bit lower than the number of vertices in the left and right cortical meshes.
The angular power spectrum p h of z h is the spherical equivalent of the spatial power spectrum for signals on flat space and is defined as for degrees = 0 , 1 , ⋯ . It quantifies how the energy in the activity pattern z h is distributed across spatial frequencies (i.e. degrees). Optionally, p h can be normalized by the total power to yield the angular power spectral density . Fig. 1 D shows the angular power spectral density of the propagating oscillations from Fig. 1 C. Note that most power is concentrated in a small number of low-degree harmonics. In Section 3.2 we will use angular power spectra to characterize reconstructed waves and to assess how they differ from the simulated waves. We will also use them to characterize signal leakage across hemispheres. This is done by computing the angular power spectra on the hemisphere that is contralateral to the hemisphere from which the waves originate. Since the simulated waves do not propagate to the contralateral hemisphere, any power in the contralateral part of the reconstructions is spurious and due to inter-hemispheric signal leakage.

Angular resolution matrices and profiles
Insight into the performance of an inverse operator can be obtained by analyzing the resolution operator that is associated with the inverse operator and a given leadfield matrix Hauk et al. (2019Hauk et al. ( , 2011 . The resolution operator R , associated with the pair ( L ♯ , L ), where L and L ♯ are given leadfield and inverse operators, respectively, is defined by The resolution operator describes the relationship between sourceactivity s ( t, f ) and its reconstruction ̂ ( , ) through In particular, the reconstruction of a dipole of unit-strength located at voxel k is equal to the k -th column of R , which is known as the k-th pointspread function (PSF). We will refer to the diagonal entry R ii of R as the gain at voxel i . The rows of R are known as cross-talk functions (CTF's) and the k -th CTF describes the contribution of all active dipole sources to the reconstruction of a dipole source at voxel k Hauk et al. (2019Hauk et al. ( , 2011 . The off-diagonal terms of R thus reflect residual mixing of source activity and cause leakage of activity across locations. Thus, if R ij ≠ 0 for i ≠ j this means that activity at voxel j leaks to voxel i . In this sense, the point-spread function at voxel j characterizes the outgoing leakage from voxel j to all other voxels and the cross-talk function at voxel i characterizes the incoming leakage into voxel i from all other voxels. The total outgoing leakage from voxel j is thus obtained by summing the entries of the point-spread function at voxel j and the total incoming leakage into a voxel i is obtained by summing the entries of the crosstalk function at voxel i . Note that the sum of outgoing leakage equals the sum of incoming leakage so that signal is conserved. We will refer to this property as global signal conservation . Whereas all resolution operators conserve signal globally, they do not necessarily conserve signal locally. A resolution matrix is said to locally conserve signal between voxels i and j if the leakage from i to j equals the leakage from j to i i.e. if = . Thus, a resolution matrix conserves signal locally for all voxels i and j if and only if it is symmetric, i.e. if = . In particular, an asymmetric resolution operator cannot preserve signal locally for all voxel-pairs and its anti-symmetric part characterizes the imbalance between outgoing and incoming leakage. The resolution matrix of MNE is symmetric for all ≥ 0, whereas the resolution matrices of wMNE, sLORETA, and Harmony are asymmetric, even in the absence of regularization (i.e. = 0 ). In this study we express the resolution operator in the basis of spherical harmonics. To do this, we first compute the left-and righthemispheric transformation matrices H L and H R , whose columns are the spherical harmonics up to some maximal degree, sampled on the leftand right cortical mesh, respectively, and normalized to unit norm. The maximal degree can be chosen based on the shortest wavelength that is of interest in a given application. Note that † and † are equal to the identity matrix so that the resolution matrix in the spherical harmonic basis is given by where D is a blockdiagonal matrix with blocks H L and H R . Let ( , ) = † ( , ) be the source-activity expressed in the Harmonic basis and let ̂ ( , ) = †̂ ( , ) be the reconstructed source-activity in the Harmonic basis, then Thus, the resolution operator in the Harmonic basis describes the relation between the source-activity expressed in the Harmonic basis and its reconstruction, also expressed in the Harmonic basis and we can hence define PSF's and CTF's in the Harmonic domain. We list several properties of R H which we will need later on. First, R H is complex-valued. Second, R H is not necessarily Hermitian. However, if R is symmetric, then R H is Hermitian and if the full harmonic basis is used, i.e. the number of harmonics equals the number of cortical vertices, than R is symmetric only if R H is Hermitian. Third, the diagonal entries of R H are of the form h † R h , for some complex-valued column vector h . Because R is real-valued, the diagonal entries of R H hence are real-valued if and only if R is symmetric. Thus, only symmetric inverse operators leave the phase of individual harmonics undistorted, up to a possible shift of 180 degrees. Also, PSF's of harmonics of the form 0 , which are the only real-valued spherical harmonics, posses the same symmetry property as the spherical harmonics, which is that − = ( ) * , where * denotes complex-conjugation. From this it follows that the corresponding reconstruction is real-valued. Therefore, since every real-valued cortical pattern can be expressed as a superposition of spherical harmonics of the form 0 , the reconstructed patterns are real-valued as well.
We define the angular gain profile of an inverse operator as the absolute value of the diagonal of R H averaged over orders within every degree. In the same way, one can define the angular phase profile by averaging the (absolute values of) the phases over orders. Together, the gain and phase profiles provide a characterization of the spatial filtering of the inverse operator, where leakage, i.e. off-diagonal terms in R H are ignored. We define the angular resolution matrix of an inverse operator by setting the diagonals of R H to zero and subsequently averaging the absolute values over orders within every degree. The averaging is done separately for each of the four quadrants of R H , corresponding to the left and right cortices (on-diagonal blocks) and the combinations leftright and right-left (off-diagonal blocks). Thus, the on-diagonal blocks are averaged over orders within degrees to obtain the intra-hemispheric angular resolution matrix and in the same way the off-diagonal blocks are averaged to obtain the inter-hemispheric angular resolution matrix. Although the angular profile and resolution matrix contain less information than R H itself, we will see that they provide insight into the properties of linear inverse operators.

Comparative performance of linear inverse operators
In this section we assess the performance of the inverse operators MNE, wMNE, sLORETA, and Harmony in reconstructing cortical propagating activity from MEG/EEG sensor-data. As described in Section 2.3 , we focus on the reconstruction of oscillations that are initiated locally and propagate isotropically with constant speed and consider wavelengths in the range 40-250 mm in steps of 5 mm. This corresponds to activity at frequency f that propagates with speed v and satisfies = ∕ (see Section 2.3 ). Other parameters were set as in Section 2.3 . Because it is largely unknown to what extent propagating activity can be reconstructed using linear inverse operators, our aim is to obtain an upper bound on performance by considering ideal conditions, i.e. no measurement noise, no background activity, no errors in the MEG and EEG forward models, and known dipole orientations. In Section 3.4.2 we consider the effects of these sources of uncertainties. Since no measurement noise is present, we set the regularization parameter of the inverse operators to a very small value ( = 10 −20 ). Performance was measured by the coherence between the true and the reconstructed waves for every ignition location and wavelength and for both MEG and EEG.
We first assess average performance by averaging the coherence amplitudes and (absolute) phases, respectively, over all ignition locations. Fig. 2 A shows the average coherence amplitude of the MEG reconstructions obtained from the different inverse operators as a function of the wavelength. We make a number of observations. First, MNE, wMNE, and sLORETA perform roughly the same, although sLORETA performs slightly better for most wavelengths. In particular, their coherence amplitudes are more or less independent of wavelength and are rather low. In contrast, Harmony performs substantially better over the entire range of wavelengths. Its performance increases for wavelengths up to about 150 mm and then saturates at a value of about 0.85. Fig. 2 B shows the coherence amplitudes for EEG. The curves are more or less the same as for MEG with the difference that the coherence amplitude of Harmony is about 0.2 lower. These observations demonstrate that, under ideal conditions and for all wavelengths considered, MNE, wMNE and sLORETA do not yield accurate reconstructions of isotropically propagating focal oscillations and Harmony performance substantially better. As described in Section 2.4 , coherence amplitude quantifies the similarity between the true and reconstructed waves while ignoring a possible difference in timing. Timing differences are measured by the (absolute) coherence phase. Fig. 2 C shows the average coherence phase of the MEG reconstructions as a function of wavelength. Notice that MNE has zero phase for all wavelengths, which, as discussed in Section 2.6 , directly follows from the symmetry of MNE's resolution matrix. The other operators do distort wave timing to some extent, but the phases are rather small; less than 4% of the oscillation period of the simulated oscillations. So, for example, for alpha oscillations with a frequency of 10 Hz, the difference in timing between the true and the reconstructed waves does not exceed 4 ms. Fig. 2 D shows that the same holds for EEG.
As an illustration, Fig. 2 E shows homologue synchronized propagating waves with a wavelength of 150 mm originating from primary visual cortex. Its reconstructions are also shown. Note that due to their noisy appearance, it is difficult to infer from the MNE and wMNE reconstructions that the true activity pattern is in fact a propagating wave. This also holds for sLORETA but to a less extent. This can be problematic in practice, especially for noisy recordings. In contrast, Harmony clearly suggests the presence of a propagating wave. It does seem, however, that all reconstructions preserve wave timing, which is also indicated by the small average timing differences in Fig. 2 C and D. Thus, even though the reconstructions of MNE, wMNE, and sLORETA have low coherence amplitude, it might still be possible, at least under ideal conditions, to estimate wave features such as wavelength and ignition location.
We now assess the cortical distribution of coherence amplitude and phase, which provides insight into how reconstruction quality depends on the ignition location of the cortical wave. As a measure of the spatial variation in coherence amplitude, we use its standard deviation. Fig. 3 A and B show the spatial variation as a function of wavelength for the different inverse operators for MEG and EEG, respectively. Roughly speaking, for all inverse operators and for both MEG and EEG, variability decreases with increasing wavelength, particularly for Harmony, which implies that for short waves, reconstruction quality is sensitive to wave location. Furthermore, for long waves that propagate undamped throughout large parts of the cortex, location has practically no impact on reconstruction quality. Fig. 3 E shows the cortical distribution of coherence amplitude for MNE and Harmony for waves with wavelength 40 mm. The topographies of the other inverse operators are similar to that of MNE and are not displayed. To emphasize the variability, we subtracted the average values and scaled the adjusted values to fill the entire range of the colormap. Note that the MEG and EEG topographies are rather similar, although MEG has higher values in buried tissue such as the insular cortex. Concerning the difference between MNE and Harmony, note that their topographies are similar in that coherence is relatively high in dorsal regions and relatively low in deep regions such as the insular cortex and the ventral temporal lobe. Furthermore, for MNE the coherence amplitude closely follows the cortical convolutions, being relatively low on the bottom of sulcal floors and relatively high on top of gyral crowns. This can be seen particularly clearly in the intraparietal and central sulci and on top of the pre-and post central gyri. Note that this is a manifestation of surface bias, i.e. the tendency of linear inverse operators to concentrate estimated source power in locations that strongly project to the sensor array. The figure also shows that Harmony displays less surface bias, which is reflected in the relatively high coherence amplitudes within cortial sulci.
To assess how accurate wave timing is reconstructed, we computed the minimum and maximum timing differences between the true and the reconstructed waves. These differences are plotted as a function of wavelength in Fig. 3 D (MEG) and E (EEG). We make several observations. First, MNE has no timing-difference with the true wave, which is again a consequence of the symmetry of the MNE resolution operator. Second, for all inverse operators and for both MEG and EEG, timing differences decrease with increasing wavelength. Third, Harmony severely distorts the timing of short waves, say with wavelengths up to 50 mm, and the distortions range up to about 45% of the oscillation period in the case of EEG. However, for larger waves, the distortions are neglectable.

Angular power spectra of isotropically propagating waves
In Section 3.1 we found that MNE, wMNE, and sLORETA yield lowquality reconstructions, irrespective of the wave's wavelength, whereas Harmony performed substantially better. The reconstructions in Fig. 1 E indicate that the poor performance is due to the noisy appearance of the reconstructions. That is to say, the reconstructions of MNE, wMNE, and sLORETA seem to have too much energy in high spatial frequencies.
The spatial frequency content of the true and reconstructed waves can be conveniently characterized by their angular power spectra, which we will do in this section. We first consider the spectra of the true waves and subsequently analyze in which aspects they differ from those of the reconstructed waves. Fig. 4 A shows the angular power spectra of the simulated waves for each wavelength in the considered range. They were obtained by averaging the logarithms of the angular power spectra of 100 waves with random locations. Such averaging is needed because, in general, angular power spectra are not invariant under rotations of the cortical sphere, so that changing the location of a wave, changes its spectrum to some extent. This is particularly the case for short waves and less or not at all for long waves. Fig. 4 A shows that the spectra have a global maximum whose location (i.e. degree) decreases with increasing wavelength (red ridge). These spectral peaks can be seen more clearly in Fig. 4 B, which shows the spectra of waves with wavelengths of 40, 100, and 200 mm. The peaks are expected and reflect the fact that the waves have a single wavelength. In fact, the theoretical relation between degree and wave-length is = 4 ∕(2 + 1) , where l is the degree, the wavelength, and r the radius of the spherical cortex (see Section 2.5 ). By setting = 8 cm we obtain the white curve in Fig. 4 A. The relationship only holds approximately, however, because the mapping between the spherical and convoluted cortex only approximately preserves distances. Since the fit of the white curve to the spectral peaks is reasonably accurate, we can use it to identify degrees with wavelengths. Also note that degrees are roughly proportional to spatial frequencies since 1/ ≈ l , up to a proportionality constant.
To compare the spectra of the simulated waves with those of the MEG and EEG reconstructions, we averaged the logarithms of reconstructed waves over the same 100 locations as above. The locations were selected from the same hemisphere from which the wave originated. This was done separately for all wavelengths. The obtained spectra are shown in Fig. 4 C for MEG and in Fig. 4 D for EEG. We first consider MNE, wMNE, and sLORETA. Comparing Fig. 4 C and Fig. 4 A shows that for all wavelengths, MNE and wMNE severely overestimate the wave's power at all but the lowest degrees (say < 10). For short waves, overestimation starts at a higher degree (say > 20). These observations hold for sLORETA as well, although to a lesser extent. Also notice that the energy at very high degrees is slightly underestimated for all wavelengths. For Harmony, overestimation of high-degree power is practically absent. All four methods, however, seem to preserve the relation between wavelength and degree (red ridges). Fig. 4 D shows that the spectra of the EEG reconstructions are very similar to those of the MEG reconstructions, but do not preserve the relation between wavelength and degree as well as in the case of MEG.
In the above analysis we considered the ipsilateral spectra of the reconstructions, i.e. the spectra of the hemisphere in which the wave originated. We now consider their contralateral spectra . Note that the contralateral spectra of the true waves are zero because the waves are confined to a single hemisphere. Hence, any power in the contralateral part of the reconstructions is due to inter-hemispheric signal leakage and is therefore spurious. Because the contralateral spectra of the reconstructions were found to be approximately independent of wavelength (results not shown), we averaged them over wavelengths to obtain a more concise representation. Fig. 4 E and F show the contralateral spectral densities for all inverse operators for MEG and EEG, respectively. Instead of using the total contralateral power to normalize the spectra, we used the total ipsilateral power as this allows to assess the severity of signal leakage. The negative values on the vertical axis thus imply that the reconstructions contain more power in the ipsilateral than in the contralateral hemisphere; a conclusion that could have been drawn more straightforwardly by inspecting the reconstructions in the dipole bases. The spectra in Fig. 4 E and F, however, provide information about the spatial frequencies through which signal leaks away to the other hemisphere. For example, for both MEG and EEG and for all inverse operators, signal leakage is strongest for low spatial frequencies and gradually weakens with increasing frequency. There are also notable differences between the operators. First, for low frequencies, sLORETA is more susceptible to signal leakage than MNE and wMNE. Second, Harmony is least susceptible to signal leakage over the entire frequency range and spurious power decreases much faster with increasing frequency than that of the other inverse operators.

Angular profiles
In the previous two sections we discussed aspects related to the quality of reconstructions in the case of a particular generative model of propagating cortical activity, namely isotropic propagation of local oscillations. We now turn to the inverse operators themselves and analyze their angular profiles.
We first consider the gain profiles for MEG, which are shown in Fig. 5 A. The horizontal axis corresponds to spatial frequency in cy-  cles/m, the values of which were obtained by exploiting the correspondence between spatial frequency and harmonic degree (see Section 3.2 ). Notice that the profiles of MNE and wMNE are very similar, implying that normalization of the leadfields, which is how wMNE is obtained from MNE, has almost no effect on the gain profile and therefore does not alter the filter characteristics of MNE. The profiles show that MNE and wMNE are sensitive to cortical activity at all spatial frequencies, at least to some extent, but are particularly sensitive to activity with spatial frequency at about 30 cycles/m. In contrast, Harmony acts more like a lowpass spatial filter with a cut-off frequency of about 25 cycles/m. This implies that cortical activity with spatial frequencies exceeding 25 cycles/m is suppressed, whereas activity with frequencies below 25 cycles/m passes through with high gain. This lowpass character makes that Harmony reconstructions are so much smoother than those of the other inverse operators. The gain profile of sLORETA is similar to that of MNE and wMNE, but is a bit tilted: it is higher for low frequencies and lower for high frequencies. The figure suggests that the gain profile of sLORETA is somewhere in between that of MNE and wMNE on the one hand, and Harmony on the other hand. Notice also that sLORETA is relatively sensitive to very low frequencies, say below 5 cycles/m. Fig. 5 B shows the gain profiles for EEG. The figure shows that the relation between the profiles of the different operators is the same as in the case of MEG in that MNE and wMNE have roughly the same profile and that sLORETA, and particularly Harmony, are more sensitive to lower frequencies and less sensitive to higher frequencies. Distinct differences between MEG and EEG can also be observed, however. First, compared to MEG, EEG is less sensitive to high-frequency activity than MEG, which is evident from the steeper roll-off of the profiles above 30 cycles/m. A practical consequence of this is that in the presence of broadband cortical background activity or measurement noise, high-frequency activity will be more difficult to reconstruct from EEG sensor data. A second difference between MEG and EEG is that the peak-sensitivity of MNE, wMNE, and sLORETA is a bit lower than for MEG: about 25 cycles/m for EEG against about 30 cycles/m for MEG. The cut-off frequency of Harmony seems to be shifted to lower frequencies by the same amount.
The differences between the gain profiles for MEG and EEG noted above apply to all considered inverse operators and are therefore likely to reflect the spatial filtering characteristics of the MEG and EEG forward models. That this is indeed the case can be seen in Fig. 5 E, which shows the gain profiles of the MEG (blue curve) and EEG (red curve) leadfield matrices. The profiles were obtained by expressing the MEG and EEG leadfield matrices in the harmonic basis and subsequently taking the Euclidean norms of the columns (i.e. leadfields). The profiles show that EEG is more sensitive to very low frequencies and less sensitive to high frequencies and that the peak-sensitivity of MEG is about 30 cycles/m, whereas that of EEG is about 25 cycles/m, in line with Fig. 5 A and B. These frequencies correspond to wavelengths of about 3-4 cm, which is about the typical depth of cortical sulci. We suspect that MEG and EEG are less sensitive to cortical activity with wavelengths slightly shorter than 3-4 cm because the activity can be completely buried within sulci, in which case the activity will be erroneously projected to surrounding gyri due to surface bias. Furthermore, MEG and EEG might be less sensitive to activity with wavelengths slightly longer than 3-4 cm because of cancellation of synchronous activity on opposite walls of cortical sulci Ahlfors et al. (2010) ;Hauk et al. (2011) ;Hindriks et al. (2014) . Thus, the peak-sensitivities of MEG and EEG might be directly related to cortical geometry.
Fig. 5 C and D show the angular phase profiles in the case of MEG and EEG, respectively. For each spatial frequency, the phase was obtained by averaging the absolute phases over all orders within the corresponding harmonic degree and were subsequently averaged over the left and right hemispheres (see Section 2.6 ). Angular phase profiles hence measure the average phase distortion of true cortical activity by a linear inverse operator, as a function of spatial frequency. For oscillatory waves, phase distortion translates into timing differences between true and reconstructed waves. Note that MNE has zero phase-distortion for all spatial frequencies. As discussed in Section 2.6 , this is a direct consequence of the symmetry of MNE's resolution matrix. Having asymmetric resolution matrices, wMNE and sLORETA do display phase-distortion, but these distortions are small: they do not exceed 0.05 radians, which, for a spherical cortex with radius 8 cm, corresponds to 4 mm. In contrast, Harmony severely distorts the phase of cortical activity for spatial frequencies > 20cycles/m, i.e. above its cut-off frequency. Below 20 cycles/m, its distortion is comparable to that of wMNE and sLORETA. In the next section we will see that the level of distortion reflects the extent to which the resolution operators are asymmetric.

Angular resolution matrices
In this section we analyse the intra-hemispheric angular resolution matrices of the inverse operators of MNE, wMNE, sLORETA, and Harmony. Fig. 6 A shows these matrices in the case of MEG. They were obtained by averaging over hemispheres, dividing the entries by the largest absolute entry, and subsequently taking the logarithm. The normalization allows to compare the resolution matrices across inverse operators as well as between MEG and EEG. The ( i, j )-th entry of the resolution matrix measures the leakage from cortical activity with spatial frequency f j to cortical activity with frequency f i , where f i and f j are the spatial frequencies corresponding to the i -th and j -th columns (or rows), respectively. In particular, the diagonal entries measure the leakage of activity within the different frequencies. The i -th column of the resolution matrix, i.e. the i -th point-spread function, thus contains the leakage of activity out of frequency f i and into all other frequencies, whereas the i -row, i.e. the i -th cross-talk function, measures the leakage of activity out of the different frequencies and into frequency f i . Fig. 6 A shows, for example, that wMNE leaks activity across frequencies more than MNE does, which is evident from the red region surrounding the diagonal, which is wider than for MNE. It also shows that sLORETA leaks strongest for low spatial frequencies and that, for all inverse operators, there is almost no leakage from high to high spatial frequencies. This is because there is relatively less activity to leak, due to the low gain at high frequencies.
We also note that Harmony leaks activity from all spatial frequencies into low spatial frequencies.
To better understand how cortical activity leaks across frequencies, we separately consider the symmetric and anti-symmetric parts of the resolution matrices as shown in Fig. 6 B and C, respectively. Note that MNE has no anti-symmetric part because its resolution matrix is symmetric. Its symmetric part is hence equal to the resolution matrix itself. The symmetric part of a resolution matrix characterizes that part of the leakage that is locally conserved, i.e. the leakage from f i to f j equals the leakage from f j to f i for all i and j . This implies that there is no net leakage into or out of any frequency. The part of the leakage that is not locally preserved is contained in the anti-symmetric part of the resolution matrix. Fig. 6 C shows that wMNE, sLORETA, and Harmony have such an anti-symmetric part, which reflects the asymmetry of their resolution matrices. The figure shows that wMNE, sLORETA, and Harmony leak broadband activity to lower frequencies, although to different extents (compare the color bars). For example, cross-frequency leakage is about 10 times stronger in Harmony than in wMNE. A convenient way of summarizing how signal is leaked across spatial frequencies is by plotting the column-averages of the anti-symmetric part of the resolution matrix as a function of spatial frequency. This gives the net incoming leakage (positive value) or net outgoing leakage (negative value) as a function of spatial frequency. Fig. 7 shows the net leakage for the different inverse operators for MEG (left panel) and EEG (right panel). Since the total leakage is zero (see Section 2.6 ), the net leakage, when summed over all frequencies, is zero as well. Note that MNE has no net leakage, again reflecting that it only leaks activity within frequencies. The other three inverse operators have net leakage into low frequencies, say up to about 40 cycles/s and net leakage out of frequencies higher than 40 cycles/s. This holds for both MEG and EEG. The net leakage of Harmony, however, is much stronger than that of MNE, wMNE and sLORETA.

Harmonic inverse-modeling of propagating waves
In the previous sections we found that among the tested linear inverse operators, Harmony behaved rather differently from the other operators and yielded markedly better wave reconstructions. In this section we provide a more detailed performance assessment of Harmony. In particular, we assess to what extent the performance of Harmony can be further increased by matching its spatial filtering characteristics to the frequency content of the simulated cortical activity (Section 3.4.1), quantify the effects of measurement noise, cortical background activity, and inaccuracies in the MEG/EEG leadfield matrices on Har-  mony's performance (Section 3.4.2), and use Harmony to reconstruct simulated complex spatiotemporal dynamics comprising multiple overlapping propagating sources (Section 3.4.3). In Sections 3.4.2 and 3.4.3 we use the same prior as used in earlier sections and in order to reduce the computational load, we restricted the left-and right-hemispheric Harmonic transformation matrices to the first 400 columns. The latter confines the wave reconstructions to the linear subspace spanned by the first 400 eigenvectors of the graph Laplacian operator.

Harmony with optimal prior angular power spectrum
In the original formulation of Harmony Petrov (2012) , the prior source-space covariance matrix is expressed in the basis of real-valued spherical harmonics and was chosen to be a diagonal matrix. The entries on the diagonal are constant within degrees and equal to 1∕(1 + ) , where l denotes harmonic degree and where ≥ 0 is a hyperparameter. For = 0 this yields the MNE inverse operator and for > 0, the inverse operator implements a spatial low-pass filter with cut-off frequency that is controlled by . As described in Section 2.2 , we used a graph-theoretical version of Harmony with diagonal entries (− ∕100) , where j indexes the graph harmonics. With this prior, an angular power spectrum can be associated by expressing it in the basis of complex-valued spherical harmonics. This spectrum is shown in 8A (blue curve). The same can be done with the original formulation of Harmony. In fact, because the prior 1∕(1 + ) is constant within degrees, it can be interpreted as an angular power spectrum: Harmony can hence be implemented with any prior angular power spectrum so that, generally, the prior spectrum can be written as p ( l ), where is a vector of hyperparameters. Alternatively, the prior power spectrum can be treated non-parametrically using suitable basis functions or obtained from an independent source of data, if available. In particular, we can choose it to be equal to the angular power spectrum of the true cortical activity. In practice, however, the true spectrum is generally unknown and hence the prior spectrum needs to be chosen in a data-driven way, for example using cross-validation. Moreover, it is unclear if this can be done at all, because information is lost in the forward mapping from source-to sensor-space. Indeed, in the original study Petrov (2012) p was set by hand so this question remains open even in the case of a single hyperparameter. Irrespective, we can assess to what extent reconstructions of propagating cortical activity improve if the prior angular power spectrum is chosen optimally, i.e. equal to the true angular power spectrum. Below, we address this question for the special case of isotropically propagating focal oscillations.
We focus on waves with wavelength = 40 mm because the performance was lowest for this wavelength (see Fig. 2 A and B). The angular power spectrum of such waves is shown in Fig. 8 A (green curve) and was obtained by averaging the spectra of 100 waves with random locations. Using this spectrum to construct the prior source covariance matrix in Harmony, we computed the coherence amplitudes of the reconstructions for waves at all cortical locations. In Fig. 8 B, the resulting values are plotted against those obtained in Section 3.1 , i.e. using the regular prior. The figure shows that for both MEG and EEG, the optimal prior yields higher coherence amplitudes for most locations, although the increase is modest. Specifically, average coherence amplitudes for regular and optimal Harmony are 0.43 and 0.52, respectively, for MEG and 0.19 and 0.26, respectively, for EEG. The spatial correlation between the regular and optimal coherence amplitudes was 0.88 for both MEG and EEG, which shows that performance increases more or less uniformly throughout the cortex. For MEG, however, the largest increases are observed around the occipital pole, where the leadfields are strong (results not shown). As an illustration, we selected the location which diplayed the largest increase in performance. The wave and its reconstructions are shown in Fig. 8 C. It is clear that the optimal prior improves the reconstruction. Fig. 9 A shows the gain profiles of Harmony with regular (blue curve) and optimal (green curve) prior in the case of MEG. The profile corresponding to the regular prior we have already encountered in Section 3.3 (see Fig. 5 A). The green curve shows that Harmony with optimal prior is highly sensitive to cortical activity at those frequencies that have high prior variance, i.e. activity with wavelengths about = 40 mm. Fig. 9 B shows the angular resolution matrices of Harmony with regular and optimal prior. Note that Harmony with optimal prior leaks activity from all frequencies into frequencies around the frequency of the simulated waves. Together, these spatial filtering and leakage properties account for the improved performance observed in Fig. 8 B.

Effects of measurement noise and leadfield errors
The results of the previous sections were obtained under ideal conditions in which no cortical background activity, measurement noise, or inaccuracies in the leadfields were present. Although these assumptions enable quantifying the theoretical performance of the different inverse methods and thereby allowing to establish an upper limit on performance, quantifying the effect of errors on performance provides valuable information for practical applications. As described in Sections 2.1 and 2.3, the strength of sensor noise and cortical background activity is controlled by the parameters sensor and source , respectively. Furthermore, we consider radial and tangential inaccuracies in the leadfields, the magnitude of which is controlled by the parameters rad and tan , respectively. In this section we assess the performance of Harmony as a function of these parameters in reconstructing hemispherically-synchronized oscillations with a wavelength of = 150 mm that originate from primary visual cortex (see Fig. 2 E).
To assess the effect of measurement noise, we considered values of sensor between 0 and 1 in steps of 0.1 and computed 100 reconstructions using Harmony as well as their coherence with the true waves. The values of the regularization parameter of the Harmony inverse operator (see Section 2.2 ) were determined by generalized cross-validation. Fig. 10 A shows the mean and standard error of the coherence amplitude as a function of the level of measurement noise sensor for both MEG and EEG. For both EEG and MEG, the coherence amplitude rapidly decreases when noise is present and keeps gradually decreasing for increasing noise-levels. The relatively large initial decrease from sensor = 0 to sensor = 0 . 1 is due to the large increase in the regularization parameter of the inverse operator, which introduces additional bias in the wave reconstructions. The gradual decrease for higher levels of sensornoise is due to the combination of further increased regularization and higher noise-levels. Fig. 10 A also shows that the standard error for EEG is higher that that for MEG. For each noise-level, the value of the regu- larization parameter, as selected by generalized cross-validation, turned out to be almost constant across realizations for MEG, but relatively variable for EEG (results not shown), which contributes to the relative high standard error of EEG. Another contribution comes from the relatively larger variability of the EEG reconstructions themselves, which is due to the larger condition number of the EEG leadfield matrix and therefore reflects the different physical properties of the EEG and MEG forward models. Fig. 10 B shows the mean and standard errors of the (absolute value of) coherence phase as a function of noise-level for both EEG and MEG. For both EEG and MEG, the coherence phases gradually increase with noise-level to about 10% of the oscillation period (on average), although for EEG the phases are more variable across realizations. We note, however, that the estimated phases are only useful if the coherence amplitude is high and otherwise are randomly distributed across the cortex. The observations in Fig. 10 A and B indicate that in the presence of measurement noise, propagating waves can be more accurately reconstructed with MEG than with EEG, both with respect to fidelity ( Fig. 10 A) and wave timing ( Fig. 10 B). To assess the effect of cortical background activity, we selected values of source between 0 and 1 in steps of 0.1 and computed 100 reconstructions using Harmony as well as their coherence amplitudes with the true waves. Fig. 10 C shows the mean and standard error of the coherence amplitude of the level of background activity for both MEG and EEG. Observe that for both EEG and MEG, the initial decrease in performance with increasing level of background activity is more gradual than in the case of measurement noise ( Fig. 10 A). Another difference is that EEG seems to be more robust against background activity than MEG. Fig. 10 D shows that for both EEG and MEG, the coherence phase remains small even for high levels of background activity. The observations in Fig. 10 C and D indicate that in the presence of moderate to strong levels of cortical background activity, propagating waves can be more accurately reconstructed with EEG than with MEG.
To assess the effect of radial errors, we varied their strength rad between 0 and 1 in steps of 0.1, and for each value, computed 100 reconstructions using Harmony, together with the corresponding coherence amplitudes. Fig. 11 A shows the mean and standard errors of the coherence amplitude as a function of the level of radial noise for both EEG and MEG. The MEG reconstructions are clearly more sensitive to radial errors than the EEG reconstructions. In fact, the performance for EEG is nearly constant across the entire range of noise-levels and the same is true for the coherence phases as can be seen in Fig. 11 B. This difference in performance is likely due to the stronger lowpass filtering characteristics of the EEG forward model, which makes the reconstructions relatively robust to model errors. Similar observations can be made about the tangential errors, as shown in Fig. 11 C and D, which were obtained by varying the tangential noise-level tan between 0 and 0.01 in steps of 0.001. Note that the levels of radial and tangential errors cannot be directly compared since the former depends on the norm of the leadfields, whereas the latter only depends on the leadfields' orientation (see Section 2.1 ). Since tangential errors are related to errors in the estimates of dipole orientations, these observations likely reflect the relative sensitivity of MEG to dipole orientations. The observations in Fig. 11 C indicate that the wave reconstructions obtained using EEG are more accurate and robust than those obtained using MEG.

Inverse-modeling of complex spatiotemporal dynamics
In the preceding sections we have modeled cortical activity either by one or two propagating sources. Although these scenarios might be relevant for certain types of task-based experiments, generally, and certainly in non-structured conditions such as sleep, general anesthesia, coma, and the resting-state, a multitude of propagating sources is likely to be active. It will therefore be interesting to quantify the performance of Harmony in case of multiple propagating sources. We simulated cortical activity by placing ten isotropically and critically-damped propagating sources in each hemisphere in random locations. Keeping these locations fixed, we generated 100 cortical activity patters by randomly varying the phase off-sets of each of the sources. This was done by multiplying the frequency-domain representation s ( f ) of each source by a factor exp ( i ), where is uniformly distributed in the interval [0, 2 ]. Wavelengths and propagation speeds of the sources were set to 60 mm and 0.3 m/s, respectively. The superposition of these randomly-timed propagating sources gives rise to complex spatiotemporal patterns in which activity propagates between multiple sources and sinks.
We conducted two sets of simulations: one in the absence of and one in the presence of measurement errors, cortical background activity, and leadfield inaccuracies. We refer to these simulations as "ideal " and "noisy ", respectively. In the noisy case, we set sensor = 0 . 2 , source = 0 . 2 , rad = 0 . 1 , and tan = 0 . 001 . The values of the regularization parameter of the Harmony inverse operator (see Section 2.2 ) were determined by generalized cross-validation. In the ideal case, the mean and standard error of the coherence amplitude over the 100 realizations were 0.70 ± 0.02 for MEG and 0.40 ± 0.03 for EEG. We computed the means and standard errors for 20 different sets of cortical locations and found that the values diplayed above are fairly representative. In the noisy case, the mean and standard error of the coherence amplitude over the 100 realizations were 0.35 ± 0.03 for MEG and 0.22 ± 0.04 for EEG. These observations indicate that the reconstructions obtained with MEG are more sensitive to noise and errors than those obtained with EEG. Fig. 10. Effect of measurement noise and background activity. A. Mean and standard errors of the coherence amplitude as a function of sensor for both MEG (blue) and EEG (red). B. Mean and standard errors of the absolute value of the coherence phase as a function of sensor for both MEG (blue) and EEG (red). C. Mean and standard errors of the coherence amplitude as a function of source for both MEG (blue) and EEG (red). B. Mean and standard errors of the absolute value of the coherence phase as a function of source for both MEG (blue) and EEG (red). The curves in panels A, B, C, and D were obtained by simulating 100 independent realizations of the EEG/MEG forward model, obtaining the source reconstructions using Harmony, and computing their spatial coherence with the true waves. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Wave-timing, as quantified by the coherence phase, is a global quantity in the sense that the phases of the true and reconstructed waves at all cortical locations contribute to its value. The simulated activity patterns, however, cover the entire cortex and hence have a timing (i.e. a phase) at every location and it is therefore interesting to see how well local timing can be estimated. For this, we pooled the local phase-differences between the true and reconstructed waves of all cortical locations and 100 realizations, resulting in the overall distribution of phase-differences shown in Fig. 12 A. The figure shows that under ideal conditions and for both MEG and EEG, the phase-differences concentrate around zero, which means that, for most cortical locations, local timing is accurately estimated. However, both densities cover the entire circle and hence for a substantial number of cortical locations the timing is inaccurately estimated. This is particularly the case for EEG. The figure also shows that, under noisy conditions, the distributions become more uniform, indicating that local wave-timing is not accurately estimated. Figs. 12 B and C provide illustrations for the ideal and noisy cases, respectively. Shown are the cortical distributions of the real-part, absolute value, and phase of a randomly simulated wave (top row), together with their reconstructions obtained with MEG (middle row) and EEG (bottom row). The difference in accuracy between MEG and EEG of local wave-timing under ideal conditions is clearly visible (third column of Fig. 12 B). Also notice the lack of correspondence between the true and the estimated local phases in the noisy case, particularly for EEG (third column of Fig. 12 C).
Figs. 12 B and C seem to suggest that the relatively poor performance of EEG is due, at least partially, to a "blurring " of the true activity patterns in the EEG reconstructions. This is particularly clear when inspecting the real parts and absolute values of the reconstructions (first and second columns of Figs. 12 B and C). To confirm this, we computed the angular power spectral densities of the true and the reconstructed waves. They are displayed in Fig. 12 D. Notice that whereas the true activity pattern has a characteristic frequency of about 15 cycles/m, the EEG reconstructions in the ideal case (solid red line) have a lower char- Fig. 11. Effect of leadfield errors. A. Mean and standard errors of the coherence amplitude as a function of rad for both MEG (blue) and EEG (red). B. Mean and standard errors of the absolute value of the coherence phase as a function of rad for both MEG (blue) and EEG (red). C. Mean and standard errors of the coherence amplitude as a function of tan for both MEG (blue) and EEG (red). B. Mean and standard errors of the absolute value of the coherence phase as a function of tan for both MEG (blue) and EEG (red). The curves in panels A, B, C, and D were obtained by simulating 100 independent realizations of the EEG/MEG forward model, obtaining the source reconstructions using Harmony, and computing their spatial coherence with the true waves. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) acteristic frequency (between 5 and 10 cycles/m). The MEG reconstructions, on the other hand, have a characteristic frequency that matches that of the true activity patterns. Furthermore, under noisy conditions, the characteristic frequencies of both the MEG and EEG reconstructions decrease, which is due to the increased regularization of Harmony's inverse operator. These observations reflect the different spatial filtering properties of the EEG and MEG forward models. Coming back to Section 3.4.1 , the reconstructions might be improved by choosing a different prior angular power spectrum for the inverse operator. As mentioned in Section 3.4.1 , an interesting direction for future research will be to investigate to what extent prior angular power spectra can be adaptively chosen, i.e. based on the observed data, since this has the potential of substantially increasing the fidelity of the reconstructions.

Discussion
In this study we provided a framework for the analysis and performance evaluation of MEG/EEG linear inverse methods with respect to the reconstruction of propagating cortical activity. We conducted numerical simulations to assess the performance of several linear inverse methods, namely Hämäläinen et al. (1993) , wMNE Lin et al. (2006b) , sLORETA Pascual-Marqui (2002 , and Harmony Petrov (2012) , in reconstructing local propagating oscillations. Furthermore, we have introduced angular power spectra, angular gain and phase profiles, and angular resolution matrices to explain their performance and to provide a general framework for analyzing linear inverse operators in the spatial frequency domain. With respect to performance, our main conclusion is that only Harmony yields reconstructions with high fidelity. We adopted the proposed framework to gain insight into these findings and found that MNE, wMNE, as well as sLORETA, albeit to a lesser extent, overestimate the power of short-wavelength cortical activity, thus yielding reconstructions that appear noisy. The angular gain profiles of MNE, wMNE, and sLORETA showed that these methods are sensitive to cortical activity with wavelengths in a wide range and hence act as broadband spatial filters. Furthermore, their angular resolution matrices showed that long-wavelength cortical activity tends Fig. 12. Harmonic inverse-modeling of complex spatiotemporal dynamics. A. Distribution of local phase-differences between the true and reconstructed waves for MEG (blue) and EEG (red) obtained under ideal (solid lines) and noisy (dashed lines) conditions. B. True (top row) and reconstructed (middle and bottom rows) activity patterns obtained by randomly selecting one of the 100 realizations under ideal conditions. Shown are the real-part (first column), the absolute value (middle column), and phase (right column) of the activity patterns. C. Same format as B. but onder noisy conditions. D. Angular power spectral densities of the true (black line) activity patterns and their reconstructions (MEG (blue) and EEG (red)) under ideal (solid lines) and noisy (dashed lines) conditions. The densities were obtained by averaging over all 100 realizations. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) to leak into short-wavelength activity. These two features are responsible for the noisy-appearance and low fidelity of the MNE, wMNE, and sLORETA reconstructions. In contrast, Harmony does not overestimate short-wavelength power and yields reconstructions that appear smooth. Its angular gain profile showed that Harmony is sensitive to long-wavelength activity, but strongly suppresses short-wavelength activity and hence acts as a lowpass spatial filter. Thus, MNE, wMNE, and sLORETA yield low-fidelity reconstructions because their spatial filter characteristics do not match the spatial frequency content of the simulated propagating activity. The spatial filter characteristics of Harmony better match the frequency content of the simulated activity and hence yields high-fidelity reconstructions.
In the evaluation we have applied the inverse methods to both MEG and EEG and found that the reconstructions based on MEG have higher fidelity than those based on EEG. Although it remains an open question to what extent this is due to the relatively low number of EEG electrodes, we observed that the filter characteristics of the inverse methods largely reflect those of the respective forward models. To further illustrate the use of the framework, we showed that if the angular power spectrum of the cortical activity is known, or can be inferred from the sensor data, it can be used to further increase the performance of Harmony. This holds at least theoretically and for isotropically propagating focal oscillations and should therefore only be considered a proof of principle. Note that in constructing this prior, nowhere it is assumed that the location of the wave is known. In sum, our study provides a framework that generalizes resolution analysis Hauk et al. (2019Hauk et al. ( , 2011 to the spatial frequency domain, and might be helpful in the development and assessment of source-modeling methods for propagating cortical activity. There are several issues that have not been addressed in this study, but which are nevertheless relevant for practical applications. A first issue is that by measuring the spatial coherence between the true and reconstructed cortical propagation patterns, we have addressed only one aspect of performance. Being a complex-valued quantity, the spatial coherence takes both amplitudes and phases into account. Thus, if the amplitudes of the true and reconstructed patterns are uncorrelated, coherence will be zero, even if the patterns are completely phase-locked. Likewise, if their phases are uncorrelated, coherence will be zero, even if the patterns have the same amplitude distributions. We do remark, however, that if the amplitudes are not accurately reconstructed, one might expect the phases to be inaccurately reconstructed as well, because phases are more susceptible to signal leakage than amplitudes are Hindriks et al. (2016) . In any case, if the reconstruction of amplitudes or phases is the main interest, instead of spatial coherence, spatial amplitude correlation or the phase-locking value, respectively, could be used for evaluation. A second issue that we did not address is that, in practice, particular features of cortical activity are often of interest, such as the number and location of sources, propagating speeds, wavelengths, and timing relationships. Such features might still be inferred from the reconstructions, even though the reconstructions themselves have low fidelity. Thus, although MNE, wMNE, and sLORETA yield low-fidelity reconstructions, the absence of phase-distortion suggests that they might still enable accurate estimation of some of these features. Other features that might be of interest are sources and sinks and flow patterns Takahashi et al. (2015) ; Townsend et al. (2015) . Currently, it is unknown to what extent such features can be reconstructed from MEG/EEG sensor-data. The most pressing issue with respect to practical applications, however, is the estimation of dipole orienta-  Murphy et al. (2009) . Knowing the dipole orientations will give access to phase-information, which is necessary to analyze timing aspects. Dipole orientations can be estimated in several ways, for example, directly from the vector-valued time-courses of the reconstructed sources Uutela and Ha (1999) or by exploiting anatomical information Dale et al. (2000) ; Lin et al. (2006a) , among several other methods (see Rubega et al. (2019) and references therein). These methods, however, are local in the sense that the estimated orientations are not explicitly required to be consistent across the cortex, hence leading to spatial discontinuities in the estimates. This is especially problematic in the context of propagating cortical activity, which is characterized by continuously varying phase-fields.
The framework laid out in this study is based on the spherical Fourier transformation and its associated basis functions, i.e. the complexvalued spherical harmonics. There are three potential drawbacks of using this basis. First, the use of spherical harmonics is enabled by exploiting the one-to-one mapping between the convoluted cortex and the sphere. Although this mapping preserves topology, it generally does not preserve distances, although such mappings are designed to minimize metric distortions Fischl (2012) . Consequently, although the spatial frequencies associated with spherical harmonics are exact when viewed on the spherical cortex, they are only approximate when viewed on the convoluted cortex. Second, the distortions are subject-specific and this could create inter-subject variation in angular quantities that does not reflect variation in function but merely relates to cortical geometry. We remark, however, that inter-subject differences in anatomy can confound any group-level functional analysis. Third, due to sampling of the sphere, the harmonics are not guaranteed to be exactly orthogonal and the deviations increase with increasing degrees. This potentially leads to inaccuracies in angular quantities such as angular power spec-tra, gain and phase profiles, and resolution matrices. This issue is not fundamental, however, and can be dealt with by appropriate sampling of the sphere Khalid et al. (2014) .
An alternative to using the spherical Fourier transformation is to define basis functions directly on the convoluted cortex. Since the Fourier basis functions are defined as the eigenvectors of the Laplace operator, and any smooth manifold admits a Laplace operator, any smooth manifold has a Fourier basis. In practice, the convoluted cortex is sampled and triangulated, yielding a graph so that one can take the eigenvectors of the graph Laplace operator, i.e. the graph harmonics , as basis functions Shuman et al. (2013) . The advantage of graph harmonics is that, irrespective of how the cortex is sampled, they are exactly orthogonal and do not require a one-to-one mapping to a reference space. Using graphharmonics has drawbacks as well, however. First, it is not clear how angular quantities can be defined, because there is no natural grouping of graph harmonics by degree. Second, if the harmonics are to have a physical interpretation, the weighted graph-Laplacian needs to be used, because the unweighted, i.e. combinatorial, graph-Laplacian is invariant under metric distortions and hence only contains topological information. As the weighted graph-Laplacian depends on the geometry of the cortex, however, it is not immediately clear how a group-level analysis can be carried out. Third, since graph harmonics are real-valued, straightforward analysis of timing aspects in source-modeling might be difficult.
We have shown that, at least for isotropically propagating focal oscillations, the performance of Harmony can be improved if the prior angular power spectrum is chosen to match that of the true wave. In practice, the true angular power spectrum is usually unknown and therefore has to be inferred from the sensor-data. Once a parametric model for the angular spectrum is assumed, the problem reduces to that of selecting the hyperparameters of the model. There exists a variety of techniques for hyperparameter tuning Grech et al. (2008) ; Hansen (2010) , the most principled of which are generalized cross-validation (GCV) Arlot and Celisse (2010) and parametric Empirical Bayesian inference. The Empirical Bayesian approach is to place a prior probability distribution on the hyperparameters and to subsequently estimate the hyperparameters by applying Bayes' rule. It can be straightforwardly applied to linear inverse operators Mitsuhata (2004) . It is unclear, however, if the MEG/EEG sensor-data contain sufficient information about the hyperparameters to update the prior distribution, i.e. if the hyperparameters are identifiable. This is unknown even for the single hyperparameter p in the original formulation of Harmony, whose value was set by hand Petrov (2012) . Identifiability can be investigated theoretically using appropriate statistical techniques such as Fisher information or empirically by using simulated of real MEG/EEG data. Also unknown is which tuning method is most suitable in a given situation and to what extent it increases the variability of the source-reconstructions. Interesting directions for future research are to conduct a comparative evaluation of the tuning properties of different inverse operators and their effects on reconstruction quality and to compare different tuning methods.
A possible generalization of the proposed framework is the use of different penalization functions in the optimization function ( Eq. (5) ), most notably L 1 -penalization to encourage sparseness of the source reconstructions in the basis of spherical Harmonics Grech et al. (2008) ; Vega-Hernandez et al. (2008) . In contrast to L 2 -regularization, which yields source reconstructions with continuous angular spectra, sparse angular spectra correspond to cortical activity with a small number of well-separated wavelengths and it is currently unclear what type of cortical activity possesses such features. The combination of L 1 -and L 2regularization functions (i.e. elastic nets) might be promising for modeling certain types of cortical activity. Although the ensuing optimization problems are non-convex, reliable iterative algorithms are available to obtain locally optimal solutions Vega-Hernandez et al. (2008) ; Vega-hernández et al. (2019) . We make two remarks about the use of elastic nets in this context. First, more flexible regularization functions require more hyperparameters and it is currently unknown if these pa-rameters can actually be estimated from MEG/EEG sensor data. Second, the angular resolution analysis as presented in Section 2.6 has to be modified, since the source-reconstructions obtained from invoking L 1 -penalization functions are non-linear in the sensor data. Another promising direction concerns recent developments in non-cryogenic detection of magnetic fields, based on optically pumped magnetometers (OPM's) Baillet (2017) . Because OPM sensors are substantially closer to the skull than traditional MEG coils, the spatial resolution is potentially higher, given that a sufficient number of sensors is used. This implies that small-wavelength cortical activity can be more accurately reconstructed, which is particularly relevant for higher-frequency oscillations in the beta and gamma frequency bands and local epileptic discharges.