Resolving the mesoscopic missing link: Biophysical modeling of EEG from cortical columns in primates

Event-related potentials (ERP) are among the most widely measured indices for studying human cognition. While their timing and magnitude provide valuable insights, their usefulness is limited by our understanding of their neural generators at the circuit level. Inverse source localization offers insights into such generators, but their solutions are not unique. To address this problem, scientists have assumed the source space generating such signals comprises a set of discrete equivalent current dipoles, representing the activity of small cortical regions. Based on this notion, theoretical studies have employed forward modeling of scalp potentials to understand how changes in circuit-level dynamics translate into macroscopic ERPs. However, experimental validation is lacking because it requires in vivo measurements of intracranial brain sources. Laminar local field potentials (LFP) offer a mechanism for estimating intracranial current sources. Yet, a theoretical link between LFPs and intracranial brain sources is missing. Here, we present a forward modeling approach for estimating mesoscopic intracranial brain sources from LFPs and predict their contribution to macroscopic ERPs. We evaluate the accuracy of this LFP-based representation of brain sources utilizing synthetic laminar neurophysiological measurements and then demonstrate the power of the approach in vivo to clarify the source of a representative cognitive ERP component. To that end, LFP was measured across the cortical layers of visual area V4 in macaque monkeys performing an attention demanding task. We show that area V4 generates dipoles through layer-specific transsynaptic currents that biophysically recapitulate the ERP component through the detailed forward modeling. The constraints imposed on EEG production by this method also revealed an important dissociation between computational and biophysical contributors. As such, this approach represents an important bridge between laminar microcircuitry, through the mesoscopic activity of cortical columns to the patterns of EEG we measure at the scalp.


Introduction
Identifying the neural sources of EEG is a substantial challenge facing the human neuroscience community. Many basic and clinical research programs rely on EEG as a window into human sensation, perception, cognition, and action. However, some insights into the mechanisms of these neural operations are invisible through this method. While the macroscopic EEG contains information, some of it is lost relative to the microscopic activations it reflects. As such, researchers are developing methods to extract as much information regarding the neural sources of EEG to strengthen this staple of neuroscientific approaches. EEG inverse source localization methods have offered a tool to estimate plausible brain sources generating such signals and, hence, make inferences about the underlying cortical mechanisms (Grech et al., 2008;Michel et al., 2004). Unfortunately, the results are indefinite because a given cranial polarization pattern can arise from multiple current configurations, requiring additional assumptions (Grech et al., 2008;Nunez et al., 2019). Researchers constrain the number of possible brain source configurations by imposing biophysical models to address this problem. The most widely used models are based on the equivalent current dipole, with each dipole representing the activity of a small cortical region (Nunez et al., 2019). The equivalent dipole model has been implemented in different inverse solutions, from a few discrete dipoles (Scherg et al., 2019) to more distributed inverse solution algorithms (e.g., sLORETA, MNE, MxNE) (Gramfort et al., 2012;Grech et al., 2008;Nunez et al., 2019).
Despite this progress, the biophysical relationship between the mesoscopic laminar LFP/CSD and the microscopic cellular sources has not been evaluated. This gap of knowledge renders the validity of predicting spatial patterns of cranial EEG voltage from current sources derived from laminar LFP questionable. Hence, we describe a biophysically plausible forward modeling approach linking mesoscopic CSD derived from LFP to macroscopic ERPs.
Based on earlier demonstrations that macaque monkeys produce homologues of human cognitive ERP components (Cohen et al., 2009;Heitz et al., 2010;Purcell et al., 2013;Reinhart et al., 2012;Sajad et al., 2019;Westerberg et al., 2020b;Woodman et al., 2007), it is now possible to obtain data necessary to measure current dipoles in the cerebral cortex under conditions generating cognitive ERPs. Recently, Westerberg and colleagues (Westerberg et al., 2022) reported current source density maps of the laminar distributions of transmembrane currents in extrastriate cortical area V4 associated with the cognitive ERP component known as the N2pc. We use these data to demonstrate the utility of our detailed forward modeling approach linking intracranial patterns of laminar activity to the overlying EEG. We show that V4 generates dipoles through layer-specific transsynaptic currents that biophysically contribute to ERP generation. Forward modeling this cortical activity renders EEG at the scalp consistent with that of previous reports of the representative ERP. Moreover, in evaluating the potential contributions of other cortical areas computationally involved in the ERP-indexed operation, this approach revealed that these computational contributors need not biophysically contribute to EEG production. In establishing a mesoscopic link between microscopic neural currents and macroscopic cranial voltages, this finding represents the first, definite forward model of a cognitive ERP from current dipoles derived from neural activity in primates. ketamine (5-25 mg/kg). Monkeys were then catheterized and intubated. Surgeries were conducted aseptically with animals under isoflurane (1-5%) anesthesia. EKG, temperature, and respiration were monitored. Postoperative antibiotics and analgesics were administered. Further detail is documented elsewhere (Westerberg et al., 2020a(Westerberg et al., , 2020b.

Magnetic resonance imaging
Anesthetized animals were placed in a 3 Tesla Magnetic Resonance Imaging (MRI) scanner (Phillips) at the Vanderbilt University Institute of Imaging Science. T1-weighted 3D MPRAGE scans were acquired with a 32-channel head coil equipped for SENSE imaging. Images were acquired using 0.5 mm isotropic voxel resolution with parameters: repetition 5 s, echo 2.5 ms, flip angle 7°2 .

Cognitive task -visual search
Monkeys performed a color pop-out search. Search arrays were presented on a CRT monitor at 60 Hz, at 57 cm distance. Stimulus generation and timing were done with TEMPO (Reflective Computing). Event times were assessed with a photodiode on the CRT. We used isoluminant red and green disks on a gray background. Target feature varied within a session. Trials were initiated by fixating within 0.5° of visual angle (dva) of a fixation dot. The time between fixation and array onset was between 750 and 1250 ms. A nonaging foreperiod function was used to determine the fixation period on a trial-by-trial basis. Arrays comprised of 6 items. Array item size scaled with eccentricity at 0.3 dva per 1 dva eccentricity so that they were smaller than the average V4 receptive field (RF) (Freeman and Simoncelli, 2011). The angular position of items relative to fixation varied session to session so that 1 item was positioned at the center of the RF. Items were equally spaced relative to each other and located at the same eccentricity. In each trial, one array item was different from the others. Monkeys made a saccade to the oddball within 1 second and maintained fixation within 2-5 dva of the target for 500 ms. Juice reward was administered following the successful completion of the trial. The target item had an equal probability of being located at any of the 6 locations. Eye movements were monitored at 1 kHz using a corneal reflection system (SR Research Eyelink). If the monkey failed to saccade to the target, they experienced a timeout (1-5 s).

Laminar CSD Recording Procedure
Laminar V4 neurophysiology was acquired at 24 kHz using a PZ5 and RZ2 (Tucker-Davis). Signals were filtered between 0.1 Hz -12 kHz. V4 data was collected from 2 monkeys (monkey Ca: left hemisphere; He: right) across 30 sessions (monkey Ca: 21; monkey He: 9) using 32-channel linear electrode arrays with 0.1 mm interelectrode spacing (Plexon) introduced through the intact dura mater each session. Arrays spanned layers of V4 with a subset of electrode contacts deliberately left outside of cortex. We computed the CSD signal using the spline-iCSD method (Pettersen et al., 2006) as implemented in the CSDplotter toolbox (https://github.com/espenhgn/CSDplotter) with custom MATLAB (R2021b, The MathWorks) scripts. In Fig. 1, we compare our estimated CSD maps with those reported by Westerberg and colleagues (Westerberg et al., 2022) using the standard CSD method (Nicholson and Freeman, 1975). For the standard CSD method, CSD was computed from the raw signal by taking the second spatial derivative along electrodes (Mehta et al., 2000;Nicholson and Freeman, 1975;Schroeder et al., 1998;Westerberg et al., 2019) and converting voltage to current (Logothetis et al., 2007). We computed the CSD by taking the second spatial derivative of the LFP: where x is the extracellular voltage at time t measured at an electrode contact at depth d, z is the inter-electrode distance, and σ is conductivity. CSD was baseline corrected at the trial level by subtracting the average activation during the 300 ms preceding array onset from the response at all timepoints. CSD was clipped 10 ms before saccade at the trial level to eliminate the influence of eye movements.

Laminar Alignment
Orthogonal array penetrations were confirmed online through a reverse-correlation RF mapping procedure Nandy et al., 2017;Westerberg et al., 2019). RFs were found to represent portions of visual space consistent with previous reports of V4 (Gattass et al., 1988). An expanded description of the RF mapping procedure for this dataset has been reported previously (Westerberg et al., 2022. Positions of recording sites relative to V4 layers were determined using CSD (Nandy et al., 2017;Schroeder et al., 1998). Current sinks following visual stimulation first appear in the granular input layers of the cortex, then propagate to extragranular compartments. We computed CSD and identified the granular input sink session-wise. Sessions were aligned by this input sink. 'L4' refers to the granular input layer, 'L2/3' to the supragranular layers, and 'L5/6' to the infragranular layers.

Boundary Element Model
The monkey's head was modeled as an isotropic and piecewise homogenous volume conductor comprised of the scalp, inner and outer skull, and the cortex surface. For the forward modeling of the experimental data, we employed the surfaces of the cortex, and the scalp and skull compartments obtained from the segmentation of the T1-weighted MRI of monkey Y in SPM12 (Penny et al., 2011) and Brain-Suite (Shattuck et al., 2001), respectively. In the detailed biophysical simulations, we utilized the symmetric surfaces provided in the NIMH Macaque Template version 2.0 (Jung et al., 2021) to construct the head model ( Fig. 2A). The position of the EEG electrodes for both models was defined employing the EEG 10-10 system and the monkey's scalp surface as described elsewhere (Giacometti et al., 2014). The scalp, skull, and brain conductivities were set as 0.43, 0.0063, and 0.33 S/m (Lee et al., 2015), respectively.

EEG Forward Model
The EEG scalp potential V e ( r e , t) at any position r e evoked by a continuous field of microscopic electric currents I( r , t) inside the brain R can be represented by the following inhomogeneous Fredholm integral equation of the second kind (Eq. (2)) (Riera et al., 2012): V e r e , t = V 0 r e , t + with j k (I, r ) = (σ k + 1 − σ k )v k (I, r ) n k ( r )/Δ1 representing the secondary currents defined for each elemental volumetric shell Ω k (i.e., a surface S k of thickness Δl → 0). σ k and v k (I, r ) denote the conductivity and surface potential of the k-th compartment in the head model (i.e., brain (σ b ), skull, and scalp)., and n k ( r ) the normal vector to the surface (S k ) of the k-th compartment at location r .
We assume I( r , t) = s( r , t) for r ∈ V and I( r , t) = 0 otherwise, where V is the volume of the brain region of interest, centered at r m . If the location of the EEG electrode ( r e ) is far enough from the center r m , then V 0 ( r e , t) can be calculated as a function of the multipolar moments (Riera et al., 2012). Under the assumption of the dipolar model, the EEG forward model can be represented by Eq. (2a) and the following equation for V 0 ( r e , t): V 0 ( r e , t) = 1 where d ( r , t) denotes the current dipole moment. The theoretical framework and numerical strategies used to compute the potentials (v k (I, r ), V 0 ( r e , t)) to obtain the EEG scalp potential V e ( r e , t) are detailed elsewhere (Hämäläinen and Sarvas, 1989).
The current dipole moment d ( r , t) is assumed to originate from post-synaptic currents caused mainly by the activation of pyramidal cells perpendicular to the cortical surface (Hämäläinen et al., 1993) and from nonlinear processes taking place in their dendrites (Herrera et al., 2020;Suzuki and Larkum, 2017). Therefore, we can write d ( r , t) as μ ( r ) ⋅ d(t), with μ ( r ) and d(t) representing its orientation and time-varying amplitude, respectively (Riera et al., 2012). The activation waveform d(t) was estimated from the volumetric current sources (CSD) calculated from the laminar recordings, and the orientation of the dipoles from a weighted average of the normal to the cortical surface at the location of the dipole and the normal of its neighboring triangles. Considering the assumptions made to calculate the CSD, the amplitude d(t) of the current dipole moment from the CSD in a volume of interest V is given by (Riera et al., 2012): The z-axis is defined as the direction perpendicular to the cortex with positive and negative values toward the supragranular and infragranular layers, respectively. We modeled dipoles located in areas V4, lusV4, LIP, 7a, and FEF and assumed that the cortical columns in those areas were perfect cylinders of 3 mm radius (r c ). The integrals were calculated using the trapezoidal method, where each subinterval corresponds to a particular grid point in the corresponding CSD method.

Biophysical Simulations
We generated a synthetic data set using detailed biophysical models of neurons to evaluate the equivalent current dipole derived from the CSD. We simulated the responses to a current pulse of 30 ms duration with noisy random amplitude in two populations of independent (unconnected) L3 and L5 pyramidal cells described previously by (Eyal et al., 2018) (ModelDB, accession #238347, 2013_03_06_cell03_789_H41_03, active model cell0603_08_model_602) and (Hay et al., 2011) (ModelDB, accession #139653, "cell #1"). The somas of the neurons were distributed uniformly in a 3 mm diameter cylinder with height corresponding to the vertical extent in area V4 of lower L3 (675-750 μm below the pia matter) and of L5 (1250-1750 μm) (Fig. 2C). We determined the number of simulated neurons from each population, maintaining the ratio of neurons in L3 and L5 reported in the anatomical literature (Vanni et al., 2020) and finding the minimal number of neurons in a cortical column of 3 mm diameter under this stimulation paradigm necessary to reproduce the approximate amplitude of the CSD observed in the V4 experimental data (Fig. 1D). We considered this approach since not all neurons are active at the same time during a task in a real cortical column composed of millions of neurons (e.g., see (Westerberg et al., 2022)). We simulated 2,200 L3 and 1000 L5 pyramidal cells. The amplitude of the current pulse was sampled from a Gaussian distribution with a standard deviation of 0.3 nA. To simulate the asymmetry observed in the experimental recordings from V4, we used mean current amplitudes of 1.90 nA and 1.85 nA for neurons in the left and right hemispheres, respectively. Reducing the mean current amplitude for neurons in the right hemisphere decreased the probability of eliciting dendritic Ca 2+ -spikes in L5 pyramidal cells, reducing the amplitude of the late sink/sources in the CSD map ( Fig. 2E) and creating a weaker dipole current source. The beginning of stimulation of each neuron was randomly sampled from a uniform distribution between 10 and 20 ms after the simulation was initiated. All simulations were performed in Python using LFPy (Hagen et al., 2018), which builds on NEURON (Hines et al., 2009).
Simulating measurements with a linear microelectrode array, we calculated the LFP produced by the activity of the neurons at 17 equally spaced vertically aligned points located at the center of the cortical column. As in the experiments, the inter-electrode distance was 100 μm. We employed the point-source approximation in LFPy to compute the extracellular potentials. The LFP was obtained by low-pass filtering at 100 Hz. The CSD patterns of the synthetic data sets were calculated using the spline-iCSD method (Pettersen et al., 2006) with the custom MATLAB (R2021b, The MathWorks) scripts used for the experimental data.
We estimated the EEG produced by the simulated population of pyramidal neurons using three approaches: compartment-based, summed transmembrane currents single-dipole, and CSD single-dipole approaches. The compartment-based approach considers the transmembrane currents entering and escaping the extracellular medium as current sources/ sinks. It calculates the EEG mapping the transmembrane currents of all neurons into the scalp potentials (Herrera et al., 2020). Thus, for the compartment-based approach, Eq. (2b) for V 0 ( r e , t) becomes (Herrera et al., 2020): summed transmembrane currents single-dipole and CSD single-dipole approaches, the EEG is estimated using the EEG dipolar model explain in section 2.7 (Eq. (2a) and Eq. (3)), which assumes the activity of a cortical column can be accurately represented by an equivalent current dipole moment at the center of the column. The summed transmembrane currents single-dipole approach calculates the equivalent current dipole moment as the sum of all dipoles between the center of mass or center of the cortical column and every compartment of each neuron (Eq. (6)). Thus, the amplitude and orientation of the equivalent current dipole moment are given by the sum (Hagen et al., 2018;Naess et al., 2021): with r m representing the center of mass of the cortical column, which is the coordinates of the dipole. In the CSD single-dipole approach, the orientation of the dipoles is fixed, corresponding to that of the cortical column, and the amplitude is calculated from Eq. (4). For our simulations, N n = 2, N p = {1,2} with N 1 = N L3 = 2200 and N 1 = N L5 = 1000. The estimated current dipole moments for each neuronal population and their combined activity were low-pass filtered at 100 Hz.
To validate the EEG estimated using the CSD single-dipole (Eq. (4)) and the transmembrane currents single-dipole (Eq. (6)) approaches, we compared them to the compartment-based EEG, considered as the "ground-truth" EEG (Eq. (5)). We quantified the mismatch between the transmembrane currents single-dipole and the CSD single-dipole estimations (y) and the ground-truth EEG (y) by calculating the relative magnitude (MAG) (Eq. (7)) and the relative difference (RDM) (Eq. (8)) (Meijs et al., 1989): where a relative difference (RDM) closer to 0 and magnitude (MAG) closer to 1 indicates a higher degree of similarity between a dipolar approach and the ground truth. To account for different depths and orientations of the cortical column, we calculated these measures at 15 randomly selected cortical columns throughout the monkey brain. Furthermore, given the relationship between electric field strength relative to dipole distance, we evaluated the dependence of these measures on the depth of the cortical column relative to the scalp. We defined the depth of the cortical column as the shortest distance between the center of mass of the cortical column and the scalp surface utilized in the EEG forward model. We calculated this distance using the Matlab point2trimesh function.

10-20 EEG Recordings
One monkey (monkey Z) was implanted with an array of electrodes approximating the human 10-20 system locations (FpFz, Fpz, F3, F4, Fz, Cz, C3, C4, Pz, P5, P6, Poz, O1, O2, Oz) (Jasper, 1958). Referencing was done through linked ears. The impedance of the individual electrodes was confirmed to be between 2-5 kOhm at 30 Hz, resembling electrodes used for human EEG. EEG was recorded using a Multichannel Acquisition Processor (Plexon) at 1 kHz and filtered between 0.7-170 Hz. The monkey performed a visual search task necessitating directed spatial attention. Data was aligned to array onset and baseline corrected by subtracting the average activity during the 50 ms preceding the array onset from all timepoints. Data was clipped 20 ms prior to saccade to eliminate eye movement artifacts.

Comparing 10-20 EEG Recordings and Forward Models
To understand the configuration of plausible neural sources contributing to the generation of the N2pc, the N2pc measured in the 10-20 recordings was compared to different compositions of neural sources. 31 different configurations were generated with all possible combinations of V4, lusV4, LIP, 7a, and FEF sources. A single source model was simply the EEG modeled from a single source localized to that dipole location. Combination of source locations were computed as the sum of voltages from each source location for each simulated electrode site. In these models, we only considered the 15 electrode sites present in the 10-20 recording dataset rather than the full 10-10 configuration modeled elsewhere.
Comparisons we performed as a Pearson correlation between each of the 15 empirical recording sites and the 15 simulated sites for each of the 31 model configurations. Therefore, each correlation comprised 15 × 30 data points. Whether a model was significantly correlated with the empirical measurement was evaluated with a Bonferroni corrected p value.

Data Availability Statement
The CSD data used to forward model EEG can be found through Data Dryad (https://doi.org/10.5061/dryad.9ghx3ffm4) or by request from the corresponding author (jacob.a.westerberg@vanderbilt.edu) or senior author (jrieradi@fiu.edu). Code specific to the methods documented here can be obtained from the GitHub repo (https://github.com/ beaherrera/CSDtoEEG-repo) or the corresponding author.

Results
We introduce an approach for estimating the mesoscopic cortical columnar current dipoles from laminar in vivo field potential recordings to determine the contribution of distinct areas to a particular macroscopic EEG signal through forward modeling. The approach consists of the following steps: First, the current sources across depth were determined using CSD methods from laminar field potential recordings. Second, current dipole moments were calculated from the estimated CSD values. Third, the estimated current dipoles were used as current sources of forward models, incorporating their location and orientation. We validated this approach by simulating synthetic field potential recordings generated from detailed biophysical models of cortical pyramidal neurons. Spatial EEG patterns derived from the two approaches were compared to the spatial EEG pattern obtained from a compartmentbased approach (Eq. (5)) (Herrera et al., 2020). We will refer to this as the ground-truth EEG. The first approach derived a single dipole from the simulated CSD (Eq. (4)); we will refer to this as the CSD single-dipole. The second approach derived a single dipole from the transmembrane currents summed over all the simulated L3 and L5 pyramidal cells within a 3 mm diameter cylinder; we will refer to this as the summed transmembrane currents (STC) single-dipole. This methodological framework was then applied to elucidate the neuronal generators of the ERP component known as the N2pc.

Intracranial current sources can be estimated accurately from LFPs
To evaluate the validity of intracranial current sources derived from CSD in the identification of cortical areas contributing to a spatial EEG pattern, we applied the new forward modeling approach to synthetic local field potentials generated from biophysical models of pyramidal cell activity. We first simulated the activity evoked by a noisy current pulse with a mean amplitude of 1.90 nA and a standard deviation of 0.3 nA in 2,200 L3 and 1,000 L5 pyramidal cells randomly distributed within a 3 mm diameter cylindrical cortical column but unconnected from each other or any other neuron (Fig. 2C-D). The current pulses elicited soma Na + spikes in L3 and L5 pyramidal cells followed by Ca 2+ spikes arising in the apical dendrites of L5 pyramidal cells (Fig. 2D). As expected, the L3 and L5 pyramidal cell populations produced distinctly different laminar LFP and CSD patterns (Fig. 2E).
The current dipole moments based on the STC and on the CSD agree reasonably well for each pyramidal cell population separately (Fig. 2F). The summation of the activity across the L3 and L5 pyramidal cells produced LFP and CSD patterns that resembled observed patterns (Fig. 2G), and the current dipole moments calculated through STC or CSD exhibited excellent agreement when (RDM = 0.69; MAG = 0.94) (Fig. 2H).
To simulate the lateralized N2pc component, we modeled two columns of independent, unconnected L3 and L5 pyramidal cells, one in each hemisphere placed on the lunate gyrus surface of extrastriate visual area V4 (Fig. 2B). To simulate the asymmetry observed in the experimental recordings from V4 associated with the focus of attention directed to one hemifield, the mean amplitude of the input to neurons in the left hemisphere was 1.90 nA (representing the focus of attention in the right hemifield), while that to neurons in the right hemisphere was reduced to 1.85 nA. This decreased the probability of eliciting dendritic Ca2+-spikes and reduced the amplitude of the late CSD sink/sources. With this more complete simulation, we calculated the spatial distribution of voltages at two time points based on the ground-truth, the STC, and the CSD approaches and observed virtually indistinguishable patterns (Fig. 2I).
To account for different depths and orientations of cortical columns in these estimates derived from combined L3-L5 activity, to the groundtruth EEG we compared the EEG derived from the STC single-dipole method and from the CSD single-dipole method placing the cylinders of simulated pyramidal cells at 15 randomly selected locations throughout the monkey brain. Relative to the ground-truth EEG, the EEG derived from the STC singledipole corresponded better (RDM = 0.21 ± 0.01; MAG = 1.27 ± 0.02) than that derived from the CSD single-dipole (RDM = 0.52 ± 0.02; MAG = 0.83 ± 0.01) (Fig. 2J). Given the relationship between electric field strength relative to dipole distance, we examined whether RDM or MAG values varied with cortical column depth relative to the scalp. A linear regression of these measures as a function of depth showed no variation for either the STC (R 2 for RDM = 0.23, for MAG = 0.00) or the CSD (R 2 for RDM = 0.04, for MAG = 0.00) (Fig. 2K). Thus, we demonstrate that the synthetic application of this approach is sufficient to yield expected spatio-temporal EEG patterns. We now proceed to employing these methods through forward modeling of empirical data.

Experimental design
Before forward modeling, three objectives must be accomplished. First, monkeys must be trained to perform the cognitive task to proficiency. The representative ERP component used do demonstrate the novel forward modeling approach was the N2pc. The N2pc manifests as a function of the deployment of selective visual attention (Eimer, 1996;Hillyard, 1994, 1990;Woodman and Luck, 1999). To investigate the N2pc, we trained macaque monkeys to perform a visual search task ( Fig. 1A-B) requiring selective attention's rapid deployment. Previous work demonstrates this task elicits an N2pc in macaque monkeys (Purcell et al., 2013). Two macaque monkeys (designated Ca, He) performed visual search for an oddball color target (red or green), presented within an array of 5 uniform distractors (green or red) (N sessions: monkeys Ca, 21; He, 9). Each animal performed above chance [chance level = 0.166] (behavioral accuracy: monkeys Ca, M = 0.88, SEM = 0.021; He, M = 0.81, SEM = 0.024) indicating they understood the task and were selectively deploying attention to accomplish it. Second, the ERP must be obtained to confirm it is being elicited. We sampled extracranial voltage fluctuations from an electrode placed outside the brain and above area V4. An N2pc was observed ( Fig. 1B) with this electrode. Third, suitable and replicable intracranial recordings must be obtained to measure local field potentials and calculate the current dipole during task performance. We obtained neural samples across the cortical layers with linear electrode arrays. After establishing that each penetration was orthogonal to the cortical surface and restricted to a cortical column, the synaptic activation was measured during task performance. The spatiotemporal profiles of synaptic currents are displayed in Fig. 1C for both the attended condition (when the search target was present in the column receptive field) and the unattended condition (when the target was positioned outside the column receptive field). Recordings were restricted to area V4 on the prelunate gyrus, a hypothesized contributor to the N2pc (Westerberg et al., 2022;Woodman et al., 2007) and location where laminar activity orthogonal to the cortical surface can be reliably measured (Nandy et al., 2017). Coincident with the N2pc, we observed differences in synaptic currents across the layers of V4 (Fig. 1D).

Modeling application to in vivo cortical activity
We employed forward modeling to compute the voltage distribution on the cranial surface caused by the translaminar currents in V4 (Nunez et al., 2019;Riera et al., 2012). Acknowledging that the currents generated by one column within area V4 are probably too weak to produce observable cranial voltages, we treated the current measured in one column as representative for the broader cortical tissue within the 3 mm cylinder contributing to the EEG. Based on the average V4 current density (Fig. 3A), we derived the dipole moment ( Fig. 3B) embedded in a boundary element model (BEM) of a macaque monkey head (Fig.  3C). That dipole was strongest when the target was within the receptive field of the column. Lead fields were generated for a dipole on the convexity between the lunate and superior temporal sulci, where neurophysiological samples were taken (Fig. 3C). Forward model solutions were calculated during the N2pc (150-190 ms following array presentation). The resulting spatial distribution of voltages was maximal posterolaterally with higher values contralateral to the target (Fig. 3D). The spatial distribution of the difference in voltages when the target versus distractor was in the receptive field exhibited negativity contralateral to the target that was maximal posterolaterally. These results of the forward model were effectively indistinguishable from previous reports of the N2pc. Moreover, the forward model voltage distributions were robust to data thinning. That is, they were qualitatively similar when using the CSD from a single session (Fig. 4A) as well as the CSD from each monkey individually ( Fig. 4B-C). Thus, currents measured from columns in V4 are sufficient to produce a voltage distribution over the scalp that emulates the observed N2pc.
To investigate the specificity of this relationship, we calculated forward models produced by dipoles placed at other locations in V4. First, placing the dipole at two different locations on the convexity of the prelunate gyrus resulted in cranial voltage distributions that were qualitatively similar to one another as well as to the original dipole location (Fig. 5A-B). This result confirms the robustness of the relationship between current dipoles in V4 and the N2pc. Second, the dipole was placed within V4 in the anterior bank of the lunate sulcus (Gattass et al., 1988) (Fig. 5C). The dipole in this location is not oriented perpendicular to the cranial surface. Although LFP samples have not been measured in this region, there is no a priori reason to assume that the laminar profile of sulcal CSD would be different from that on the gyrus within the same cortical region. The dipole in the sulcus, hereafter referred to as lusV4, resulted in a more posterolateral spatial voltage distribution and similar lateralization relative to target hemifield. However, the forward model voltages were weaker than those derived from dipoles on the gyrus. These findings show that the summation of currents from multiple parts of V4, both gyral and sulcal, contribute to the N2pc. The spatially extensive appearance of the observed ERP thus can be partly understood as arising from variation in the orientation of contributing dipoles.
It seems unlikely that an ERP would arise from just one cortical area. Certainly, other cortical areas contribute to target selection and attention allocation during visual search. Therefore, using this forward modeling tool, we explored the possible contributions of other cortical areas to the N2pc. Lacking a strong a priori rationale for variability off columnar currents between cortical areas, we adopted provisionally the most parsimonious assumption that the dipole measured in V4 is representative of that in other cortical areas.
First, we modeled the possible contribution of two areas in parietal cortex (Fig. 5D). The lateral intraparietal area (LIP), located on the lateral bank of the intraparietal sulcus, is known to contribute to target selection during visual search and is a well-recognized node in the attention network (Bisley et al., 2011;Goldberg et al., 2006;Ipata et al., 2006;Tanaka et al., 2015;Thomas and Paré, 2007). Area 7a, occupying the inferior parietal gyrus, also shows robust attention-related modulation Steinmetz, 2001a, 2001b;Rawley and Constantinidis, 2010;Steinmetz et al., 1994;Steinmetz and Constantinidis, 1995). Unlike LIP, area 7a is oriented to be more conducive to the biophysical generation of potentials measured at the scalp. The forward model of a dipole in LIP produced a voltage pattern with centromedial distribution and weak lateralization relative to target location. This indicates that neural processes in LIP contribute weakly if at all to the N2pc in macaques. Conversely, forward models of a dipole positioned in 7a produced robust potentials measurable at the scalp. However, these potentials were concentrated more centromedially than the posterolateral N2pc.
Next, we modeled the contribution of the frontal eye field (FEF) in the arcuate sulcus of the frontal lobe (Fig. 5E). FEF plays a significant role in visual search and visual attention more generally (Schall, 2015), in part through direct influences on V4 (Moore and Armstrong, 2003;Moore and Zirnsak, 2017). Simultaneous recordings of spikes and LFP in FEF with occipital EEG have demonstrated correlations between LFP in FEF and N2pc polarization (Cohen et al., 2009;Purcell et al., 2013). Nevertheless, the forward model of the dipole in FEF produced a cranial voltage distribution different from the observed N2pc but sharing a modest degree of lateralization relative to target hemifield. The results of forward modeling dipoles in parietal and frontal areas demonstrate that cortical areas contributing little biophysically to an ERP can contribute neuro-computationally to the process indexed by that ERP.
A given cortical area is unlikely to be the only generator of an EEG signal. Thus, the locations and various dipole orientations of diverse sources could explain the spatial distribution of the EEG (Fig. 6). Therefore, we modeled and quantitatively compared all combinations of plausible dipole sources to the N2pc measured from a monkey (identified as Z) performing the same visual search task across 18 recording sessions with an array of 15 EEG electrodes positioned across the cranium (Fig. 7A). Crucially, this monkey did not have a craniotomy that would affect the spatial distribution of the EEG. Computing the forward solution for each of the 5 dipole locations and unique orientations, we linearly summed the solutions for each of the 31 combinations of dipoles across the areas. During the time period of the N2pc in these data (150-190 ms after the visual search array based on data in monkeys Ca and He) we observe a significant N2pc at the posterior EEG recording sites P5 and P6 (t(35) = 2.42, p = 0.02) (Fig. 7B).
For each of the 31 combinations of source locations, we measured the correlation between the empirical data and the modeled data at each of the matching 15 electrode positions. Note, a distinct forward model result was generated for each of the 30 laminar recording sessions meaning each correlation was between the empirical EEG voltage distribution and the 30 sessions of data for each of the 31 model configurations. A single example is shown in the inset of Fig. 7C. Results of all correlations are summarized in Fig. 7C with their corresponding configurations displayed below each bar. While many source configurations resulted in significant correlations with the empirical data, two stand out as superior to the others as indicated by their largest R-squared value -V4 + lusV4 and V4 + lusV4 + LIP. Note both of these models incorporate sources across both gyral and sulcal V4. Also, the better of the two also incorporates a parietal source. We found that the models have highest correlation with empirical values included gyral and sulcal V4 sources (Fig. 7D). The most representative model exhibits the lateralized and posteriorly-distributed signature of the N2pc (Fig. 7E).

Discussion
ERPs have been a powerful tool for investigating brain mechanisms of human cognition for decades. While further insights can be gained from EEG alone, its utility as a biomarker of disorder and index of mechanism is limited by uncertainty about its generation. In other words, knowing the composition of the neural circuitry generating specific cognitive signals in the EEG will improve understanding of what the cognitive signal represents and constrain models of the cortical circuitry enabling the cognitive process indexed by the EEG signal. A systematic search for EEG generators with invasive methods in humans is not possible. Fortunately, as a common animal model for human cognition, macaque monkeys exhibit EEG signals like the N2pc. Even so, key blocks in the bridge from neural activity within cortical columns to the macroscopic EEG signal are missing. Here we provided an approach by which measurement of neural activity across the cortical layers in macaque monkeys during cognitive tasks is forward modeled to infer its impact on the global EEG signal. To our knowledge, this is the first documentation of forward modeling synaptic currents across cortical layers to the EEG signal in awake, behaving primates. It is also the first experimental confirmation in primate cortical tissue of theoretical models developed with data from rodents.
In investigating the neural sources of the ERP component known as the N2pc as a demonstration for this approach, our results provide strong support for the hypothesis that it is generated from neural activity present in occipital and parietal areas. Early reports of the N2pc hypothesized that areas such as V4 contributed to its generation (Luck and Hillyard, 1994) with findings from an MEG study (Hopf et al., 2000) lending more direct support for occipital and parietal sources, specifically. Our comparison of various configurations of plausible neural sources with an empirically measured N2pc distribution suggests that the combination of sources in occipital area V4 and parietal area LIP best models the N2pc. This demonstration serves to highlight the importance of considering the biophysical geometry in EEG production. It is important to note that the dipole used for the modeling was measured in V4 and applied for all other areas in addition to V4. Therefore, neural recordings in LIP are necessary to determine whether a pattern of laminar activation sufficient to produce a dipole at the time of the N2pc is present. While that caveat for LIP remains, we can more confidently conclude that the laminar activation in V4 contributes to the N2pc as our measured activity comes specifically from that area.
While significant insight is gained from evaluating the biophysical plausibility of this variety of brain areas, the results highlight a potentially counterintuitive relationship. That is, FEF, LIP, and 7a-areas definitely contributing computationally to attentional control-do not contribute much biophysically to the N2pc indexing that attention. In particular, a source model of a single dipole in FEF is the least correlated with empirical data, and FEF is not present in the top 9 best-correlated models. This highlights an important consideration -locations of neural computation (e.g., target enhancement, distractor suppression, attention deployment) and locations of the biophysical sources of an ERP indexing those computations need not be the same. That is, while an area like FEF has activation paralleling the N2pc (Cohen et al., 2009;Purcell et al., 2013), its physical position is not sufficient to contribute to its manifestation. In contrast, V4, an area showing more modest attentional modulation than FEF, is oriented in such a way that it is conducive to the generation of electric fields which can be measured at the scalp. However, it is entirely possible that the modulation of activity yielding the difference in dipole present in V4 producing the N2pc originates from FEF . Previous work has established that a robust set of connections exist between V4 and FEF (Anderson et al., 2011;Pouget et al., 2009;Schall et al., 1995;Stanton et al., 1995;Ungerleider et al., 2008) and modulation of FEF activity impacts V4 (Moore and Armstrong, 2003;Noudoost et al., 2014). Work remains to better understand the configuration of neural circuitry generating cognitive EEG signals and importantly, the origins of these signals -not just the areas producing the signal measured at the scalp.
While this study represents an important step in mapping microscopic through mesoscopic to macroscopic signals, this research program is by no means complete. Recent work highlights the strength of investigating the neural mechanisms of EEG at an even finer scale. We can better understand the generation of these potentials by understanding the underlying neural processes at the level of individual synapses through simulations (Herrera et al., 2020;Naess et al., 2021). In our modeling, we compute a dipole from the CSD, which is then used to calculate the definite voltage distribution at the scalp. This process sacrifices laminar information about the functional architecture generating EEG. For example, it is difficult to discern whether an observed source/sink in supragranular layers reflects synaptic activation onto the apical dendrites of L5 pyramidal cells, L2/3 pyramidal cells, or some combination (Pesaran et al., 2018). Additionally, CSD methods obscure any variation of neuronal activity in the radial plane of the column (Pettersen et al., 2006), parallel to the cortical layers.
Several decomposition techniques have been proposed to distinguish the contribution of distinct populations of neurons across layers to the LFPs (Di et al., 1990;Einevoll et al., 2007;Głąbska et al., 2014;Korovaichuk et al., 2010;Łęski et al., 2010;Makarov et al., 2010). Some methods are based on ICA and PCA decomposition methods that assume zero correlation or independence between the brain sources and do not account for the biophysical properties of area-specific cortical circuits (Di et al., 1990;Korovaichuk et al., 2010;Łęski et al., 2010;Makarov et al., 2010). Other techniques, such as the laminar population analysis (Einevoll et al., 2007) and dynamic causal modeling (Pinotsis et al., 2017), fit generative models, constructed from knowledge about local circuit processing and architecture to LFP recordings to estimate distinct intracranial sources. Yet, a systematic translation of these different neuronal sources to individual current dipoles and from them to macroscopic EEG signals is still lacking. Riera et al. (Riera et al., 2012) proposed a biophysical modeling strategy to link these scales. Their approach consisted of estimating the characteristic dynamic equations of intracranial currents from phenomenological models of principal neurons, pyramidal cells, and employing them to develop a generalized inverse solver that accounted for distinct neuronal population dynamics. The estimated brain sources would then be used to reconstruct the dynamics of principal neuronal populations through their phenomenological models. Alternatively, recent studies have inferred the circuit-level dynamics underlying cognitive ERPs by fitting the exogenous drives of a predefined canonical cortical microcircuit model to replicate the event-related current dipoles estimated through inverse modeling (Jones et al., 2009(Jones et al., , 2007Kohl et al., 2022;Law et al., 2022;Naess et al., 2021). The canonical cortical microcircuit model was constructed based on anatomical and electrophysiological findings from sensory cortical areas (Jones et al., 2009(Jones et al., , 2007. New data from the agranular frontal cortex indicates that this microcircuit model does not generalize (Godlove et al., 2014;Ninomiya et al., 2015).
In sum, current research is rapidly generating more and more refined models of EEG production. Increasingly, the role of activity in cortical columns is becoming clearer. Although much of this is done computationally, experimental validation is necessary. Experimental elaboration of cortical columnar activity including the relative contributions of different neuron types is an important next step in the computational work. However, the forward modeling we demonstrate here empirically establishes the missing link between theorized mesoscopic activity patterns to EEG production.

Conclusions
Theoretical models have explained how patterns of activity across the layers of cortical columns can produce EEG. Previous testing of these models has focused on lissencephalic rodent brains. With the more complex macaque brain, we determined that the observed pattern of laminar activation in a cortical area yields spatial distributions of EEG indistinguishable from previously reported ERP components. In so doing, we demonstrated that the dipole derived from the mesoscopic CSD within a cortical column is a more than adequate representative of the transmembrane currents of pyramidal cells. This work resolves a missing link between microscopic neural signals and macroscopic EEG through mesoscopic columnar dipoles. Establishing this link reveals a cortical source for an important EEG component and offers tools to explore how other patterns of laminar activity contribute to different EEG signals.

Fig. 1.
Experimental Design and Laminar Current Source Density. A. Monkeys viewed a fixation point and following a variable delay, a six-item visual search array appeared where one item was saliently different that the others (red among green or green among red). Monkeys shifted gaze to the oddball to receive a juice reward. Magenta circle indicates position of gaze. The oddball sometimes appeared in the receptive field (1/6 of trials) resulting in the attended condition. Blue circle indicates the position of the RF. All other trials were considered unattended. B. Timing of events in each trial. A dot appeared at the centered of the screen and once monkeys successfully made fixation, a 750-1250 ms delay ensued. The array appeared and monkeys made a saccade to the target as rapidly as possible (~125 -500 ms). Monkeys maintained fixation of the target for 500 ms to receive a juice reward. C. Proxy extracranial signal from an electrode placed outside the brain, above area V4. N2pc serves as our representative cognitive ERP indicating directed selective attention and can be observed ~150-190 ms following search array onset and is highlighted in orange. D. Laminar current density computed across electrodes positioned along V4 layers averaged across sessions following alignment relative to the layer 4/5 boundary for both the iCSD method (left) used in this study and Nicholson and Freeman method (right) used in previous reports of these data. Profiles shown for the attention condition (top) and unattended condition (bottom) and are plotted relative to the search array onset. E. Difference in current density profiles between attention conditions (attended -unattended) showing a difference in extragranular currents during the N2pc window. Biophysical forward modeling of synthetic data. A. Boundary element model composite with electrode locations and color-coded surfaces obtained from the NIMH Macaque Template version 2.0 derived from 31 macaque brains. B. Illustration of simulated stronger (black) and weaker (gray) dipoles in area V4 between the superior temporal sulcus and lunate sulcus. C-D. Top, Simulated cortical column of 3 mm diameter comprised of 2200 L3 (C) and 1000 L5 pyramidal cells (D). The somas were randomly distributed in the cylinder at depths corresponding to the vertical extent of deep L3 and of L5 in area V4. Bottom, Raster plots and associated poststimulus time histograms of soma Na + action potentials (black) and dendritic Ca 2+ spikes (red) produced by 100 randomly selected L3 pyramidal cells (C) and L5 pyramidal cells (D) in response to brief suprathreshold stimulation with a noisy current pulse (mean amplitude: 1.9 nA). Each neuron received stimulation with uniform probability between 10 and 20 ms after the launch of each simulated testing trial. E. LFP (left) and CSD (right) evoked by the suprathreshold stimulation of the collection of L3 (top) and L5 (bottom) pyramidal cells. Each neuron received stimulation with uniform probability between 10 and 20 ms after the beginning of each simulated testing trial. F. X, Y, and Z components of the current dipole moment calculated from the summed transmembrane currents of all neurons -Eq. (6) -(STC, black) and from the CSD -Eq. (4) (5)), the STC single-dipole (middle, Eq. (6)), and the CSD single-dipole (right, Eq. (4)). The spatial distributions of cranial voltages derived from each approach (bottom) are difficult to distinguish from one another. J. Comparison of the EEG estimated from the STC single-dipole (black) and the CSD single-dipole (red) to the ground-truth compartment based EEG at 15 random locations in the monkey's brain based on the relative difference (RDM) (top) and magnitude (MAG) (bottom) measures. K. RDM (top) and MAG (bottom) measures comparing ground-truth compartment-based EEG to that derived from STC (black) and CSD (red) as a function of the depth of the center of mass of cylinder of 2200 L3 and 1000 L5 simulated pyramidal cells relative to the scalp. Across 15 random cortical column locations in the macaque brain neither estimate derived from either method varied with depth.   EEG forward models for example session CSD and individual monkey CSD. EEG voltage distribution over the scalp of the BEM simulated with a 61-channel electrode EEG array computed using CSD from an example session (A). and from CSD measured and averaged across sessions from monkey Ca (B) and monkey He (C). Both 3D posterolateral views and 2D disk views are shown for each subset of data. 3D head images are oriented such that the hemisphere in which the recordings were conducted are highlighted (left hemisphere, monkey Ca and example session; right hemisphere, monkey He). Note the difference plots are also computed such that the attention target ipsilateral to the recording chamber is subtracted from the contralateral presentations.

Fig. 5.
Forward models from different bilaterally symmetric dipole locations. The dipole was calculated from the current source density measured 150-190 ms following array presentation. Positions are shown on the enlarged cortex boundary element model surface with filled circle. In visual cortex panels, an additional unfilled circle denotes the position of the dipole used in Fig. 3. Anatomical landmarks are indicated for reference in the inset images with the following abbreviations: sts -superior temporal sulcus, lus -lunate sulcus, plg -prelunate gyrus, cs -central sulcus, ips -intraparietal sulcus, ps -principal sulcus, as -arcuate sulcus. Cranial voltage heatmaps are displayed on posterolateral and 2D disk views for target in contralateral hemifield (left), target in ipsilateral hemifield (middle), and their difference (right). A-B. Dipoles placed at two other sites on the prelunate gyrus generate cranial voltages that were effectively indistinguishable from that observed with the first dipole location. C. Dipoles placed in part of V4 in the lunate sulcus generate a cranial voltage pattern with similar lateralization but different spatial distribution due to dipole orientation. D. The dipoles were also positioned in the lateral intraparietal area (LIP) on the lateral bank of the intraparietal sulcus, a previously hypothesized generator of the N2pc. E. Another location in parietal cortex (area 7a) was chosen as it also demonstrates robust attentional modulation and is located on a gyrus. F. Lastly, FEF was chosen as it has been related to the N2pc previously and also shows robust attentional modulation. It is important to note that the dipoles used throughout all panels were those measured in prelunate V4 and may or may not be a sufficient representation of the empirically measured currents in that area.  Comparison of N2pc voltage distribution over the scalp between empirical measurement and forward models. A. EEG was recorded in a subset of 10-10 configuration electrodes (n = 15) from a monkey with intact skull. B. Voltage potentials for contra-vs. ipsilateral presentations of an attentional target in the visual search task. Differences are observed primarily at posterior sites about 150-190 ms following array presentation. C. Pearson correlation between EEG signal measured empirically (EEG emp ) and the forward modeled EEG signal (EEG mod ) for the difference between a contra-and ipsilateral target presentation at each of the electrode sites present in the empirical recording (n = 15). Negatively correlated configurations were not considered and are displayed as 0. Data was taken as the average difference between 150 and 190 ms following array presentation for both empirically measured and modeled data. EEG model was generated for each laminar recording session (n = 30). Correlation was performed for each possible combination of V4, lusV4, LIP, 7a, and FEF sources (n = 31). R-squared was computed and the bar plot reflects these values sorted from lowest to highest with significant correlation (following Bonferroni correction) indicated above with an asterisk. The combination of sources yielding each bar's value is indicated below by dots with the legend at left. Inset shows the correlation observed with data for a single bar (V4+lusV4+LIP). D. R-squared values were grouped into 0.1 bins. Histogram plotted at top for the frequency of R-squared values. For each of those bins, the prevalence of each source was evaluated where the size of each colored bar indicates the total number of times that source was present in a model resulting in a R-squared value in the respective bin. For example, only V4, lusV4, and LIP sources were present in R-squared value measurements in the 0.3-0.4 bin. E. EEG distribution from the modeled data using the most representative model (V4+lusV4+LIP).