Beyond deep versus superficial: true laminar inference with MEG

Multilayer MEG source modelling enables depth-resolved inference across all six cortical laminae in humans | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Multilayer MEG source modelling enables depth-resolved inference across all six cortical laminae in humans View ORCID Profile Maciej J. Szul , Ishita Agarwal , View ORCID Profile Quentin Moreau , View ORCID Profile Solène Gailhard , View ORCID Profile Carolina Fernandez Pujol , View ORCID Profile Yunkai Zhu , Matteo Maspoli , View ORCID Profile Danila Shelepenkov , View ORCID Profile Bassem Hiba , View ORCID Profile Sebastien Daligault , View ORCID Profile Franck Lamberton , View ORCID Profile Denis Schwartz , View ORCID Profile Mathilde Bonnefond , View ORCID Profile Andrew R. Dykstra , View ORCID Profile Sven Bestmann , View ORCID Profile Gareth R. Barnes , View ORCID Profile James J. Bonaiuto doi: https://doi.org/10.1101/2025.05.28.656642 Maciej J. Szul 1 Institut des Sciences Cognitives Marc Jeannerod , CNRS UMR5229, Lyon, France 2 Université Claude Bernard Lyon 1, Université de Lyon , France 3 Department of Psychiatry and Psychotherapy, Faculty of Medicine, University of Tübingen , Germany Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Maciej J. Szul Ishita Agarwal 1 Institut des Sciences Cognitives Marc Jeannerod , CNRS UMR5229, Lyon, France 2 Université Claude Bernard Lyon 1, Université de Lyon , France 4 Department of Psychological Sciences, Purdue University , West Lafayette, Indiana, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Quentin Moreau 1 Institut des Sciences Cognitives Marc Jeannerod , CNRS UMR5229, Lyon, France 2 Université Claude Bernard Lyon 1, Université de Lyon , France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Quentin Moreau Solène Gailhard 1 Institut des Sciences Cognitives Marc Jeannerod , CNRS UMR5229, Lyon, France 2 Université Claude Bernard Lyon 1, Université de Lyon , France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Solène Gailhard Carolina Fernandez Pujol 5 Department of Neuroscience, Carney Institute of Brain Science, Brown University , Providence, RI, USA 6 Department of Biomedical Engineering, University of Miami , Coral Gables, FL, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Carolina Fernandez Pujol Yunkai Zhu 6 Department of Biomedical Engineering, University of Miami , Coral Gables, FL, USA 7 School of Communication Sciences and Disorders, University of Central Florida , Orlando, FL, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Yunkai Zhu Matteo Maspoli 1 Institut des Sciences Cognitives Marc Jeannerod , CNRS UMR5229, Lyon, France 2 Université Claude Bernard Lyon 1, Université de Lyon , France 8 Lyon Neuroscience Research Center , CRNL, INSERM U1028, CNRS UMR5292, Lyon, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Danila Shelepenkov 2 Université Claude Bernard Lyon 1, Université de Lyon , France 8 Lyon Neuroscience Research Center , CRNL, INSERM U1028, CNRS UMR5292, Lyon, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Danila Shelepenkov Bassem Hiba 1 Institut des Sciences Cognitives Marc Jeannerod , CNRS UMR5229, Lyon, France 2 Université Claude Bernard Lyon 1, Université de Lyon , France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Bassem Hiba Sebastien Daligault 9 CERMEP-Imagerie du Vivant , Lyon, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Sebastien Daligault Franck Lamberton 9 CERMEP-Imagerie du Vivant , Lyon, France 10 Lyon-East Health Research Federation, CNRS UAR3453, INSERM US7, Université Claude Bernard Lyon 1 , Lyon, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Franck Lamberton Denis Schwartz 8 Lyon Neuroscience Research Center , CRNL, INSERM U1028, CNRS UMR5292, Lyon, France 9 CERMEP-Imagerie du Vivant , Lyon, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Denis Schwartz Mathilde Bonnefond 2 Université Claude Bernard Lyon 1, Université de Lyon , France 8 Lyon Neuroscience Research Center , CRNL, INSERM U1028, CNRS UMR5292, Lyon, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Mathilde Bonnefond Andrew R. Dykstra 6 Department of Biomedical Engineering, University of Miami , Coral Gables, FL, USA 7 School of Communication Sciences and Disorders, University of Central Florida , Orlando, FL, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Andrew R. Dykstra Sven Bestmann 11 Department of Imaging Neuroscience, UCL Queen Square Institute of Neurology, University College London , London, UK 12 Department of Clinical and Movement Neurosciences, UCL Queen Square Institute of Neurology, University College London , London, UK Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Sven Bestmann Gareth R. Barnes 11 Department of Imaging Neuroscience, UCL Queen Square Institute of Neurology, University College London , London, UK Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Gareth R. Barnes James J. Bonaiuto 1 Institut des Sciences Cognitives Marc Jeannerod , CNRS UMR5229, Lyon, France 2 Université Claude Bernard Lyon 1, Université de Lyon , France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for James J. Bonaiuto For correspondence: james.bonaiuto{at}isc.cnrs.fr Abstract Full Text Info/History Metrics Preview PDF Abstract Neural dynamics at the laminar level are critical for cortical computation. However, in humans, non-invasive methods to probe such dynamics have been limited to coarse distinctions between deep and superficial layers. Here, we demonstrate that under certain conditions, magnetoencephalography (MEG) can achieve laminar inference by localizing sources at the level of individual cortical laminae. Using a multilayer source reconstruction approach, we systematically assess the limits of MEG depth resolution, and show that laminar precision is achievable under realistic signal-to-noise ratios and co-registration accuracy. We show that accurate laminar inference depends critically on aligning forward model dipole orientations with true cortical column orientations, and that regional variations in cortical anatomy influence reconstruction fidelity. We then apply this approach to empirical data from three independent datasets, revealing the expected laminar patterns of activity during event-related fields in primary visual, somatosensory, and motor cortices. These findings position MEG as a powerful tool for investigating lamina-specific neural dynamics in cognition and behavior, and offer new opportunities to bridge invasive electrophysiology and human neuroimaging. Teaser Multilayer MEG enables depth-resolved localization across cortical laminae under achievable conditions. Introduction The morphology, connectivity, and computations of cortical laminae have been central topics in neuroscience for over a century 1 . Invasive laminar electrophysiology has provided fundamental insights into sensory and sensorimotor computations - e.g., laminar specific receptive fields 2 and the directional flow of information in feed-forward and feedback pathways 3 - that would be otherwise inaccessible. In cognitive neuroscience, canonical microcircuit models and other laminar interaction frameworks have generated mechanistic hypotheses about lamina-specific contributions to perception, cognition, and conscious processing 4 – 9 . Clinically, assessment of laminar dynamics has improved our understanding of epileptogenic propagation 10 and holds promise for refining our knowledge of disorders such as Parkinson’s disease 11 , but these investigations have been mostly limited to invasive methods. Invasive electrophysiological recordings are widely used in animal models, whereas they are restricted to specific clinical populations in humans 12 . Traditional invasive recordings provide direct access to neuronal events, capturing single-unit and local population activity with high spatial and temporal resolution, though newer technologies such as Neuropixels now extend this access across entire cortical columns and brain structures 13 . However, their use requires surgical implantation, necessitating ethical considerations in animal research, and medical justification in human studies. Furthermore, although invasive recordings enable laminar-level recordings in cortex, they remain extracellular and thus are still subject to an inverse problem 14 . This limitation has motivated recent efforts to develop magnetrodes, which aim to directly measure intracellular currents 15 . Nevertheless, these invasive methods are not feasible for broader applications in cognitive neuroscience. In contrast, for both practical and ethical considerations, non-invasive methods like electroencephalography (EEG), magnetoencephalography (MEG), and functional magnetic resonance imaging (fMRI) dominate human neuroscience 16 , 17 . Their whole brain coverage enables the study of network-level processing, but each modality has inherent limitations. EEG and MEG provide high temporal resolution, but suffer from an inverse problem and are thought to lack the spatial specificity needed for laminar inference 17 – 19 , whereas fMRI has excellent spatial, but poor temporal resolution 16 , 20 . Recent advances in high resolution fMRI have enabled laminar-level analysis, but its inherently limited temporal resolution precludes the study of laminar dynamics 21 – 23 . MEG has traditionally been considered incapable of resolving sources with a precision at or below the thickness of the cortex 18 , 24 . However, this assumption has been challenged by studies indicating that there is no fundamental theoretical barrier to achieving high spatial resolution with MEG 25 , 26 ; the primary limitation is the practical issue of signal-to-noise ratio (SNR). With sufficiently high SNR, MEG can, in principle, resolve sources with much finer precision than traditionally assumed 25 , 26 . Achieving high signal quality in actual recordings requires minimizing co-registration errors, between-session head position variability, and within-session movement, all of which can be addressed using customized head-casts or 3D-printed helmets that stabilize the head relative to the sensor array, an approach called high precision MEG (hpMEG) 27 – 29 . hpMEG spurred the development of early laminar MEG approaches, which involved Bayesian model comparison between separate forward models for two cortical depths: deep and superficial, or analysis of source signals using a forward model combining surfaces for the two layers 25 , 26 , 30 . These approaches have been used to successfully detect laminar differences predicted by theoretical models and invasive recordings, such as lamina-specific spectral activity 30 and distinct laminar signatures of beta bursts 31 . However, previous studies have been limited to coarse deep versus superficial distinctions. Here, we demonstrate that hpMEG can resolve sources across individual cortical laminae. We expand the source space to include an arbitrary number of intermediate layers. We refer to the equidistant depth surfaces used in our forward models as layers (denoted with Arabic numerals, i.e. 1-11), while we use laminae (denoted with Roman numerals) to refer to cytoarchitectonically defined cortical layers (i.e. I-VI). We first use simulations to systematically probe the limits of data quality required for true laminar precision. We show that accurate laminar inference critically depends on aligning the dipole orientations in the forward model with the true cortical column orientations, and we demonstrate that variations in anatomical features across the cortex leads to regionally varying reconstruction accuracy, consistent with prior findings 26 . Finally, we demonstrate the method’s applicability to real data by reconstructing laminar activity profiles from visual and motor event-related fields (ERFs). In the visual cortex, hpMEG recovers the expected feedforward activation sequence consistent with laminar models of sensory processing. In the somatosensory and motor cortices, analysis of two independent datasets reveals reproducible deep-superficial activation patterns characteristic of motor execution. These results position hpMEG as a powerful tool for interrogating neural circuits at a laminar level, offering new opportunities to study how computations within and across cortical laminae shape network-level interactions underlying cognition and behavior 7 . Methods Empirical Data Acquisition We validated multilaminar MEG inference using three independent human datasets comprising visual and motor paradigms acquired with high-precision MEG. Across datasets, healthy adult participants underwent structural MRI for individual cortical surface reconstruction and MEG recordings acquired with whole-head CTF systems at high sampling rates, with head position stabilized using participant-specific head-casts and continuously monitored using fiducial coils. The datasets span visually evoked responses in primary visual cortex, as well as self-paced and visually cued button presses engaging precentral and postcentral cortices, thereby sampling distinct cortical regions, task demands, and acquisition sites. Despite differences in experimental design, all datasets shared a common preprocessing and analysis framework, enabling consistent laminar source reconstruction and cross-dataset comparison of laminar activation patterns. Full details of participants, tasks, MRI and MEG acquisition protocols, preprocessing pipelines, and analysis parameters for each dataset are provided in the Supplementary Materials. Empirical Data Preprocessing For all three empirical datasets, individual cortical surfaces were reconstructed from structural T1- and T2-weighted MRI using FreeSurfer recon-all (v6.0.1; see Supplementary Materials for dataset-specific MRI protocols) 32 . We then used FreeSurfer’s mris_expand function to generate 9 equidistant intermediate surfaces between the pial boundary surface and white matter boundary surface, yielding 11 cortical surfaces in total ( Figure 1b ). An arbitrary number of surfaces could be used, but we chose 11 for computational efficiency. These intermediate surfaces are not meant to represent the 6 cortical laminae, but rather to sample the cortical depth like a laminar electrode with multiple, evenly spaced contacts. Freesurfer operates on hemispheres independently, resulting in deep vertices and mesh faces cutting through subcortical structures. These vertices and associated faces were removed from each hemisphere of each layer mesh. The left and right hemisphere meshes of each intermediate surface were then combined by concatenation of their vertices and faces (left then right). Download figure Open in new tab Figure 1. Overview of the simulation strategy and analysis. a) Pial and white matter boundaries surfaces are extracted from anatomical MRI volumes. b) Intermediate equidistant surfaces are generated between the pial and white matter surfaces (labeled as superficial (S) and deep (D) respectively). c) Surfaces are downsampled together, maintaining vertex correspondence across layers. Dipole orientations are constrained using vectors linking corresponding vertices (link vectors). d) The thickness of cortical laminae varies across the cortical depth 40 – 42 , which is evenly sampled by the equidistant source surface layers. e) Each colored line represents the model evidence (relative to the worst model, ΔF) over source layer models, for a signal simulated at a particular layer (the simulated layer is indicated by the line color). The source layer model with the maximal ΔF is indicated by “˄”. f) Result matrix summarizing ΔF across simulated source locations, with peak relative model evidence marked with “˄”. g) Error is calculated from the result matrix as the absolute distance in mm or layers from the simulated source (*) to the peak ΔF (˄). h) Bias is calculated as the relative position of the peak ΔF(˄) to the simulated source (*) in layers or mm. The resulting bihemispheric meshes each contained approximately 300,000 vertices, exceeding computational feasibility for source reconstruction. FreeSurfer generates the pial mesh by expanding the white matter boundary mesh until the pial boundary is reached. Thus, every pial mesh vertex has a corresponding white matter mesh vertex. The vectors linking these corresponding vertices, known as link vectors 33 , have been shown to provide the best estimate of cortical column orientation out of available methods 34 . One advantage of link vectors for constraining dipole orientation is that the vectors are equivalent across corresponding vertices in all intermediate surfaces. This ensures that model-based layer comparison is based solely on dipole location. However, most methods for mesh downsampling disrupt the correspondence between vertices when applied independently to each layer as they alter vertex coordinates rather than simply removing vertices. We therefore downsampled the pial boundary surface by a factor of 10, using the vtkDecimatePro function of the VTK library (v9.3.0) 35 because it is purely subtractive (i.e. it only removes vertices), resulting in a mesh with approximately 30,000 vertices. We then downsampled the other 10 meshes by removing the same vertices that were removed from the pial surface, thus preserving vertex correspondence, and copying the pial mesh face structure (i.e. edges between vertices) 34 . This enabled a straightforward computation of link vectors ( Figure 1c ), which were used to constrain dipole orientations in all of the main simulations and analyses. MEG data were preprocessed separately for each dataset, time-locked to the event of interest (visual grating onset or button press), and averaged over trials to obtain event-related fields used for laminar source reconstruction (details in Supplementary Materials). Simulations We tested the efficacy of each analysis method using synthetic datasets based on the cortical surfaces and MEG data of a single participant from the visually cued button-press dataset (see Supplementary Materials). The MEG data was only used to define the sensor layout, sampling rate (600 Hz), number of trials (531), and number of samples (1201) for the simulations; the MEG sensor data itself was discarded. All simulations and analyses were implemented using the laMEG software package (v0.0.7; https://github.com/danclab/laMEG ), built on a custom version of SPM 36 compiled as a python library ( https://github.com/danclab/spm ), and are available at http://github.com/danclab/multilaminar_sim . For each simulation, we specified a source centered at a vertex on one of the 11 cortical surfaces spanning the depth of the cortex from the pial to the white matter surface. We simulated a localized patch of current density with a Gaussian-shaped temporal profile over a 2 s time window with a peak at the center of the window and full width at half maximum (FWHM) of 400 ms (time course displayed in Figure S1 ). The spatial extent of the simulated source was set to 5 mm FWHM (corresponding to a mean patch size of 14.93 vertices). To generate synthetic MEG data, we used a single-shell forward model 37 based on a surface composed of the 11 cortical layer meshes concatenated into a single mesh. The dipole moment for each simulation was set to 10 nAm unless otherwise specified. We selected 100 random cortical locations (vertices on the pial surface) and conducted separate simulations for different sensor-level signal-to-noise ratio (SNR) levels (−500, −50, −35, −20, −10, −5, 0, and 5 dB) and co-registration error levels (0, 0.5, 1, 2, 3, 4, and 5 mm). White noise was added at the sensor level, scaled to achieve the desired per-trial SNR (computed as the ratio of signal power to noise power across all sensors) 26 . Download figure Open in new tab Figure S1. Effect of sliding window size on laminar inference accuracy. Laminar model evidence over time (ΔF, relative to worst model at each time point) for simulated dipolar sources across 11 cortical depths, shown for sliding window sizes of 15 ms, 25 ms, 50 ms, and 100 ms (rows). Columns correspond to sources simulated at superficial, middle, and deep cortical depths. Each trace reflects the average ΔF across 100 simulated cortical locations at a fixed sensor-level SNR of −20 dB. Colors indicate cortical evaluation depth (S = superficial surface layer, D = deep surface layer). The black dashed line shows the simulated signal time course. Inset panels in the first row (a-c) show model evidence at t = 0ms, with peak ΔF values aligning with the simulated source depth (dashed vertical line). The depth-specific pattern of model evidence was consistent across window sizes for each simulated depth (columns). Results confirm that laminar model discrimination is consistent across a wide range of window sizes, with all tested windows successfully recovering the correct source depth. To simulate MEG-MRI co-registration error, we applied a random rigid-body perturbation to the MRI-defined nasion (NAS), left preauricular (LPA), and right preauricular (RPA) fiducials. Let f i denote the three fiducial coordinates and their barycenter; we first centered the fiducials by , then drew (i) a random translation t with isotropic direction and normally distributed signed magnitude with scale equal to the target error level (mm), and (ii) a random rotation R about an isotropically distributed axis with rotation-vector magnitude drawn from a zero-mean normal distribution scaled by the target error level. We applied the transform about the barycenter, , which preserves inter-fiducial distances (i.e., a rigid triangle) but yields different displacements at each fiducial due to rotation. Because the combined translation+rotation does not, in general, make every fiducial move by the same amount, we operationalized the “error level” as the maximum fiducial displacement, , and used rejection sampling to redraw ( t , R ) until | d max - d target | < 0.05 mm. The resulting perturbed fiducials ( NAS′ , LPA′ , RPA′ ) were then used for MEG-MRI co-registration, thereby simulating misalignment between the subject’s head position in the MEG scanner and their anatomical MRI 26 . To evaluate the accuracy of laminar source reconstruction in the presence of multiple simultaneous sources, we conducted a separate set of simulations in which two sources of equal strength were generated at different cortical depths. For each of the 100 selected cortical locations, we fixed one source in the middle layer while systematically varying the second source from the most superficial to the deepest cortical surface. The simulated sources followed the same Gaussian-shaped temporal profile as in the single-source simulations, with a dipole moment of 10 nAm and a spatial extent of 5 mm FWHM. Synthetic MEG data were generated using a single-shell forward model, with white noise added to achieve a sensor-level SNR of −20 dB. To evaluate the impact of different cortical column orientation estimation methods on laminar source reconstruction accuracy, we systematically varied the approach used to define dipole orientations during both simulation and reconstruction. We tested seven distinct methods, including surface normal vectors computed from the original high-resolution mesh 34 and a downsampled surface mesh 38 , cortical patch statistics 39 , and link vectors 33 , with additional conditions for whether dipole orientations were fixed across layers. For each method, we simulated sources at the 100 selected cortical locations, assigning simulated dipole orientations based on the selected approach and then reconstructing sources using each of the different orientation definitions. The downsampled surface normal method estimated dipole orientation at each vertex of the decimated cortical mesh by averaging the normal vectors of all adjacent triangular faces. This was implemented using the spm_mesh_normals function in SPM 36 . The original surface normal method followed the same approach but was applied to the full-resolution cortical surface before being mapped onto the decimated mesh for source localization 34 . The cortical patch statistics method calculated the mean normal vector at each vertex of the downsampled mesh using all neighboring vertices from the original high-resolution surface 39 . In contrast, the link vector method defined dipole orientation as the vector linking each pial surface vertex to its corresponding white matter vertex 33 , leveraging the one-to-one correspondence maintained during mesh decimation. Because different orientation estimation methods can yield subtly different vector fields, we assessed the angular discrepancies between these definitions and their effects on laminar inference accuracy. Dipole orientation vectors derived from each method were incorporated into the lead field matrix used for source inversion 34 in order to evaluate how orientation differences impact model evidence and depth localization accuracy across simulations. For each of these methods, we tested versions where the orientation was computed independently for each layer and versions where it was computed on the pial surface and fixed across layers (except for the link vectors method which is equivalent across layers by definition). To examine how anatomical features influence the accuracy of laminar source reconstruction, we generated a large-scale set of simulations in which sources were placed at every vertex of a layered cortical mesh (one source per simulation). The mesh consisted of 29,130 locations at 11 cortical depths, resulting in a total of 320,430 simulated sources. Each simulation used a Gaussian-shaped temporal profile with a full width at half maximum (FWHM) of 400 ms and a dipole moment of 10 nAm ( Figure S1 includes the time course of this simulated signal). Synthetic MEG data were generated using a single-shell forward model, and white noise was added to yield a sensor-level SNR of −20 dB. All of these simulations were conducted with zero co-registration error. Simulation Analysis For each simulated dataset, we applied multiple sparse priors (MSP) inversion 43 to each simulated dataset within a Hann windowed 50 ms analysis window centered around the peak of the simulated Gaussian temporal profile. These data were projected into 274 orthogonal spatial modes and 4 temporal modes, with a patch size of 5 mm (unless otherwise specified) and a single prior corresponding to the vertex the activity was simulated at. This was done to test the feasibility of a sliding time window approach previously applied to bilaminar MEG dynamics 31 . Because each time window is evaluated independently in this framework, a single 50 ms window was sufficient for our analysis. We confirmed the robustness of this window size in supplementary simulations testing 15, 25, 50, and 100 ms windows ( Figure S1 ), which showed consistent reconstruction performance across a broad temporal range. Rather than performing source reconstruction on a single multilayer cortical mesh, we conducted separate MSP inversions for each of the 11 cortical depth surfaces, treating each as an independent forward model. For each simulation, we computed the model evidence (free energy) for each depth-specific source reconstruction and identified the cortical surface with the highest model evidence ( Figure 1e ). This allowed us to infer the most likely depth of the simulated source by comparing the relative free energy across reconstructions at different cortical depths ( Figure 1f ). We quantified reconstruction accuracy by identifying the surface with the highest model evidence and comparing its depth to the known simulated source depth. Reconstruction error was computed as the absolute difference between the simulated and inferred depths ( Figure 1g ), measured both in number of surfaces and in millimeters by scaling surface indices using local cortical thickness. Reconstruction bias was defined as the signed difference between simulated and inferred depths ( Figure 1h ). To assess statistical significance, we performed permutation tests using 10,000 iterations, shuffling the free energy matrices for each simulated cortical depth and recomputing error and bias distributions under the null hypothesis of no depth sensitivity. P-values were calculated by comparing observed reconstruction errors and biases against these shuffled distributions, with statistical significance determined at p < 0.05. Effect sizes were computed using Cohen’s d, defined as the difference between observed and shuffled means divided by the standard deviation of the shuffled distribution. We conducted bootstrap resampling (5,000 iterations), drawing random subsets of the data to compute 95% confidence intervals for error and bias. We also computed a diagonal dominance score to quantify the extent to which model evidence matrices favored reconstructions at the correct depth. This score, representing the ratio of diagonal values to total values in the free energy matrix, was compared against a permuted distribution ( N = 10,000) to assess whether observed diagonality exceeded chance levels. We accounted for cortical laminar thickness variability by incorporating the BigBrain atlas 44 , 45 . Because cortical laminae differ in thickness across the brain ( Figure 1d ), mapping simulated and inferred sources to cytoarchitectonic laminae rather than fixed-depth surfaces improves biological interpretability. To achieve this, we mapped the BigBrain atlas to the subject’s anatomy via FreeSurfer’s surface-based registration through fsaverage space 46 . This allowed us to align the proportional cortical laminae boundaries from the atlas with the subject’s individual cortical geometry. For each simulated cortical location, we extracted the proportional laminae boundaries (six depth values defining the transitions between laminae) using precomputed laminar thickness estimates from the atlas. These boundaries were then scaled by the local cortical thickness to derive laminar depths in millimeters. The final mapping step assigned each of the 11 cortical surface layers used for source reconstruction to a corresponding cytoarchitectonic laminae. Reconstruction errors were then computed in laminar space rather than in fixed surface coordinates; specifically, an inferred source depth was only considered an error if it fell outside the lamina of the simulated source, regardless of its precise position in cortical depth. Significance was evaluated using the same permutation-based approach used for error in terms of layers and millimeters ( N = 10,000). To evaluate how local anatomy influences laminar inference, we quantified reconstruction error at each cortical location and related it to anatomical features. For each vertex, we computed absolute reconstruction error in millimeters. We then analyzed the relationship between reconstruction error and several anatomical properties: cortical thickness (measured between the pial and white matter surfaces), distance to the scalp surface, orientation of the cortical column relative to the scalp, and variability in lead field magnitude across depth. Orientation was computed as the absolute dot product between the simulated dipole orientation and the scalp normal at each vertex. Lead field variability was quantified as the root mean square difference in lead field magnitude across surface layers, relative to the superficial (pial) surface. These features were then related to reconstruction error using kernel density estimation and receiver operating characteristic (ROC) analysis to identify anatomical thresholds predictive of accurate inference. To assess whether anatomical features jointly predict laminar inference accuracy better than single features, we performed multivariate logistic regression and random forest classification. We evaluated model performance using 5-fold cross-validated ROC analysis. Logistic regression models were tested both with and without pairwise interaction terms to quantify potential nonlinear feature interactions. Random forest hyperparameters were set to 500 trees, unlimited depth, and balanced class weights. Empirical Data Analysis Laminar source reconstruction used the laMEG software package (v0.0.7; https://github.com/danclab/laMEG ), built on a custom version of SPM 36 compiled as a python library ( https://github.com/danclab/spm ). For each participant and session, candidate sources were first localized by co-registering the averaged epochs to the combined 11-layer mesh (formed by concatenating all depth surfaces into a single source space), and then using an empirical Bayes beamformer (EBB) with a 5-mm patch size, four temporal modes, and a number of spatial modes equal to the rank of the preprocessed MEG data, within a task- and component-specific time window (see Supplementary Materials for dataset- and component-specific details). This EBB step was used solely for spatial localization and vertex selection; all laminar inference and model comparison were performed using multiple sparse priors (MSP), as described below. To select a single vertex for laminar inference within each region of interest (primary visual, motor, or somatosensory cortex depending on the ERF and component), we combined functional and anatomical criteria. Source time series were extracted from the middle cortical surface, and vertices were ranked by mean signal magnitude within the component-specific analysis window relative to baseline (defined per dataset and component; see Supplementary Materials). Among vertices in the anatomically defined region of interest, we retained the top 1% by signal magnitude and selected the vertex with the highest composite anatomical suitability score. This score was defined as the sum of z-scored cortical thickness, deep-superficial lead-field RMS difference, alignment between cortical column orientation and local scalp normal, and inverse distance to the scalp, thereby favoring locations with both strong signal and favorable anatomical properties for laminar discrimination. At the selected vertex, sliding-window Bayesian model comparison was then performed across all 11 cortical depths. We applied multiple sparse priors (MSP) inversion to each dataset within overlapping 25 ms windows (unless otherwise specified), with Hann windowing applied within each window, a number of spatial modes equal to the rank of the preprocessed data, 4 temporal modes, a patch size of 5 mm (unless otherwise specified), and a single prior corresponding to the selected vertex. Separate sliding window MSP inversions were performed for each of the 11 cortical depth surfaces, treating each surface as an independent forward model. This yielded model evidence (free energy) time series for each depth-specific source reconstruction, for each subject and session. In order to interpret the cortical depth surfaces in terms of laminae, we mapped the BigBrain atlas to each subject’s anatomy 46 , computed proportional laminae boundaries from the atlas, scaled these boundaries by the local cortical thickness, and assigned each of the 11 cortical surface layers to a corresponding cytoarchitectonic lamina. Free energy values were then averaged within laminae. The resulting model evidence was then compared in two ways: i) by computing model evidence relative to the worst model (ΔF) at each time point, and ii) using Bayesian model selection to identify the most probable family of models at the group level at each time point. For Bayesian model selection we used the VBA (Variational Bayesian Analysis) toolbox 47 . The free energy measures of the six laminae were grouped in three families: supragranular (laminae I-III), granular (lamina IV) and infragranular (laminae V-VI) families. A group-level random effect analysis (with participant as a random effect) was used to compute protected exceedance probabilities (PEPs) 48 at each time point, quantifying the likelihood of one family model being more probable than the others. Results SNR-dependent patterns of laminar inference accuracy To determine the SNR required for accurate laminar reconstruction across multiple cortical depths, we selected 100 random cortical locations (vertices on the pial surface). At each location, we ran 11 simulations, each simulating a source signal with a Gaussian temporal profile at the corresponding vertex on one of 11 surface meshes ranging from the white matter (deepest) to pial (most superficial) surface (see Figure S1 for simulated time-course). For each simulation, we ran source reconstruction (MSP, with the prior set to the simulated cortical location) within a 50 ms central window 11 times, each time using a different surface mesh to construct the forward model, and compared the resulting model evidence, approximated by free energy. In the absence of noise (e.g. sensor noise, co-registration error etc), the model based on the surface corresponding to the depth of the simulated signal will have greater model evidence, or free energy (i.e. if a signal is simulated on the pial surface, the model based on the pial surface should have the greatest free energy). This set of simulations was repeated at each selected cortical location with different levels of simulated sensor-level SNR, varying from −∞ (pure noise) to 5 dB. For each location, this yielded a matrix of model evidence (free energy) with a column for simulation at each cortical depth, and a row for each forward model. For each simulated cortical depth (matrix column), we transformed free energy values for each forward model by subtracting the lowest value, thus representing model evidence relative to the worst model. For each matrix, we identified the model with the peak free energy for each simulated cortical depth, and generated a probability density of these peaks for each SNR level. We also computed a measure of diagonal dominance, or the ratio of diagonal values to total values (perfect reconstruction would yield a purely diagonal matrix), and compared this to a null distribution obtained by shuffling the matrices. At low levels of SNR (< −50 dB), there was no pattern in the mean free energy (averaged over all simulated locations) over simulated cortical depths ( Figure 2a ). However, the probability density of surfaces with the maximal free energy demonstrates that in most locations, the peak free energy corresponded to either the most superficial or the deepest surface model ( Figure 2e ). At −50 dB, a step change in the mean free energy was observed, with the greatest mean free energy for roughly the surface corresponding to the simulated depth ( Figure 2b ), and a diagonal dominance score greater than expected by chance ( d = 3.22, p = 0.001; Figure S2a ). A corresponding gradient in the probability density of the peak free energy model was observed, with a greater probability for superficial surface model peaks for superficial simulations, and for deep surface model peaks for deep simulations ( Figure 2f ). At higher levels of SNR (≥- 35 dB), both the mean free energy ( Figure 2c,d ) and the probability density ( Figure 2g,h ) exhibited significant diagonality (all d > 3.5, all p < 0.001; Figure S2a ), indicating near perfect inference of the depth of the simulated signal. Download figure Open in new tab Figure S2. Diagonal-dominance scores for laminar model-evidence matrices across different SNR and co-registration error values. (a) Varying sensor-level SNR, and (b) varying co-registration error. Blue curves show observed diagonal dominance (trace / sum of each matrix), and red curves show the mean from permuted matrices with shaded 95% CI. Black stars mark SNR or co-registration levels where observed diagonal dominance differed significantly from chance (p < 0.05). Download figure Open in new tab Figure 2. Laminar source-reconstruction accuracy as a function of sensor-level SNR. a-d) Mean model evidence (free energy, ΔF) matrices (row-normalized to the worst model; S = superficial surface layer, D = deep surface layer) for simulated signals placed on 11 cortical surfaces at four SNR levels (−∞ dB, −50 dB, −35 dB, −20 dB). Accurate laminar inference begins to emerge around −35 dB ( c ), with near-perfect identification at −20 dB ( d ). Triangular markers denote the forward model with peak free energy for each simulated depth. e-h) Probability density of each forward model being the maximum-evidence solution, aggregated across cortical locations: a diagonal pattern appears from −35 dB, indicating accurate depth localization. i, j) Laminar reconstruction error (in layers and mm) falls below chance at −35 dB and improves at higher SNR; shuffled controls show ∼1 mm error due to chance matches between simulated and evaluated layers. k, l) Laminar reconstruction bias (in surface layers and millimeters), again comparing simulated and shuffled data. Shaded bands represent 95% confidence intervals, and asterisks mark SNR levels at which observed error or bias differed significantly from chance (p < 0.05). At each cortical location, we computed bias as the mean distance, over the 11 simulations, between the simulated source and the surface model with the highest free energy (in terms of surfaces and millimeters), and error as the mean absolute distance between the simulated surface and the best surface models (also in surfaces and millimeters). For each metric, we compared the observed values across simulated SNR levels with the null distribution obtained by shuffling the surface model free energy for each simulated depth. At low SNR levels (−50 dB and below), both error and bias were at chance levels ( p > 0.05 and p > 0.8, respectively). A significant reduction in error emerged at −35 dB ( d = −4.17, p 4.8, all p 0.8), though extreme deep or superficial biases were observed at some locations below −50 dB ( Figure 2k,l ). These results were nearly identical when error and bias were computed in millimeters rather than surface layers ( Figure 2k,l ). Consistent with previous simulation results 26 , inference was biased at low SNR levels when using the empirical Bayesian beamformer algorithm (EBB; Figure S3 ) 49 rather than MSP with a spatial prior, as used in this study. We also examined error separately for each simulated surface layer and found no differences between layers ( Figure S4a ). Laminar inference at a coarse scale (deep versus superficial) is therefore possible in principle at −50 to −35 dB SNR, with fine-scale inference possible at −20 dB and above. Download figure Open in new tab Figure S3. Laminar source-reconstruction accuracy with EBB as a function of SNR. a-d) Mean model evidence (free energy, ΔF) matrices (row-normalized to the worst model; S = superficial surface layer, D = deep surface layer) for simulated signals placed on 11 cortical surfaces at four SNR levels (−∞ dB, −50 dB, −35 dB, −20 dB). Triangular markers denote the forward model with peak free energy for each simulated depth. e-h) Probability density of each forward model being the maximum-evidence solution, aggregated across cortical locations. i, j) Laminar reconstruction error (in layers and mm). k, l) Laminar reconstruction bias (in surface layers and millimeters), again comparing simulated and shuffled data. Shaded bands represent 95% confidence intervals, and asterisks mark SNR levels at which observed error or bias differed significantly from chance (p < 0.05). Download figure Open in new tab Figure S4. Layer-by-layer source-reconstruction error. (a) Mean error as a function of SNR for each simulation surface layer, from superficial (S, light cyan) to deep (D, magenta). (b) Corresponding layer-by-layer error over increasing co-registration error. Impact of co-registration error on laminar inference Previous work suggests that accurate laminar inference is feasible with up to 2 mm of co-registration error 25 , 26 . To evaluate whether this holds in our multilayer framework, we ran another set of similar simulations at the same cortical locations, this time fixing the SNR at 0 dB, and varying the simulated co-registration error from 0 to 5 mm. At levels of co-registration less than or equal to 2 mm, the surface model with the highest mean free energy corresponded to the simulated cortical depth ( Figure 3a-c ), and the peak free energy model probability density exhibited strong diagonality (all d > 4, all p < 0.001: Figure 3e-f ; Figure S2b ), indicating accurate inference across the range of simulated depths. At higher levels of co-registration error (3-5 mm), the mean free energy still exhibited significant diagonality (all d > 3.9, all p < 0.001; Figure 3d ; Figure S2b ), but the probability density only indicated a gradient in free energy peaks in the deepest and most superficial surfaces ( Figure 3h ). Download figure Open in new tab Figure 3. Laminar source-reconstruction accuracy as a function of co-registration error at fixed 0 dB SNR. a-d) Mean free energy (ΔF) matrices (row-normalized to the worst model; S = superficial surface layer, D = deep surface layer) for simulated signals placed on 11 cortical surfaces at 0, 1, 2, and 5 mm registration error. Triangles mark the best-fitting forward model for each simulated depth. e-h) Probability density of each surface model being the maximum-evidence solution, aggregated across cortical locations. i, j) Laminar reconstruction error (in surface layers and millimeters) for simulated versus shuffled data at varying registration errors. k, l) Corresponding bias (in surface layers and millimeters). Shaded bands denote 95% confidence intervals, and asterisks mark significance at p < 0.05. We computed bias and error for the co-registration error simulations using the same approach as in the SNR simulations and compared the observed values against the null distribution obtained by shuffling. Error was lowest at 0 mm ( d = −4.84, p < 0.001) and progressively increased with greater co-registration error, though it remained significantly lower than chance levels up to 2 mm ( d = −2.77, p = 0.002; Figure 3i,j ). At 3 mm and above, error was at chance level (all p > 0.06). There was no significant bias at any level of co-registration error (all p > 0.7). These results were nearly identical when error and bias were computed in millimeters rather than surface layers ( Figure 3k,l ). When split by simulated surface layer, the accuracy of laminar inference for simulated deep sources was less affected by co-registration error than superficial sources ( Figure S4b ). Accurate laminar inference therefore requires a co-registration error of 2 mm or less, though a coarse deep-versus-superficial distinction may be possible in some locations at errors up to 5 mm, provided there is sufficient SNR. SNR and co-registration error required for true laminar inference Given that cortical laminae vary in thickness both within and across regions 40 – 42 , we assessed laminar reconstruction accuracy in terms of cytoarchitectonic laminae rather than fixed-depth surfaces. Specifically, surface layers that mapped to the same lamina were treated as equivalent (e.g., errors between adjacent surface layers were not considered errors if both corresponded to lamina VI). To achieve this, we re-analyzed the effects of SNR and co-registration error using the BigBrain atlas 44 , 45 , which we mapped to the subject’s anatomy via FreeSurfer’s surface-based registration through fsaverage space 46 . This allowed us to estimate proportional laminar thickness at each simulated cortical location ( Figure 4a-f ), assign surface layers to their corresponding laminae, and evaluate reconstruction accuracy in terms of cytoarchitectonic rather than geometric depth. Download figure Open in new tab Figure 4. Laminar source-reconstruction accuracy as a function of SNR, co-registration error, and cortical thickness. a-f) Cortical surfaces showing laminar thickness across the brain, derived from the BigBrain atlas and mapped onto the surface models (representing laminae I-VI respectively) used in the simulations. These maps illustrate the variability of laminar thickness across cortical regions, which influences inference accuracy. g) Mean laminar reconstruction error (in laminae) across SNR levels, computed as the absolute difference between the simulated and inferred lamina. Lines indicate the mean error for each lamina, with shaded bands representing 95% confidence intervals. Lighter shaded bands represent shuffled null distributions. Asterisks mark SNR levels where reconstruction error was significantly lower than chance. h) Mean laminar reconstruction error as a function of co-registration error. i) Relationship between laminar thickness and reconstruction error. Each line represents a different lamina (I-VI), showing how increased laminar thickness is associated with reduced reconstruction error. We re-analyzed both SNR and co-registration error simulations under this framework. Error was computed as the absolute difference between the simulated and inferred lamina, and significance was assessed using permutation tests ( N = 10,000) against a shuffled null distribution. At very low SNR levels (−100 dB and below), reconstruction error was indistinguishable from chance for all laminae (all p > 0.8). At −50 dB, inference remained at chance for middle laminae (laminae III and IV; both p > 0.07) but was significantly better than chance for superficial and deep laminae (all p < 0.001). By −35 dB, reconstruction accuracy improved across all laminae (all p < 0.001), with effect sizes ranging from d = −11.36 to d = −23.15 ( Figure 4g ). The most pronounced improvements were observed at −20 dB and above, where errors reached minimal levels (all d < −14.57, all p < 0.001). These results confirm that true laminar inference is possible at sufficiently high SNR, even when accounting for regional variability in cortical laminae thickness. Co-registration errors had a similar impact, with small misalignment (≤ 2 mm) having minimal effects on reconstruction accuracy across all laminae (all p < 0.001), while increasing error progressively degraded performance. At 2 mm, reconstruction accuracy remained high (all d < −4.45, all p 0.4; Figure 4h ). At 5 mm, only the most superficial (lamina I) and deepest laminae (laminae V and VI) remained significantly different from chance (lamina I: d = −8.51, p < 0.001; lamina V: d = −2.10, p = 0.018; lamina VI: d = −8.45, p < 0.001). These findings suggest that fine-scale laminar inference is robust to small co-registration errors but becomes unreliable at misalignment exceeding ∼3 mm, particularly for middle-layer sources. At errors of 5 mm, only deep and superficial sources retain some discriminability, indicating that accurate reconstruction at this scale requires co-registration precision within at least approximately 2 mm. Finally, we assessed how laminar reconstruction accuracy depends on laminar thickness, focusing on the results of simulations with SNR of at least −20 dB and co-registration errors less than or equal to 2 mm. Across simulated locations, thicker laminae exhibited systematically lower reconstruction error ( Figure 4i ), reflecting a greater spatial margin for accurate depth localization. This effect was most pronounced in deep and superficial laminae, where increased cortical thickness improved inference reliability. Middle laminae, in contrast, showed greater variability in error, consistent with the higher misclassification rates observed at moderate SNR and co-registration error levels. Cortical laminae thickness therefore impacts laminar reconstruction accuracy, with thicker laminae providing a greater margin for accurate depth localization. The increased variability in middle laminae suggests that reconstruction accuracy is not uniform across the cortex, requiring inference strategies that consider regional differences in laminar architecture. Limitations of laminar inference for multiple simultaneous sources Having established that a single source can be accurately localized in depth given sufficient SNR and co-registration accuracy, we next examined whether two simultaneous sources of equal strength at different layers could be differentiated. To test this, we simulated two sources at each of the 100 cortical locations: one fixed in the middle layer and the other varying from the most superficial to the deepest layer surface (SNR = −20 dB, co-registration error = 0 mm). Across all simulations, rather than producing two distinct peaks in relative model evidence, a single peak consistently emerged at an intermediate depth between the two sources, both in mean free energy ( Figure 5a ) and in the peak free energy model probability density ( Figure 5b ). This suggests that the presence of the middle-layer source systematically biased inference toward the middle-layer model compared to single-source simulations ( Figure 5e,f ). However, further analysis of the free energy distribution across models ( Figure 5c ) and across simulations ( Figure 5d ) suggests that distinct contributions from both sources may still be present. Notably, model evidence for the middle-layer surface was higher relative to single-source simulations ( Figure 5g,h ), indicating that while the current approach does not fully resolve two simultaneous sources, further refinement of the method may improve discriminability. Download figure Open in new tab Figure 5. Two-source simulations at −20 dB SNR highlight the difficulty of separating simultaneous sources at different depths. a,b) When one source was fixed at the middle layer and a second source was varied from superficial to deep, the mean free-energy (ΔF) matrices and peak-probability plots showed a single intermediate peak, indicating bias toward the middle-layer model (S = superficial surface layer, D = deep surface layer). c,d) Despite this bias, plots of ΔF across simulated layers (c) and across forward-model evaluations (d) reveal that deeper and more superficial contributions still influence the model evidence. e-h) In contrast, single-source simulations at the same SNR showed stronger diagonality in ΔF (e,f) and clearer separation in layer-by-layer analyses (g,h), suggesting that additional refinements are needed to fully disentangle two concurrent laminar sources. Cortical column orientation estimation influences laminar inference accuracy We assessed the influence of dipole orientation estimation on laminar source reconstruction by evaluating four methods 34 : surface normals from the original high-resolution mesh (OSN), surface normals from the downsampled mesh (DSN) 38 , cortical patch statistics (CPS) 39 , and link vectors (LV) computed between corresponding vertices on the pial and white matter surfaces 33 . For each method except link vectors, we included variants with orientations fixed across layers (“fixed”) and orientations independently computed per layer (link vectors have the same orientation at corresponding vertices across surfaces by definition). This resulted in seven variations in dipole orientation definition. Figure 6a illustrates these methods applied across cortical depth at representative vertices. Download figure Open in new tab Figure 6. Impact of dipole orientation estimation on laminar source reconstruction accuracy. a) Example dipole orientation vectors at a subset of cortical locations across layers for four methods: link vectors (LV), original surface normals (OSN), cortical patch statistics (CPS), and downsampled surface normals (DSN). b) Relationship between orientation mismatch (angular difference in degrees) and laminar reconstruction error, over all simulation-reconstruction method pairs and cortical depths. c-e) Angular deviation (median across vertices) between orientation estimation methods at three cortical depths: superficial ( c , surface layer 1), middle ( d , surface layer 6), and deep ( e , surface layer 11). f-h) Corresponding reconstruction error (in layers) for each simulation-reconstruction method pair at the same depths, averaged across 100 cortical locations. Errors are lowest along the diagonal, where simulation and reconstruction methods match, and largest when methods differ substantially (LV = link vectors; OSN = original surface normals; CPS = cortical patch statistics; DSN = downsampled surface normals. The suffix “f” indicates dipole orientations were fixed across layers). To quantify differences between orientation vectors, we computed the angular deviation between each pair of methods at each cortical depth. These deviations were highest between LV and surface-normal-based approaches ( Figure 6c-e ). The fixed variants of the surface-normal-based approaches used the surface normals computed from the pial surface, and therefore the angular deviation between fixed and non-fixed variants of the same family increased with cortical depth ( Figure 6c-e ). We then tested how these orientation discrepancies affected reconstruction performance by simulating laminar sources using each method and reconstructing them with each of the others. Reconstruction error, defined as the deviation between the inferred and ground truth layer, was lowest when the same method was used for both simulation and reconstruction ( Figure 6f-h ), consistent with prior findings 34 . Errors were highest when simulation and reconstruction orientations were mismatched, particularly across different estimation families (e.g., LV vs DSN). Across cortical depths, we observed a strong positive and nonlinear relationship between angular discrepancies in dipole orientation and errors in laminar source reconstruction ( Figure 6b ), indicating that even modest orientation misalignment can degrade depth precision. Importantly, we found no consistent advantage for fixed-orientation strategies over those allowing depth-varying orientations, suggesting that anatomical consistency across layers is less critical than aligning dipole orientations with the true cortical column geometry. These results underscore the importance of accurately estimating dipole orientation when performing laminar inference in real data, and caution against using anatomically implausible orientation models. Anatomical features determine the accuracy of laminar inference The previous analyses demonstrate that laminar inference accuracy critically depends on signal quality, co-registration precision, and accurate estimates of cortical column orientation. However, even with optimal signal-to-noise ratio (≥ −20 dB) and zero co-registration error, reconstruction accuracy is not uniform across the cortex 26 . This variability suggests that local anatomical properties may impose inherent constraints on laminar reconstruction. To evaluate how anatomical variability influences laminar source reconstruction, we simulated, one at a time, sources at every vertex of a multilayer cortical mesh (29,130 locations × 11 layers) under these idealized conditions. Layer depth was inferred using model evidence across all cortical layers. Reconstruction error, defined as the absolute difference between the simulated and inferred depth, was zero for 24,699 out of 29,130 vertices (84.79%), indicating that laminar inference is feasible across most of the cortex under these conditions. However, reconstruction accuracy varied regionally, with errors concentrated in anatomically constrained areas ( Figure 7a ). Download figure Open in new tab Figure 7. Anatomical factors influencing laminar reconstruction error. a-b) Inflated-surface maps of laminar reconstruction error (a) and lead-field RMS difference (deep vs. superficial; b) c-h ) Inflated-surface maps and joint density and marginal histograms for each anatomical metric versus lead-field RMS difference: and dipole orientation (c,d), scalp-to-cortex distance (e,f), and cortical thickness (g-h). Joint density plots are color-coded by laminar reconstruction error (mm). Lower lead-field differences, tangential orientations, greater distance to the scalp, and thinner cortex are associated with increased laminar-reconstruction error. To understand the anatomical basis for this regional variation, we examined four candidate factors: lead field separation between deep and superficial layers ( Figure 7b ), dipole orientation relative to the scalp ( Figure 7c ), scalp-to-cortex distance ( Figure 7e ), and cortical thickness ( Figure 7g ). Although none of these features perfectly separated vertices with zero error from those with nonzero error, all were significantly associated with reconstruction error. Specifically, vertices with low laminar discriminability exhibited smaller lead field differences between depths, more tangential orientations, greater distances from sensors, and thinner cortices. Joint density plots and marginal histograms revealed that these relationships were nonlinear ( Figure 7d,f,h ). Among all features, the dissimilarity between deep and superficial lead fields emerged as the strongest predictor of error (optimal threshold = 0.09 a.u.; AUC = 0.756), suggesting that the ability to distinguish layer-specific contributions at the sensor level places a fundamental constraint on laminar inference. Cortical thickness (optimal threshold = 2.91 mm; AUC = 0.749) and scalp distance (optimal threshold = 24.43 mm; AUC = 0.689) also contributed, consistent with prior work on MEG sensitivity profiles, whereas dipole orientation had a weaker but still detectable effect (optimal threshold = 0.71 a.u.; AUC = 0.595). These findings suggest that even under favorable conditions, anatomical constraints shape the spatial distribution of inference precision. To determine whether combinations of anatomical features better predict laminar inference accuracy than single features, we performed a multivariate logistic regression analysis. The logistic regression model significantly discriminated between error-free and erroneous reconstructions (mean cross-validated AUC = 0.810 ± 0.013), with the largest positive coefficient assigned to lead field difference ( β = 10.64), followed by orientation relative to the scalp ( β = 1.01) and cortical thickness ( β = 0.76); scalp-to-cortex distance had minimal influence ( β = −0.02). Adding pairwise interaction terms provided only a negligible improvement (mean cross-validated AUC = 0.812 ± 0.015). A random forest classifier corroborated these findings, highlighting the lead field difference as the strongest predictor (importance = 32%), followed by thickness (27%), distance (22%), and orientation (19%), though it did not enhance prediction accuracy beyond logistic regression (mean cross-validated AUC = 0.804 ± 0.013). Together, these results indicate that anatomical constraints on laminar inference accuracy primarily arise from additive contributions of distinct features rather than from complex interactions. Visual ERFs show expected pattern of feedforward laminar activation Having established through simulations the conditions under which laminar inference is reliable, we next applied multilayer hpMEG source reconstruction to empirical visual event-related fields (ERFs), using these insights to guide analysis choices. Participants were cued to attend to one visual hemifield, and for each trial we analyzed the VEF separately in the hemisphere contralateral versus ipsilateral to the attended stimulus location 50 . For each condition, the peak of the early visually evoked field (VEF) component was localized to primary visual cortex (V1), and sliding-window laminar inference was performed at a single anatomically and functionally optimal vertex, selected by jointly maximizing signal strength and anatomical suitability (Methods). Vertex selection was used to maximize laminar discriminability rather than to identify biologically privileged cortical locations. The resulting laminar dynamics closely matched established feedforward models of early visual processing. At the peak of the VEF (80-100 ms; M = 90.6 ms), model evidence was dominated by superficial laminae (I-III) together with lamina IV ( Figure 8 ; ipsilateral condition in Figure S6 ), consistent with thalamocortical input to lamina IV and subsequent activation of supragranular pyramidal populations 51 , 52 . This superficial dominance was followed by a later increase in infragranular activity (laminae V-VI) between approximately 110 and 150 ms, possibly reflecting feedback or modulatory processes targeting deep laminae 53 , 54 . Download figure Open in new tab Figure S5. Effect of patch size mismatch on laminar source reconstruction accuracy across SNR levels. a,b) Laminar source reconstruction error (in layers and mm) falls below chance at higher SNR only when the patch size of the source reconstruction model matches that of the simulated data (5mm in these simulations). c,d) Over-estimation of patch size biases laminar inference toward superficial layers, whereas under-estimation slightly biases toward deep layers. Shaded bands represent 95% confidence intervals, and asterisks mark SNR levels at which observed error or bias differed significantly from chance (p < 0.05). Download figure Open in new tab Figure S6. Feedforward laminar activation in ipsilateral primary visual cortex during visuospatial attention. Left: Locations of the selected VEF source vertices for each participant, projected onto an inflated fsaverage cortical surface; each colored sphere denotes one subject. Right: Sliding-window laminar inference applied to visually evoked fields (VEFs) localized to primary visual cortex (V1) during a visuospatial attention task, analyzing responses in the hemisphere ipsilateral to the attended stimulus location. Colored traces show lamina-specific relative model evidence (ΔF, relative to the worst model at each time point) averaged across participants; the dark trace shows the corresponding middle-layer source time series. Horizontal color bars indicate, at each time point, the laminar family with the highest protected exceedance probability (PEP), grouping layers into supragranular (laminae I-III), granular (lamina IV), and infragranular (laminae V-VI) families. The ipsilateral VEF exhibits the same canonical feedforward laminar sequence observed contralaterally, with early granular dominance consistent with thalamic input, followed by supragranular dominance at the VEF peak and subsequent infragranular dominance. Download figure Open in new tab Figure 8. Feedforward laminar activation in primary visual cortex during visuospatial attention. Left: Locations of the selected VEF source vertices for each participant, projected onto an inflated fsaverage cortical surface; each colored sphere denotes one subject. Right: Sliding-window laminar inference applied to visually evoked fields (VEFs) localized to primary visual cortex (V1) during a visuospatial attention task, analyzing responses in the hemisphere contralateral to the attended stimulus location (see Figure S6 for the ipsilateral condition). Colored traces show lamina-specific relative model evidence (ΔF, relative to the worst model at each time point) averaged across participants; the dark trace shows the corresponding middle-layer source time series. Horizontal color bars indicate, at each time point, the laminar family with the highest protected exceedance probability (PEP), grouping layers into supragranular (laminae I-III), granular (lamina IV), and infragranular (laminae V-VI) families. The VEF exhibits a canonical feedforward laminar sequence, with early dominance of the granular family consistent with thalamic input, followed by supragranular dominance at the VEF peak and subsequent infragranular dominance, reflecting later deep-laminae engagement. To assess robustness at the group level, we performed Bayesian model selection 47 using family-wise comparisons of model evidence, grouping laminae into supragranular (I-III), granular (IV), and infragranular (V-VI) families. Protected exceedance probabilities (PEPs) were computed at each time point to quantify the likelihood that one laminar family dominated over the others while accounting for chance-level equivalence 48 . This analysis confirmed a canonical feedforward sequence: early dominance of the granular family (around 60 ms), consistent with initial thalamic input to V1; a subsequent supragranular dominance during the VEF peak (79-104 ms); and a later infragranular dominance (146-155 ms), corresponding to deep-laminae engagement following feedforward processing 6 , 54 – 56 . An equivalent pattern was observed in the ipsilateral hemisphere ( Figure S6 ). These results demonstrate that multilayer hpMEG recovers the expected laminar sequence of early visual processing in human V1, validating the empirical applicability of the multilayer framework and its ability to resolve feedforward laminar dynamics non-invasively. Motor ERFs reveal characteristic laminar activation sequences We next applied multilayer hpMEG source reconstruction to motor event-related fields (ERFs) using two independent datasets: one involving visually cued button presses and another involving self-paced button presses. In each dataset, we localized the motor field (MF) component to the hand area of primary motor cortex (M1) and the motor-evoked field I (MEFI) component to the corresponding region of primary somatosensory cortex (S1), and performed sliding-window laminar inference at each site. Qualitative inspection revealed highly similar laminar dynamics across the two datasets for both components ( Figure S7 ), motivating their combination to enable group-level Bayesian model selection, which was not feasible within either dataset alone due to limited sample size. Download figure Open in new tab Figure S7. Reproducible laminar activation patterns in motor and somatosensory cortices across tasks. a-d) Sliding-window laminar inference applied separately to visually cued ( a,b ) and self-paced ( c,d ) button presses. Motor field (MF) activity in primary motor cortex (M1) is shown in (a) and (c), and motor-evoked field I (MEFI) activity in primary somatosensory cortex (S1) in (b) and (d). Left panels: Locations of the selected source vertices for each participant, projected onto an inflated fsaverage cortical surface; each colored sphere denotes one subject. Right panels: Lamina-specific relative model evidence (ΔF) and the corresponding middle-layer source time series. Across both tasks and both cortical regions, qualitatively similar laminar dynamics are observed, with dominant infragranular activity surrounding movement onset in M1 followed by delayed supragranular activity, and predominantly infra- and then supra-granular responses in S1 at movement onset with additional superficial transients at later latencies. Using the combined dataset, we performed family-wise Bayesian model selection for the MF in M1 and the MEFI in S1 ( Figure 9 ). In M1, laminar inference revealed a reproducible sequence dominated by infragranular activity (laminae V-VI) beginning approximately 150 ms before movement onset, consistent with the contribution of deep-laminae pyramidal neurons projecting to the corticospinal tract to the generation of the MF component 57 – 59 . This deep-laminae dominance recurred in two additional transients, coinciding with the button press and shortly thereafter. A later emergence of supragranular activity (laminae I-III), approximately 300-350 ms after the button press, followed the deep-laminae transients. This temporal pattern suggests successive phases of motor preparation, corticospinal output, and post-movement reactivation, with the delayed superficial response potentially reflecting feedback or efference-copy-related processing within local cortical circuits. Download figure Open in new tab Figure 9. Laminar activation sequences of motor and somatosensory ERFs during movement. a) Motor field (MF) activity localized to the hand area of primary motor cortex (M1). b) Motor-evoked field I (MEFI) activity localized to primary somatosensory cortex (S1). Left panels: Locations of the selected source vertices for each participant, projected onto an inflated fsaverage cortical surface; each colored sphere denotes one subject. Right panels: Sliding-window laminar inference applied to the combined dataset of visually cued and self-paced button presses. Colored traces show lamina-specific relative model evidence (ΔF, relative to the worst model at each time point) averaged across participants; dark traces show the corresponding middle-layer source time series. Horizontal color bars indicate, at each time point, the laminar family with the highest protected exceedance probability (PEP), grouping layers into supragranular (laminae I-III), granular (lamina IV), and infragranular (laminae V-VI) families. MF responses in M1 are characterized by dominant infragranular activity preceding, during, and shortly after movement, followed by delayed supragranular activation, whereas MEFI responses in S1 show predominant infra- and then supra-granular activity at movement onset and subsequent deep and superficial transients. Laminar inference of the MEFI in S1 showed a distinct profile, characterized by predominant infa- and then supra-granular activity at movement onset and additional deep and superficial transients at later latencies. This pattern is consistent with prior work identifying the MEFI as a short-latency (∼40 ms) response to proprioceptive and cutaneous reafference in the postcentral gyrus 57 , 60 . The close correspondence between laminar dynamics observed in the two original datasets ( Figure S7 ), together with their convergence in the combined analysis, demonstrates that multilayer hpMEG robustly reproduces established sensorimotor activation motifs and extends them into the laminar domain, revealing consistent deep-superficial activation sequences in human motor circuits. Discussion Building upon previous models of laminar inference in high-precision MEG (hpMEG) that compared deep versus superficial sources, we introduce a multilayer source reconstruction framework that enables depth-resolved inference across the full cortical column. Simulations establish feasibility limits for laminar inference, anatomical constraints guide vertex selection, and Bayesian model comparison provides statistical inference over cortical depth. Our simulations reveal that this precision requires a sufficiently high signal-to-noise ratio (approximately −35 dB for deep-versus-superficial distinctions and −20 dB for full laminar inference) and minimal co-registration error (< 2mm), conditions that can be met in practice through the use of individualized head-casts and optimized recording protocols 27 , 28 . Although coarse deep-versus-superficial distinctions remain feasible at lower SNRs with two-layer source models, true six-laminae resolution requires fine-grained depth coverage to account for variations in laminar thickness across the cortex. Crucially, we confirm that correctly estimating cortical column orientation is essential for accurate laminar inference, as errors in dipole orientation estimation systematically degrade localization accuracy. Corroborating previously established results 26 , we show that anatomical features of the cortex such as thickness, dipole orientation, scalp-to-cortex distance, and lead field separation across depths nonlinearly shape the spatial distribution of inference accuracy, with lead field dissimilarity emerging as the strongest predictor. Beyond simulations, we demonstrate the empirical viability of multilayer hpMEG by recovering canonical laminar activation sequences in human visual and sensorimotor cortex. Notably, visual and motor ERFs exhibited complementary laminar sequences, with early granular dominance in sensory cortex and early infragranular dominance in motor cortex, reflecting their distinct circuit architectures and afferent/efferent nature, respectively. In primary visual cortex, visually evoked fields exhibited a feedforward laminar progression, with early granular dominance consistent with thalamocortical input, followed by supragranular activity at the response peak and subsequent infragranular engagement, with the same feedforward laminar sequence observed in both contralateral and ipsilateral hemispheres. In motor circuits, combining two independent datasets revealed highly reproducible laminar dynamics: motor fields in primary motor cortex were dominated by infragranular activity beginning before movement onset and recurring around execution, consistent with corticospinal output, followed by delayed supragranular responses, while motor-evoked responses in primary somatosensory cortex showed predominant superficial activity at movement onset, consistent with somatosensory reafference. The convergence of these laminar patterns across cortical systems, task variants, and acquisition sites demonstrates that multilayer hpMEG not only meets theoretical requirements for laminar precision but also recovers well-established physiological motifs in empirical human data, extending non-invasive laminar inference beyond deep-versus-superficial distinctions. Despite the significant improvements offered by our multilayer approach, a key limitation is the inability of current model-based inference methods to resolve two simultaneous sources of equal magnitude at different depths. While exact temporal alignment and perfectly matched dipole moments are unlikely in biological systems, this limitation highlights the constraints of relying solely on Bayesian model comparison. Specifically, two closely spaced dipoles at different depths may generate a sensor-level field pattern indistinguishable from a single intermediate dipole, limiting the spatial discrimination achievable by current MEG systems. One potential solution would be to use sensors capable of measuring higher spatial-frequency components of the magnetic fields, such as optically pumped magnetometers (OPMs), which could increase sensitivity to subtle spatial differences between closely spaced sources 61 . Future refinements should explore both sensor technology advances and alternative inference strategies, such as Bayesian model averaging (BMA), which assigns probabilistic weights to multiple competing models rather than selecting a single best-fit model 62 . These combined approaches may enhance the simultaneous resolution of multiple sources, overcoming a fundamental limitation of the current winner-takes-all approach. A related consideration is the spatial extent of the generative source. In additional simulations, we simulated 5 mm patches of cortical activity and assessed reconstruction performance using source models with varying patch sizes (2.5 mm, 5 mm, and 10 mm FWHM). We found that mismatches between the true patch size and the extent assumed in the generative model systematically biased reconstruction: underestimation led to deeper solutions, whereas overestimation biased inference toward superficial layers ( Figure S5 ). These results align with prior findings 25 , 26 , and reinforce the importance of accurate assumptions about the spatial spread of the source, which could vary between laminae and depends on factors such as lateral cortical connectivity. In supplementary simulations using different sliding window sizes (15, 25, 50, and 100 ms; Figure S1 ), we assessed the temporal sensitivity of laminar model evidence. Across all window lengths, model-based inference consistently recovered the true laminar depth, but the temporal profile of model evidence differed. Narrow windows (e.g., 15 ms) yielded sharper but noisier evidence profiles that collapsed earlier, while wider windows (e.g., 100 ms) integrated more information, resulting in broader and more stable evidence plateaus over time. This trade-off reflects the balance between temporal resolution and sensitivity: short windows enable finer tracking of rapid changes but require high SNR, whereas longer windows are more robust to noise but may obscure transient laminar dynamics. As demonstrated above, model-based laminar inference requires careful selection of parameters and prior models, and in practice necessitates restricting analysis to a small number of vertices that are both anatomically favorable and functionally informative. A more direct, though equally technically challenging, approach would be to extract source signals directly from each layer of a multi-layer surface model across the cortex, rather than relying on model-based inference at specific vertices. Such an approach could be computationally tractable at the whole-cortex level and might resemble laminar local field potential (LFP) recordings, potentially enabling direct current source density (CSD) estimation and facilitating direct comparisons with laminar LFP data and laminar biophysical neural modeling 31 , 63 – 65 . However, even laminar LFP recordings are extracellular and thus remain subject to an inverse problem due to aggregation of diverse signal sources 66 that is typically ignored in conventional CSD analyses 14 . Moreover, applying CSD-like methods 67 , 68 to hpMEG presents significant challenges, including sign ambiguity 69 and sensitivity to noise 70 . A practical consideration for real-world applications is the effect of dynamic head motion on co-registration error. In our simulations, co-registration error was modeled as static for computational efficiency, but in real recordings, even small head movements may introduce depth-localization errors. The most effective way to mitigate motion is through individualized head-casts, which substantially reduce head displacement and can realistically achieve < 2 mm co-registration error 27 , 28 . Alternatively, wearable customized arrays of OPMs could achieve similar precision by moving with the head, thus inherently reducing motion-induced errors 71 . While head-casts or wearable OPM arrays do not eliminate movement entirely, residual movement artifacts can be addressed using additional correction strategies, including regression of motion-related signals 72 or temporally extended source space separation (tSSS), which virtually transforms the sensor array to move with the head 73 , 74 , and was applied in the analysis of the button-press datasets presented here. Combining physical stabilization with signal-level correction may be necessary to reach the precision required for robust laminar inference in practice. The ability to resolve sources with laminar precision necessitates a reevaluation of the theoretical limits of MEG resolution. The estimated number of pyramidal cells required to generate a detectable MEG signal ranges from 10,000 to 50,000, with the majority of the signal arising from laminae II/III and V 75 . Neuronal density varies substantially across the cortex, from approximately 20,000 cells/mm 2 in frontal regions to around 45,000 in sensorimotor areas and up to 90,000 in the foveal region of primary visual cortex (V1) 76 . Factors such as cortical thickness, gyrification, and orientation of cortical columns relative to the sensors further constrain spatial resolution and depth precision. However, even with conventional MEG recording (without head-casts), spatial resolution as fine as ∼1 mm 2 in V1 77 , 78 and ∼3.5 mm 2 in the sensorimotor cortex 79 has been demonstrated. Given the enhancements offered by hpMEG, achieving laminar resolution is a realistic goal. If it is possible to resolve a source on the cortical surface with millimeter precision, the same principles should, in theory, allow for distinguishing sources at different cortical depths provided the right conditions are met. Our results demonstrate that constant source orientation across layers may not be appropriate for laminar inference, opening the possibility of using piecewise estimates of cortical column curvature. Cortical columns are not strictly straight 80 , and future models could incorporate curvature by combining depth-varying orientation vectors derived from equivolumetric laminar surfaces 45 , as opposed to the equidistant surfaces used here. Estimating cortical column orientation using diffusion tensor imaging (DTI) is a promising alternative, as anisotropic diffusion within cortical gray matter can reflect the organization of cortical columns 81 . However, DTI-based estimates can be confounded by white matter fibers projecting into the cortex (the gyral bias effect), potentially biasing orientation measures. Advanced diffusion imaging techniques such as high angular resolution diffusion imaging (HARDI) and constrained spherical deconvolution, combined with surface-based analyses that sample diffusion metrics along radial cortical trajectories, can help isolate cortical diffusion signals from adjacent white matter influences 82 , 83 . Integrating these advanced imaging approaches into laminar MEG source models could yield more biologically plausible dipole orientation estimates and potentially improve depth localization, especially in highly folded cortex. The current framework provides a foundation for evaluating such models and for testing how column curvature influences laminar MEG source reconstructions. Given these constraints, it is natural to ask how hpMEG compares to ultra-high field fMRI (7 - 11.7 Tesla) for laminar investigations. While fMRI provides spatially precise laminar profiles, it is fundamentally limited by its reliance on the hemodynamic response, which is subject to vascular drainage effects 84 that can result in superficially biased laminar profiles. Efforts to mitigate this issue, such as deconvolution 85 , vascular correction techniques 86 , and contrast-based experimental designs 87 , have improved laminar fMRI but do not fully address this limitation. Another important fMRI limitation is its low temporal resolution (e.g. ∼12 s per whole-brain acquisition at 0.8 mm resolution 88 ). As a result, many studies restrict laminar fMRI acquisition to a limited part of the brain in order to achieve better temporal resolution. hpMEG, in contrast, provides millisecond-resolution measurements of neuronal population activity with whole-head coverage and is unaffected by vascular artifacts, making it uniquely suited for studying activity with rapid laminar dynamics such as beta bursts 29 , 89 , induced oscillatory activity 30 , and event-related responses. However, because MEG requires anatomical imaging to define laminar boundaries, its accuracy depends on high-quality MRI-derived cortical reconstructions. Although recent attempts have been made to generate individualized laminar maps over the whole cortex using multimodal MRI and machine learning 90 – 94 , accurately achieving this for the whole cortex remains a major technical challenge. Successfully addressing this issue will require integrating advanced multimodal imaging approaches, ultra-high-field MRI, movement correction, and sophisticated computational methods. Current hpMEG analyses must therefore rely on laminar histological atlases based on a single, ex vivo brain such as the BigBrain atlas 44 , 45 . Developing methods for individualized laminar boundary estimation in vivo would represent a major step forward in both hpMEG and laminar fMRI research. Given that lead field separation between superficial and deep cortical layers emerged as the strongest anatomical predictor of laminar inference accuracy, our results highlight the importance of accurate forward modeling. Our analyses, and all previous laminar MEG studies, relied on a simple single-shell forward model 37 , which may underestimate differences in lead fields across cortical depths. Future work could improve laminar inference precision by employing more sophisticated forward models, such as boundary element models (BEM) 95 or finite element models (FEM) 96 , which explicitly incorporate tissue-specific conductivity profiles. Recent methodological advancements improve FEM by integrating DTI-based conductivity tensors 97 , offering another compelling route to combine MEG with DTI data to refine laminar specificity in source localization. Emerging OPM technologies hold promise for further enhancing laminar inference by offering improved signal quality, allowing flexible sensor placement, and reducing operational costs 61 . While still in its early stages, simulation studies suggest that OPM-based MEG, if configured correctly, can achieve laminar resolution comparable to classic, SQUID-based MEG, at least for bilaminar models 98 . Future work should assess whether these benefits extend to true laminar inference. Both our results and those of Helbling 98 suggest that incorporating individual anatomical features such as cortical thickness and curvature could optimize OPM sensor placement to capture higher spatial-frequency components of the magnetic fields or achieve higher effective SNRs, potentially improving laminar precision within targeted cortical regions. Simulations could help refine these configurations, eliminating the need for time-consuming trial-and-error adjustments. Such targeted precision may be particularly useful for pre-surgical mapping of eloquent cortical regions or for resolving lamina-specific computations in particular regions of interest. In summary, we establish hpMEG as a viable tool for non-invasive laminar inference, defining the conditions under which true depth-resolved localization is possible. By systematically exploring the trade-offs between signal quality, anatomical constraints, and inference accuracy, we delineate best-case and marginal conditions for MEG-based laminar inference. While current implementations require restrictive priors and high SNR, future developments, including Bayesian model averaging, individualized laminar and orientation priors, more accurate forward models, and emerging sensor technologies, may further improve hpMEG’s ability to resolve cortical computations at the laminar level. These advances pave the way for a new era of laminar human neuroscience, bridging invasive electrophysiology with whole-brain, millisecond-resolution functional imaging. Funding This work was supported by a Fondation pour l’Audition grant RD-2022-03 to J.J.B. and A.R.D., a European Research Council (ERC) consolidator grant 864550 to J.J.B., a European Research Council (ERC) starting grant 716862 to M.B., and a French National Research Agency (ANR) grant TheAlphaMEG (2025-2030, ANR-20-CE17-0023) to M.B. Author contributions Conceptualization: J.J.B., S.B., G.R.B. Methodology: J.J.B., M.J.S., M.M., D.Sh., M.B., S.D., F.L., D.Sc. Software: J.J.B., M.J.S., D.Sh., I.A., Q.M., S.G. Validation: J.J.B., M.J.S., M.M., D.Sh., I.A. Formal analysis: J.J.B., M.J.S., M.M., D.Sh., I.A. Investigation: J.J.B., M.J.S., C.F.P., Y.Z., F.L., S.D. Visualization: J.J.B., M.J.S., M.M., I.A., D.Sh. Writing—original draft: J.J.B., M.J.S., M.M., D.Sh., I.A., Q.M., S.G. Writing—review & editing: All authors. Supervision: J.J.B., A.R.D., M.B. Project administration: J.J.B. Funding acquisition: J.J.B., A.R.D., M.B. Competing interests The authors declare that they have no competing interests. Data and materials availability All data and code used in this study are available at https://github.com/danclab/multilaminar_sim/ , https://github.com/danclab/laminar_motor_erf/ , https://github.com/danclab/auditory_laminar , and https://github.com/cophyteam/laminar_visual_erf . Additional information on methods and additional figures are provided in the Supplementary Materials. Supplementary Materials Visually evoked response dataset The visually evoked response dataset was collected at the CERMEP imaging platform (Lyon, France). A total of 28 healthy participants were included in the study (11 male, aged 25.7 ± 3.86 years). All participants were right-handed, as assessed by the Edinburgh Handedness Inventory 1 , and reported having either normal or corrected-to-normal vision. For each participant, structural MRI, MEG recordings, and eye-tracking data were acquired. MRI data were acquired on a 3T Siemens Sonata scanner (CERMEP, Lyon, France) using three protocols. The first protocol acquired a T1-weighted MPRAGE sequence (TR = 2100 ms, TE = 3.33 ms, TI = 900 ms, flip angle = 8°, 1mm isotropic resolution, 256 × 256 mm field of view, 192 slices) and these images were used to create participant-specific head-casts to reduce between-session co-registration error and within-session movement during MEG recordings. Head-casts were created using a 3D scalp surface extracted from the T1-weighted volume using FreeSurfer (v6.0.0) 2 , 3D-printed, placed in a dewar replica, and filled with polyurethane foam. Recesses for nasion and preauricular (LPA/RPA) fiducial coils were incorporated to enable coil locations to be recovered in the participant’s native MRI space from the head-cast CAD model 3 . The second protocol acquired a sagittal 3D T2-weighted SPACE sequence (TR = 3200 ms, TE = 299 ms, flip angle = 120°, 0.8 mm isotropic resolution, 280 × 320 mm field of view). The third protocol used a quantitative multi-parameter mapping (MPM) protocol optimized for cortical surface reconstruction 4 , comprising multi-echo 3D FLASH acquisitions with proton density (PD) weighting (flip angle = 6°, TE = 2.3-18.4 ms in eight echoes), T1 weighting (flip angle = 21°, TE = 2.3-18.4 ms in eight echoes), and magnetization transfer (MT) weighting (flip angle = 6°, TE = 2.3-13.8 ms in six echoes), each acquired with TR = 25 ms at 0.8 mm isotropic resolution (280 × 320 mm field of view) using parallel imaging (GRAPPA, factor 2). Calibration data included B1 mapping and gradient-echo field mapping to enable correction of transmit-field inhomogeneity and geometric distortions. Quantitative maps (PD, R1, MT, and R2*) were computed using the hMRI toolbox 5 . These maps, as well as the T2-weighted images, were aligned to the T1 used to create the head-cast, and a custom FreeSurfer-based (v6.0.0) 2 reconstruction procedure 6 was used to generate pial and white/grey matter boundary surfaces from the MPM-derived PD and T1 contrasts together with the T2-weighted volume. MEG data was collected using a CTF Omega 275 channel system in a magnetically shielded room at a sampling rate of 1200 Hz. Eye tracking data was recorded using an EyeLink 1000+ at a sampling rate of 1000 Hz. Participants performed a selective visual orientation-discrimination task while in the supine position. On each trial, participants were simultaneously presented with two vertical grating stimuli (diameter 3.5°, 12 cycles across) located laterally at an eccentricity of 3.7° of visual angle from the screen center. The task was to determine the orientation of the target grating relative to an implicit 45° tilted reference diagonal (not displayed on screen), while ignoring the distractor on the other side of the screen. Stimuli were projected on a screen positioned at approximately 50 cm from the subject. Participants responded using a two button button-box using either the index or middle finger of their right hand. The experiment was structured into four blocks of 300 trials each, with the side of the target (left/right) balanced across blocks. A central cue (a leftward or rightward pointing arrow) was presented at the beginning of each block, indicating the side of the target stimulus for that block. MEG data was preprocessed using MNE Python (v.1.9) 7 . Raw data was band-pass filtered between 0.1 and 80 Hz, line noise was removed by a zero-phase overlap-add FIR filter at 50 Hz, and the data were then downsampled to 600 Hz and epoched around stimulus onset (−0.5 s to 0.35 s). To ensure high data quality, epochs were excluded if the participant’s gaze deviated by more than 1.5° of visual angle from the fixation cross prior to stimulus onset, if they contained blinks or other pronounced ocular artifacts, or if they were associated with incorrect behavioral responses or premature responses (reaction time < 250 ms). Following these automated procedures, all remaining epochs were visually inspected, and additional trials contaminated by noise or movement-related artifacts were manually removed. Sensor-level SNR was computed within the analysis window as the ratio in dB of the RMS of the trial-averaged signal across channels to the RMS of trial-by-trial deviations from that mean, yielding a single SNR time series per participant. Participants whose peak SNR within this window was < −15 dB were excluded (n = 8). Trials were analyzed separately for the left- and right-attention conditions. After preprocessing, this left an average of 187.1 trials per subject (SD = 52.5), including 89.7 trials (SD = 27.7) for the left-attention condition and 97.4 trials (SD = 28.5) for the right-attention condition. Localization was performed in each hemisphere within primary visual cortex (Brodmann area 17), defined using each participant’s FreeSurfer Brodmann Areas ex vivo atlas (BA_exvivo) 8 , over 50 to 150 ms relative to the onset of the visual grating, in line with the expected visual response timing 9 . Sliding time window model comparison was performed with overlapping windows of 10 ms. All preprocessing and laminar analysis code for this dataset can be found at https://github.com/cophyteam/laminar_visual_erf . Visually cued button-press dataset The visually cued button-press dataset was collected at the University College London Wellcome Centre for Human Neuroimaging. A total of eight healthy participants were included in the study (six male, aged 28.5 ± 8.52). For each participant, structural MRI, MEG recordings, and eye-tracking data were acquired. MRI data were acquired on a 3T Siemens Magnetom TIM Trio system using the body coil for RF transmission and a 32-channel head coil for reception in two separate protocols. The first protocol generated a precise scalp image for head-cast construction 3 using a T1-weighted 3D FLASH sequence (1 mm isotropic resolution, field of view = 256 × 256 × 192 mm [A-P, H-F, R-L], TR = 7.96 ms, flip angle = 12°, single echo, bandwidth = 425 Hz/pixel, acquisition time = 3 min 42 s). A spacious 12-channel head coil was used without padding or headphones to avoid skin displacement. Participant-specific head-casts were created from a 3D scalp model extracted from the T1-weighted images using SPM 3 , 10 ; the model was 3D-printed, placed in a dewar replica, and filled with polyurethane foam. Recesses for nasion and preauricular (LPA/RPA) fiducial coils were incorporated so that coil locations could be determined directly in native MRI space from the head-cast CAD model 3 . The second protocol used a quantitative multi-parameter mapping (MPM) protocol optimized for cortical surface reconstruction 4 . Three RF- and gradient-spoiled multi-echo 3D FLASH acquisitions with predominantly PD, T1, and MT weighting were acquired at 800 μm isotropic resolution (field of view = 256 mm H-F, 224 mm A-P, 179 mm R-L). Flip angles were 6° (PD), 21° (T1), and 6° (MT). Gradient echoes were acquired with alternating readout polarity at eight equidistant echo times ranging from 2.34 to 18.44 ms (ΔTE = 2.30 ms; bandwidth = 488 Hz/pixel); six echoes were acquired for the MT-weighted volume to maintain TR = 25 ms across contrasts. MT weighting was achieved using a 4 ms Gaussian RF pulse applied 2 kHz off-resonance with a nominal flip angle of 220° prior to excitation. Parallel imaging (GRAPPA, acceleration factor 2 in both A-P and R-L directions, 40 integrated reference lines) accelerated acquisition. Additional calibration sequences mapped transmit field inhomogeneities and corrected geometric distortions due to B0 inhomogeneity. Total acquisition time was <30 min 11 , 12 . Quantitative maps (PD, R1, MT, and R2*) were computed according to established procedures 4 . These maps were aligned to the T1 used to create the head-cast, and a custom FreeSurfer-based (v6.0.1) 2 reconstruction procedure was used to generate pial and white/grey matter boundary surfaces from the PD and T1 volumes of the MPMs 6 . MEG was recorded at 1200 Hz on a 275-channel Canadian Thin Films (CTF) MEG system with super-conducting quantum interference device (SQUID)-based axial gradiometers (VSM MedTech, Vancouver, Canada) in a magnetically shielded room. Participants performed a visually cued action-selection task (random-dot kinematogram, 2 s; 500 ms delay; left/right arrow instruction; speeded left/right button press) while seated upright and wearing a participant-specific head-cast to stabilize head position 3 . Subjects completed 1-6 sessions on different days; each session comprised three 15-min blocks of 180 trials for 540-2700 trials per participant (M = 1822.5). Two sessions each were excluded from three participants due to technical issues during data acquisition, and a total of six sessions across three participants were excluded from the motor field (MF) component analysis due to having components that localized outside of the primary motor cortex. Visual stimuli were rear-projected at ∼50 cm, and responses were made on a button box with the right index (left button) or middle (right button) finger. All acquisition and task materials have been described previously 11 . Raw MEG was movement-compensated with Maxwell filtering (Signal Space Separation, SSS) 13 , 14 using continuously tracked head-position indicator (cHPI) coils. For each run, cHPI coil locations were extracted and transformed from CTF device coordinates to MNE’s device frame; time-resolved head position was then estimated, and we generated quality-control plots of coil trajectories, inter-coil distances, and head-position traces/fields. The device→head transform from the first run of each session was saved and used as the common destination so that subsequent runs were realigned to the same head coordinate frame (three blocks referenced to the first). SSS was applied in head coordinates with temporal extension (st_duration = 10 s) and a fixed spherical origin at [0, 0, 0.04] m (head frame). This procedure yields a single movement-corrected raw file per block that is co-registered across blocks within a session. Data were then downsampled to 600 Hz and line noise was removed using an iterative ZapLine procedure applied until residual 50 Hz power was eliminated 15 . Independent component analysis (ICA; Infomax, 25 components) was fit on a copy of the data bandpass-filtered between 1 and 60 Hz. Components reflecting ocular, cardiac, or muscle artifacts were manually identified and the corresponding unmixing matrix was then applied to the original, unfiltered data to remove them. The cleaned data were then bandpass-filtered between 0.1 and 30 Hz for event-related field analysis and epoched around button press events (−2 s to 2 s). Trials with missing responses were excluded. Epochs from each block were then concatenated within session, and the concatenated epochs were averaged to obtain motor ERFs for laminar inference (yielding one session-level average per subject). For the motor field (MF) component of the motor ERF, localization was performed in the left precentral gyrus, defined using each participant’s FreeSurfer Desikan-Killiany cortical parcellation, over −75 to −25 ms relative to the button press. For the motor-evoked field I (MEFI) component, localization was performed within the left postcentral gyrus, defined analogously, using a post-movement window of 15 to 55 ms. All preprocessing and laminar analysis code for this dataset can be found at https://github.com/danclab/laminar_motor_erf . Self-paced button-press dataset The self-paced button-press dataset was collected at the CERMEP imaging platform (Lyon, France). Data were collected from five neurologically healthy adults (one male; 32.4 ± 5.9 y). For each participant, structural MRI and MEG data were acquired. MRI data were acquired using either a 3T Siemens Sonata system (CERMEP, Lyon, France) or a 3T Siemens Magnetom Vida-XT-64 system (Numaris/X VA20A-04T4; University of Miami, Coral Gables, FL, USA). The structural protocol comprised: (i) a sagittal 3D T1-weighted MPRAGE (TR = 2530 ms, TI = 1100 ms, flip angle = 9°, 0.8 mm isotropic, GRAPPA acceleration factor 2, TE = 3.41 ms in Lyon and 3.42 ms in Miami); (ii) a sagittal 3D T2-weighted SPACE sequence (TR = 3200 ms, TE = 409 ms, 0.8 mm isotropic, variable flip angle mode “T2 var”, acceleration implemented using compressed sensing [factor 3.0] in Miami and GRAPPA/CAIPIRINHA in Lyon to achieve comparable scan duration); and (iii) two sagittal multi-echo 3D FLASH sequences (TR = 20 ms, eight echoes, TE = 1.87-14.59 ms) acquired at flip angles of 5° and 30° with 1 mm isotropic resolution (GRAPPA factor 2). The FLASH images were processed in FreeSurfer to synthesize a 5° FLASH volume, and the resulting high-SNR synthetic image was used to derive the scalp surface for alignment of 3D head scans. FreeSurfer (v6.0.1) 2 was used to reconstruct the cortical pial and white matter boundary surfaces from the T1- and T2-weighted volumes. MEG was recorded at 1200 Hz on a 275-channel Canadian Thin Films (CTF) MEG system with super-conducting quantum interference device (SQUID)-based axial gradiometers (CTF MEG Neuro Innovations, Inc. Coquitlam, Canada) in a magnetically shielded room. In a single session, participants completed two blocks of an auditory oddball task followed by two blocks of a multi-tone masker task, described in detail elsewhere 16 . Between trials, participants initiated the next trial by pressing a button with their right index finger on a response box. Analyses in the present study focused exclusively on these self-paced button-press events, which occurred independently of auditory stimulation. Participants wore individualized head-casts to stabilize head position during MEG acquisition 3 while in the supine position. Unlike the visually cued dataset, the fiducial coils (nasion, LPA, RPA) were taped directly to the participant’s skin to track head motion relative to the cast. To recover the coil positions in native MRI space, each participant’s head and fiducials were digitized using an Artec Eva handheld 3D scanner. The resulting 3D mesh was cropped to the facial region and co-registered to the MRI-derived scalp surface (from the synthetic FLASH image) using a combination of fiducial-based alignment and iterative closest point (ICP) registration implemented in MNE-Python. The final fiducial coordinates were used for MEG-MRI co-registration in subsequent laminar source reconstruction analyses. Each session comprised 316-332 button presses per participant (M = 328.8). Raw MEG was movement-compensated with Maxwell filtering (temporal Signal Space Separation, tSSS) 13 , 14 using continuously tracked head-position indicator (cHPI) coils. For each run, cHPI coil locations were extracted and transformed from CTF device coordinates to MNE’s device frame; time-resolved head position was estimated, and quality-control plots of coil trajectories, inter-coil distances, and head-position traces were generated. The device→head transform from the first run of each session was used as the common destination so that subsequent runs were realigned to a consistent head coordinate frame. SSS was applied in head coordinates with temporal extension (st_duration = 10 s) and a fixed spherical origin at [0, 0, 0.04] m (head frame), yielding one movement-corrected raw file per block. Data were then downsampled to 600 Hz, and line noise was removed using an iterative ZapLine procedure applied until residual 50 Hz power was eliminated 15 . Independent component analysis (ICA; Infomax, 25 components) was fit on a copy of the data bandpass-filtered between 1 and 60 Hz. Components reflecting ocular, cardiac, or muscle artifacts were identified manually and the corresponding unmixing matrix was then applied to the original, unfiltered data to remove these components. The cleaned data were subsequently bandpass-filtered between 0.1 and 30 Hz for event-related field analysis, epoched around button-press events (−2 s to 2 s), and trials with missing responses were excluded. Epochs from all blocks within a session were concatenated, and autoreject 17 was then applied to automatically detect and repair residual artifacts using cross-validated consensus thresholds (consensus = 0-1, n_interpolate = [1, 4, 32]). Cleaned epochs were averaged to obtain motor ERFs for laminar inference, yielding one session-level average per participant. For the MF component of the motor ERF, localization was performed in the left precentral gyrus, defined using each participant’s FreeSurfer Desikan-Killiany cortical parcellation, over −75 to −25 ms relative to the button press; for the MEFI component, in the left postcentral gyrus, defined in the same way, with a 15 to 55 ms post-movement window. Two participants were excluded from the MF analysis and one from the MEFI analysis due to having components that localized outside of the hand area of the primary motor or somatosensory cortex, respectively. All preprocessing and laminar analysis code for this dataset can be found at https://github.com/danclab/auditory_laminar . Acknowledgements We gratefully acknowledge support from the CNRS/IN2P3 Computing Center (Lyon - France) for providing computing and data-processing resources needed for this work. Funder Information Declared European Research Council, https://ror.org/0472cxd90 , ERC-CoG 864550 , ERC-StG 716862 Fondation pour l’Audition , RD-2022-03 French National Research Agency (ANR) , ANR-20-CE17-0023 Footnotes Application of method to event-related fields in visual, somatosensory, and motor cortices. References 1. ↵ Smith , C. U. M. A century of cortical architectonics . J. Hist. Neurosci . 1 , 201 – 218 ( 1992 ). OpenUrl CrossRef PubMed 2. ↵ Hirsch , J. A. & Martinez , L. M. Laminar processing in the visual cortical column . Curr. Opin. Neurobiol . 16 , 377 – 384 ( 2006 ). OpenUrl CrossRef PubMed Web of Science 3. ↵ Harris , K. D. & Shepherd , G. M. G. The neocortical circuit: themes and variations . Nat. Neurosci . 18 , 170 – 181 ( 2015 ). OpenUrl CrossRef PubMed 4. ↵ Bastos , A. M. , Loonis , R. , Kornblith , S. , Lundqvist , M. & Miller , E. K. Laminar recordings in frontal cortex suggest distinct layers for maintenance and control of working memory . Proc. Natl. Acad. Sci . 115 , 1117 – 1122 ( 2018 ). OpenUrl Abstract / FREE Full Text 5. Huang , S. , Wu , S. J. , Sansone , G. , Ibrahim , L. A. & Fishell , G. Layer 1 neocortex: gating and integrating multidimensional signals . Neuron 112 , 184 – 200 ( 2024 ). OpenUrl CrossRef PubMed 6. ↵ Bijanzadeh , M. , Nurminen , L. , Merlin , S. , Clark , A. M. & Angelucci , A. Distinct laminar processing of local and global context in primate primary visual cortex . Neuron 100 , 259 – 274 ( 2018 ). OpenUrl CrossRef PubMed 7. ↵ Dykstra , A. R. , et al. Testing circuit-level theories of consciousness in humans . Trends Cogn. Sci . In press, ( 2025 ). 8. Bastos , A. M. et al. Canonical Microcircuits for Predictive Coding . Neuron 76 , 695 – 711 ( 2012 ). OpenUrl CrossRef PubMed Web of Science 9. ↵ Shipp , S. , Adams , R. A. & Friston , K. J. Reflections on agranular architecture: predictive coding in the motor cortex . Trends Neurosci . 36 , 706 – 16 ( 2013 ). OpenUrl CrossRef PubMed Web of Science 10. ↵ Borbély , S. , Halasy , K. , Somogyvári , Z. , Détári , L. & Világi , I. Laminar analysis of initiation and spread of epileptiform discharges in three in vitro models . Brain Res. Bull . 69 , 161 – 167 ( 2006 ). OpenUrl CrossRef PubMed 11. ↵ Underwood , C. F. & Parr-Brownlie , L. C. Primary motor cortex in Parkinson’s disease: Functional changes and opportunities for neurostimulation . Neurobiol. Dis . 147 , 105159 ( 2021 ). OpenUrl CrossRef PubMed 12. ↵ Pigorini , A. et al. Simultaneous invasive and non-invasive recordings in humans: A novel Rosetta stone for deciphering brain activity . J. Neurosci. Methods 408 , 110160 ( 2024 ). OpenUrl CrossRef PubMed 13. ↵ Chung , J. E. et al. High-density single-unit human cortical recordings using the Neuropixels probe . Neuron 110 , 2409 – 2421 ( 2022 ). OpenUrl CrossRef PubMed 14. ↵ Gratiy , S. L. , Devor , A. , Einevoll , G. T. & Dale , A. M. On the estimation of population-specific synaptic currents from laminar multielectrode recordings . Front. Neuroinformatics 5 , 32 ( 2011 ). OpenUrl PubMed 15. ↵ Caruso , L. et al. In Vivo Magnetic Recording of Neuronal Activity . Neuron 95 , 1283 – 1291 .e4 ( 2017 ). OpenUrl CrossRef PubMed 16. ↵ Rosen , B. R. & Savoy , R. L. fMRI at 20: Has it changed the world? NeuroImage 62 , 1316 – 1324 ( 2012 ). OpenUrl CrossRef PubMed 17. ↵ Uhlhaas , P. J. et al. Magnetoencephalography as a Tool in Psychiatric Research: Current Status and Perspective . Biol. Psychiatry Cogn. Neurosci. Neuroimaging 2 , 235 – 244 ( 2017 ). OpenUrl PubMed 18. ↵ Leahy , R. M. , Mosher , J. C. , Spencer , M. E. , Huang , M. X. & Lewine , J. D. A study of dipole localization accuracy for MEG and EEG using a human skull phantom . Electroencephalogr. Clin. Neurophysiol . ( 1998 ). 19. ↵ Malmivuo , J. Comparison of the Properties of EEG and MEG in Detecting the Electric Activity of the Brain . Brain Topogr . 25 , 1 – 19 ( 2012 ). OpenUrl CrossRef PubMed 20. ↵ Logothetis , N. K. What we can do and what we cannot do with fMRI . Nature 453 , 869 – 878 ( 2008 ). OpenUrl CrossRef PubMed Web of Science 21. ↵ Lawrence , S. J. D. , Formisano , E. , Muckli , L. & De Lange , F. P. Laminar fMRI: Applications for cognitive neuroscience . NeuroImage 197 , 785 – 791 ( 2019 ). OpenUrl CrossRef PubMed 22. Self , M. W. , van Kerkoerle , T. , Goebel , R. & Roelfsema , P. R. Benchmarking laminar fMRI: Neuronal spiking and synaptic activity during top-down and bottom-up processing in the different layers of cortex . NeuroImage 197 , 806 – 817 ( 2019 ). OpenUrl CrossRef PubMed 23. ↵ Yang , J. , Huber , L. , Yu , Y. & Bandettini , P. A. Linking cortical circuit models to human cognition with laminar fMRI . Neurosci. Biobehav. Rev . 128 , 467 – 478 ( 2021 ). OpenUrl CrossRef PubMed 24. ↵ Samuelsson , J. G. , Peled , N. , Mamashli , F. , Ahveninen , J. & Hämäläinen , M. S. Spatial fidelity of MEG/EEG source estimates: A general evaluation approach . NeuroImage 224 , 117430 ( 2021 ). OpenUrl CrossRef PubMed 25. ↵ Troebinger , L. , López , J. D. , Lutti , A. , Bestmann , S. & Barnes , G. Discrimination of cortical laminae using MEG . NeuroImage 102 , 885 – 893 ( 2014 ). OpenUrl CrossRef PubMed 26. ↵ Bonaiuto , J. et al. Non-invasive laminar inference with MEG: Comparison of methods and source inversion algorithms . NeuroImage 167 , 372 – 383 ( 2018 ). OpenUrl CrossRef PubMed 27. ↵ Troebinger , L. et al. High precision anatomy for MEG . NeuroImage 86 , 583 – 591 ( 2014 ). OpenUrl CrossRef PubMed 28. ↵ Meyer , S. S. et al. Flexible head-casts for high spatial precision MEG . J. Neurosci. Methods 276 , 38 – 45 ( 2017 ). OpenUrl CrossRef PubMed 29. ↵ Szul , M. J. et al. Diverse beta burst waveform motifs characterize movement-related cortical dynamics . Prog. Neurobiol . 228 , 102490 ( 2023 ). OpenUrl CrossRef PubMed 30. ↵ Bonaiuto , J. et al. Lamina-specific cortical dynamics in human visual and sensorimotor cortices . eLife 7 , e33977 ( 2018 ). OpenUrl CrossRef PubMed 31. ↵ Bonaiuto , J. J. et al. Laminar dynamics of high amplitude beta bursts in human motor cortex . NeuroImage 242 , 118479 ( 2021 ). OpenUrl CrossRef PubMed 32. ↵ Fischl , B. FreeSurfer . NeuroImage 62 , 774 – 781 ( 2012 ). OpenUrl CrossRef PubMed Web of Science 33. ↵ Dale , A. M. , Fischl , B. & Sereno , M. I. Cortical Surface-Based Analysis: I. Segmentation and Surface Reconstruction . NeuroImage 9 , 179 – 194 ( 1999 ). OpenUrl CrossRef PubMed Web of Science 34. ↵ Bonaiuto , J. et al. Estimates of cortical column orientation improve MEG source inversion . NeuroImage 216 , 116862 ( 2020 ). OpenUrl CrossRef PubMed 35. ↵ Schroeder , W. , Martin , K. M. & Lorensen , W. E. The Visualization Toolkit an Object-Oriented Approach to 3D Graphics . ( Prentice-Hall, Inc. , 1998 ). 36. ↵ Ashburner , J. et al. SPM12 manual. Wellcome Trust Cent . Neuroimaging Lond. UK 2464 , ( 2014 ). 37. ↵ Nolte , G. The magnetic lead field theorem in the quasi-static approximation and its use for magnetoenchephalography forward calculation in realistic volume conductors . Phys. Med. Biol . 48 , 3637 – 3652 ( 2003 ). OpenUrl CrossRef PubMed Web of Science 38. ↵ Hillebrand , A. & Barnes , G. R. The use of anatomical constraints with MEG beamformers . NeuroImage 20 , 2302 – 13 ( 2003 ). OpenUrl CrossRef PubMed 39. ↵ Lin , F.-H. , Belliveau , J. W. , Dale , A. M. & Hämäläinen , M. S. Distributed current estimates using cortical orientation constraints . Hum. Brain Mapp . 27 , 1 – 13 ( 2006 ). OpenUrl CrossRef PubMed Web of Science 40. ↵ Brodmann , K. Vergleichende Lokalisationslehre Der Grosshirnrinde in Ihren Prinzipien Dargestellt Auf Grund Des Zellenbaues . ( Barth , 1909 ). 41. von Economo , C. F. & Koskinas , G. N. Die Cytoarchitektonik Der Hirnrinde Des Erwachsenen Menschen . ( J. Springer , 1925 ). 42. ↵ Wagstyl , K. , Ronan , L. , Goodyer , I. M. & Fletcher , P. C. Cortical thickness gradients in structural hierarchies . Neuroimage 111 , 241 – 250 ( 2015 ). OpenUrl CrossRef PubMed 43. ↵ Friston , K. et al. Multiple sparse priors for the M/EEG inverse problem . NeuroImage 39 , 1104 – 1120 ( 2008 ). OpenUrl CrossRef PubMed Web of Science 44. ↵ Wagstyl , K. et al. BigBrain 3D atlas of cortical layers: Cortical and laminar thickness gradients diverge in sensory and motor cortices . PLoS Biol . 18 , e3000678 ( 2020 ). OpenUrl CrossRef PubMed 45. ↵ Wagstyl , K. et al. Mapping Cortical Laminar Structure in the 3D BigBrain . Cereb. Cortex 28 , 2551 – 2562 ( 2018 ). OpenUrl CrossRef PubMed 46. ↵ Fischl , B. , Sereno , M. I. , Tootell , R. B. H. & Dale , A. M. High-resolution intersubject averaging and a coordinate system for the cortical surface . Hum. Brain Mapp . 8 , 272 – 284 ( 1999 ). OpenUrl CrossRef PubMed Web of Science 47. ↵ Daunizeau , J. , Adam , V. & Rigoux , L. VBA: a probabilistic treatment of nonlinear models for neurobiological and behavioural data . PLoS Comput. Biol . 10 , e1003441 ( 2014 ). OpenUrl CrossRef PubMed 48. ↵ Stephan , K. E. , Penny , W. D. , Daunizeau , J. , Moran , R. J. & Friston , K. J. Bayesian model selection for group studies . NeuroImage 46 , 1004 – 1017 ( 2009 ). OpenUrl CrossRef PubMed Web of Science 49. ↵ Belardinelli , P. , Ortiz , E. , Barnes , G. R. , Noppeney , U. & Preissl , H. Source reconstruction accuracy of MEG and EEG Bayesian inversion approaches . PloS One 7 , e51985 ( 2012 ). OpenUrl CrossRef PubMed 50. ↵ Kinsbourne , M. Hemineglect and hemisphere rivalry . Adv Neurol 18 , 41 – 49 ( 1977 ). OpenUrl PubMed 51. ↵ Yabuta , N. H. & Callaway , E. M. Functional streams and local connections of layer 4C neurons in primary visual cortex of the macaque monkey . J. Neurosci . 18 , 9489 – 9499 ( 1998 ). OpenUrl Abstract / FREE Full Text 52. ↵ Yoshimura , Y. , Dantzker , J. L. & Callaway , E. M. Excitatory cortical neurons form fine-scale functional networks . Nature 433 , 868 – 873 ( 2005 ). OpenUrl CrossRef PubMed Web of Science 53. ↵ Bortone , D. S. , Olsen , S. R. & Scanziani , M. Translaminar inhibitory cells recruited by layer 6 corticothalamic neurons suppress visual cortex . Neuron 82 , 474 – 485 ( 2014 ). OpenUrl CrossRef PubMed Web of Science 54. ↵ Callaway , E. M. Feedforward, feedback and inhibitory connections in primate visual cortex . Neural Netw . 17 , 625 – 632 ( 2004 ). OpenUrl CrossRef PubMed Web of Science 55. Self , M. W. , van Kerkoerle , T. , Super , H. & Roelfsema , P. R. Distinct roles of the cortical layers of area V1 in figure-ground segregation . Curr. Biol . 23 , 2121 – 2129 ( 2013 ). OpenUrl CrossRef PubMed 56. ↵ Maier , A. , Adams , G. K. , Aura , C. & Leopold , D. A. Distinct superficial and deep laminar domains of activity in the visual cortex during rest and stimulation . Front. Syst. Neurosci . 4 , 31 ( 2010 ). OpenUrl CrossRef PubMed 57. ↵ Cheyne , D. , Bakhtazad , L. & Gaetz , W. Spatiotemporal mapping of cortical activity accompanying voluntary movements using an event-related beamforming approach . Hum. Brain Mapp . 27 , 213 – 229 ( 2006 ). OpenUrl CrossRef PubMed Web of Science 58. Cheyne , D. & Weinberg , H. Neuromagnetic fields accompanying unilateral finger movements: pre-movement and movement-evoked fields . Exp. Brain Res . 78 , ( 1989 ). 59. ↵ Kristeva , R. , Cheyne , D. & Deecke , L. Neuromagnetic fields accompanying unilateral and bilateral voluntary movements: topography and analysis of cortical sources . Electroencephalogr. Clin. Neurophysiol. Potentials Sect . 81 , 284 – 298 ( 1991 ). OpenUrl 60. ↵ Oishi , M. , Kameyama , S. , Fukuda , M. , Tsuchiya , K. & Kondo , T. Cortical activation in area 3b related to finger movement: an MEG study . Neuroreport 15 , 57 – 62 ( 2004 ). OpenUrl CrossRef PubMed Web of Science 61. ↵ Brookes , M. J. et al. Magnetoencephalography with optically pumped magnetometers (OPM-MEG): the next generation of functional neuroimaging . Trends Neurosci . 45 , 621 – 634 ( 2022 ). OpenUrl CrossRef PubMed 62. ↵ Hinne , M. , Gronau , Q. F. , van den Bergh , D. & Wagenmakers , E.-J. A Conceptual Introduction to Bayesian Model Averaging . Adv. Methods Pract. Psychol. Sci . 3 , 200 – 215 ( 2020 ). OpenUrl CrossRef 63. ↵ Fernandez Pujol , C. , Blundon , E. G. & Dykstra , A. R. Laminar specificity of the auditory perceptual awareness negativity: A biophysical modeling study . PLoS Comput. Biol . 19 , e1011003 ( 2023 ). OpenUrl CrossRef PubMed 64. Neymotin , S. A. et al. Human Neocortical Neurosolver (HNN), a new software tool for interpreting the cellular and network origin of human MEG/EEG data . eLife 9 , e51214 ( 2020 ). OpenUrl CrossRef PubMed 65. ↵ Fernandez Pujol , C. , Bruce , J. , Thorpe , R. V. , Jones , S. R. & Dykstra , A. R. Neurophysiology of mismatch negativity generation: a biophysical modeling study . bioRxiv 2025 – 11 ( 2025 ). 66. ↵ Buzsáki , G. , Anastassiou , C. A. & Koch , C. The origin of extracellular fields and currents — EEG, ECoG, LFP and spikes . Nat. Rev. Neurosci . 13 , 407 ( 2012 ). OpenUrl CrossRef PubMed 67. ↵ Mitzdorf , U. Current source-density method and application in cat cerebral cortex: investigation of evoked potentials and EEG phenomena . Physiol. Rev . 65 , 37 – 100 ( 1985 ). OpenUrl CrossRef PubMed Web of Science 68. ↵ Tenke , C. E. & Kayser , J. Generator localization by current source density (CSD): Implications of volume conduction and field closure at intracranial and scalp resolutions . Clin. Neurophysiol . 123 , 2328 – 2345 ( 2012 ). OpenUrl CrossRef PubMed 69. ↵ Vidaurre , D. et al. Spectrally resolved fast transient brain states in electrophysiological data . NeuroImage 126 , 81 – 95 ( 2016 ). OpenUrl CrossRef PubMed 70. ↵ Klein , N. , Siegle , J. H. , Teichert , T. & Kass , R. E. Cross-population coupling of neural activity based on Gaussian process current source densities . PLOS Comput. Biol . 17 , e1009601 ( 2021 ). OpenUrl CrossRef PubMed 71. ↵ Wang , W. et al. Design of locally arranged sensor arrays in wearable OPM-MEG based on sensor volume constraints . Measurement 229 , 114373 ( 2024 ). OpenUrl CrossRef 72. ↵ Stolk , A. , Todorovic , A. , Schoffelen , J.-M. & Oostenveld , R. Online and offline tools for head movement compensation in MEG . NeuroImage 68 , 39 – 48 ( 2013 ). OpenUrl CrossRef PubMed Web of Science 73. ↵ Taulu , S. & Kajola , M. Presentation of electromagnetic multichannel data: The signal space separation method . J. Appl. Phys . 97 , 124905 ( 2005 ). OpenUrl CrossRef 74. ↵ Taulu , S. & Simola , J. Spatiotemporal signal space separation method for rejecting nearby interference in MEG measurements . Phys. Med. Biol . 51 , 1759 ( 2006 ). OpenUrl CrossRef PubMed 75. ↵ Murakami , S. & Okada , Y. Contributions of principal neocortical neurons to magnetoencephalography and electroencephalography signals . J. Physiol . 575 , 925 – 936 ( 2006 ). OpenUrl CrossRef PubMed Web of Science 76. ↵ Ribeiro , P. F. M. et al. The human cerebral cortex is neither one nor many: neuronal distribution reveals two quantitatively different zones in the gray matter, three in the white matter, and explains local variations in cortical folding . Front. Neuroanat . 7 , ( 2013 ). 77. ↵ Cichy , R. M. , Ramirez , F. M. & Pantazis , D. Can visual information encoded in cortical columns be decoded from magnetoencephalography data in humans? NeuroImage 121 , 193 – 204 ( 2015 ). OpenUrl CrossRef PubMed 78. ↵ Nasiotis , K. , Clavagnier , S. , Baillet , S. & Pack , C. C. High-resolution retinotopic maps estimated with magnetoencephalography . NeuroImage 145 , 107 – 117 ( 2017 ). OpenUrl CrossRef PubMed 79. ↵ Barratt , E. L. , Francis , S. T. , Morris , P. G. & Brookes , M. J. Mapping the topological organisation of beta oscillations in motor cortex using MEG . NeuroImage 181 , 831 – 844 ( 2018 ). OpenUrl CrossRef PubMed 80. ↵ Bok , S. Der Einfluss in den Furchen und Windungen auftretenden Krümmungen der Grosshirnrinde auf die Rindenarchitektur . Z. Für Gesamte Neurol. Psychiatr . 121 , 682 – 750 ( 1929 ). OpenUrl CrossRef 81. ↵ Leuze , C. W. U. et al. Layer-Specific Intracortical Connectivity Revealed with Diffusion MRI . Cereb. Cortex 24 , 328 – 339 ( 2014 ). OpenUrl CrossRef PubMed Web of Science 82. ↵ Gulban , O. F. et al. Cortical fibers orientation mapping using in-vivo whole brain 7 T diffusion MRI . Neuroimage 178 , 104 – 118 ( 2018 ). OpenUrl CrossRef PubMed 83. ↵ St-Onge , E. , Daducci , A. , Girard , G. & Descoteaux , M. Surface-enhanced tractography (SET) . NeuroImage 169 , 524 – 539 ( 2018 ). OpenUrl CrossRef PubMed 84. ↵ Havlicek , M. & Uludağ , K. A dynamical model of the laminar BOLD response . NeuroImage 204 , 116209 ( 2020 ). OpenUrl CrossRef PubMed 85. ↵ Markuerkiaga , I. , Marques , J. P. , Gallagher , T. E. & Norris , D. G. Estimation of laminar BOLD activation profiles using deconvolution with a physiological point spread function . J. Neurosci. Methods 353 , 109095 ( 2021 ). OpenUrl CrossRef PubMed 86. ↵ Chai , Y. , Li , L. , Huber , L. , Poser , B. A. & Bandettini , P. A. Integrated VASO and perfusion contrast: A new tool for laminar functional MRI . NeuroImage 207 , 116358 ( 2020 ). OpenUrl CrossRef PubMed 87. ↵ Polimeni , J. R. , Fischl , B. , Greve , D. N. & Wald , L. L. Laminar analysis of 7 T BOLD using an imposed spatial activation pattern in human V1 . NeuroImage 52 , 1334 – 1346 ( 2010 ). OpenUrl CrossRef PubMed Web of Science 88. ↵ Chai , Y. et al. Unlocking near-whole-brain, layer-specific functional connectivity with 3D VAPER fMRI . Imaging Neurosci . 2 , 1 – 20 ( 2024 ). OpenUrl 89. ↵ Rayson , H. et al. Bursting with Potential: How Sensorimotor Beta Bursts Develop from Infancy to Adulthood . J. Neurosci . 43 , 8487 – 8503 ( 2023 ). OpenUrl Abstract / FREE Full Text 90. ↵ Chen , G. , Wang , F. , Gore , J. C. & Roe , A. W. Identification of cortical lamination in awake monkeys by high resolution magnetic resonance imaging . NeuroImage 59 , 3441 – 3449 ( 2012 ). OpenUrl CrossRef PubMed 91. Trampel , R. , Bazin , P.-L. , Pine , K. & Weiskopf , N. In-vivo magnetic resonance imaging (MRI) of laminae in the human cortex . NeuroImage 197 , 707 – 715 ( 2019 ). OpenUrl CrossRef PubMed 92. Kundu , S. et al. Mapping the individual human cortex using multidimensional MRI and unsupervised learning . Brain Commun . 5 , fcad258 ( 2023 ). OpenUrl 93. Autio , J. A. et al. Charting cortical-layer specific area boundaries using Gibbs’ ringing attenuated T1w/T2w-FLAIR myelin MRI . bioRxiv 2024.09.27.615294 ( 2024 ) doi: 10.1101/2024.09.27.615294 . OpenUrl Abstract / FREE Full Text 94. ↵ Zeng , X. et al. Segmentation of supragranular and infragranular layers in ultra-high-resolution 7T ex vivo MRI of the human cerebral cortex . Cereb. Cortex 34 , bhae362 ( 2024 ). OpenUrl CrossRef PubMed 95. ↵ Stenroos , M. , Hunold , A. & Haueisen , J. Comparison of three-shell and simplified volume conductor models in magnetoencephalography . NeuroImage 94 , 337 – 348 ( 2014 ). OpenUrl CrossRef PubMed 96. ↵ Erdbrügger , T. et al. CutFEM-based MEG forward modeling improves source separability and sensitivity to quasi-radial sources: A somatosensory group study . Hum. Brain Mapp . 45 , e26810 ( 2024 ). OpenUrl CrossRef PubMed 97. ↵ Medani , T. et al. Brainstorm-DUNEuro: An integrated and user-friendly Finite Element Method for modeling electromagnetic brain activity . NeuroImage 267 , 119851 ( 2023 ). OpenUrl CrossRef PubMed 98. ↵ Helbling , S. Inferring laminar origins of MEG signals with optically pumped magnetometers (OPMs): A simulation study . Imaging Neurosci . 3 , imag_a_00410 ( 2025 ). OpenUrl References 1. ↵ Oldfield , R. C. The assessment and analysis of handedness: The Edinburgh inventory . Neuropsychologia 9 , 97 – 113 ( 1971 ). OpenUrl CrossRef PubMed Web of Science 2. ↵ Fischl , B. FreeSurfer . NeuroImage 62 , 774 – 781 ( 2012 ). OpenUrl CrossRef PubMed Web of Science 3. ↵ Meyer , S. S. et al. Flexible head-casts for high spatial precision MEG . J. Neurosci. Methods 276 , 38 – 45 ( 2017 ). OpenUrl CrossRef PubMed 4. ↵ Weiskopf , N. et al. Quantitative multi-parameter mapping of R1, PD(*), MT, and R2(*) at 3T: a multi-center validation . Front. Neurosci . 7 , 95 ( 2013 ). OpenUrl CrossRef PubMed 5. ↵ Tabelow , K. et al. hMRI–A toolbox for quantitative MRI in neuroscience and clinical research . Neuroimage 194 , 191 – 210 ( 2019 ). OpenUrl CrossRef PubMed 6. ↵ Carey , D. , et al. Quantitative MRI Provides Markers Of Intra-, Inter-Regional, And Age-Related Differences In Young Adult Cortical Microstructure . bioRxiv doi: 10.1101/139568 ( 2017 ) 10.1101/139568. OpenUrl Abstract / FREE Full Text 7. ↵ Gramfort , A. et al. MEG and EEG data analysis with MNE-Python . Front. Neuroinformatics 7 , 267 ( 2013 ). OpenUrl 8. ↵ Geyer , S. Turner , R. Fischl , B. Estimating the Location of Brodmann Areas from Cortical Folding Patterns Using Histology and Ex Vivo MRI . in Microstructural Parcellation of the Human Cerebral Cortex (eds Geyer , S. & Turner , R. ) 129 – 156 ( Springer Berlin Heidelberg , Berlin, Heidelberg , 2013 ). doi: 10.1007/978-3-642-37824-9_4 . OpenUrl CrossRef 9. ↵ Bijanzadeh , M. , Nurminen , L. , Merlin , S. , Clark , A. M. & Angelucci , A. Distinct laminar processing of local and global context in primate primary visual cortex . Neuron 100 , 259 – 274 ( 2018 ). OpenUrl CrossRef PubMed 10. ↵ Troebinger , L. et al. High precision anatomy for MEG . NeuroImage 86 , 583 – 591 ( 2014 ). OpenUrl CrossRef PubMed 11. ↵ Bonaiuto , J. et al. Lamina-specific cortical dynamics in human visual and sensorimotor cortices . eLife 7 , e33977 ( 2018 ). OpenUrl CrossRef PubMed 12. ↵ Bonaiuto , J. et al. Non-invasive laminar inference with MEG: Comparison of methods and source inversion algorithms . NeuroImage 167 , 372 – 383 ( 2018 ). OpenUrl CrossRef PubMed 13. ↵ Taulu , S. & Simola , J. Spatiotemporal signal space separation method for rejecting nearby interference in MEG measurements . Phys. Med. Biol . 51 , 1759 ( 2006 ). OpenUrl CrossRef PubMed 14. ↵ Taulu , S. , Kajola , M. & Simola , J. Suppression of Interference and Artifacts by the Signal Space Separation Method . Brain Topogr . 16 , 269 – 275 ( 2004 ). OpenUrl CrossRef PubMed Web of Science 15. ↵ de Cheveigné , A. ZapLine: A simple and effective method to remove power line artifacts . NeuroImage 207 , 116356 ( 2020 ). OpenUrl CrossRef PubMed 16. ↵ Dykstra , A. R. & Gutschalk , A. Does the mismatch negativity operate on a consciously accessible memory trace? Sci. Adv . 1 , e1500677 ( 2015 ). OpenUrl FREE Full Text 17. ↵ Jas , M. , Engemann , D. A. , Bekhti , Y. , Raimondo , F. & Gramfort , A. Autoreject: Automated artifact rejection for MEG and EEG data . NeuroImage 159 , 417 – 429 ( 2017 ). OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted February 07, 2026. Download PDF Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Multilayer MEG source modelling enables depth-resolved inference across all six cortical laminae in humans Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Multilayer MEG source modelling enables depth-resolved inference across all six cortical laminae in humans Maciej J. Szul , Ishita Agarwal , Quentin Moreau , Solène Gailhard , Carolina Fernandez Pujol , Yunkai Zhu , Matteo Maspoli , Danila Shelepenkov , Bassem Hiba , Sebastien Daligault , Franck Lamberton , Denis Schwartz , Mathilde Bonnefond , Andrew R. Dykstra , Sven Bestmann , Gareth R. Barnes , James J. Bonaiuto bioRxiv 2025.05.28.656642; doi: https://doi.org/10.1101/2025.05.28.656642 Share This Article: Copy Citation Tools Multilayer MEG source modelling enables depth-resolved inference across all six cortical laminae in humans Maciej J. Szul , Ishita Agarwal , Quentin Moreau , Solène Gailhard , Carolina Fernandez Pujol , Yunkai Zhu , Matteo Maspoli , Danila Shelepenkov , Bassem Hiba , Sebastien Daligault , Franck Lamberton , Denis Schwartz , Mathilde Bonnefond , Andrew R. Dykstra , Sven Bestmann , Gareth R. Barnes , James J. Bonaiuto bioRxiv 2025.05.28.656642; doi: https://doi.org/10.1101/2025.05.28.656642 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Neuroscience Subject Areas All Articles Animal Behavior and Cognition (7616) Biochemistry (17625) Bioengineering (13852) Bioinformatics (41825) Biophysics (21397) Cancer Biology (18524) Cell Biology (25417) Clinical Trials (138) Developmental Biology (13350) Ecology (19858) Epidemiology (2067) Evolutionary Biology (24277) Genetics (15581) Genomics (22459) Immunology (17698) Microbiology (40278) Molecular Biology (17134) Neuroscience (88400) Paleontology (666) Pathology (2823) Pharmacology and Toxicology (4812) Physiology (7632) Plant Biology (15106) Scientific Communication and Education (2042) Synthetic Biology (4281) Systems Biology (9807) Zoology (2266)