Biophysically detailed forward modeling of the neural origin of EEG and MEG signals

ContentslistsavailableatScienceDirect

NeuroImage

journalhomepage:www.elsevier.com/locate/neuroimage

Biophysically detailed forward modeling of the neural origin of EEG and MEG signals

Solveig Næss

, Geir Halnes

, Espen Hagen

, Donald J. Hagler Jr.

, Anders M. Dale

, Gaute T. Einevoll

, Torbjørn V. Ness

aDepartment of Informatics, University of Oslo, Oslo 0316, Norway

bFaculty of Science and Technology, Norwegian University of Life Sciences, 1432 Ås, Norway

cDepartment of Physics, University of Oslo, Oslo 0316, Norway

dDepartment of Radiology, University of California, La Jolla, CA 92093, USA

eDepartment of Neurosciences, University of California, La Jolla, CA 92093, USA

a b s t r a c t

Electroencephalography(EEG)andmagnetoencephalography(MEG)areamongthemostimportanttechniquesfornon-invasivelystudyingcognitionanddisease inthehumanbrain.Thesesignalsareknowntooriginatefromcorticalneuralactivity,typicallydescribedintermsofcurrentdipoles.Whilethelinkbetween corticalcurrentdipolesandEEG/MEGsignalsisrelativelywellunderstood,surprisinglylittleisknownaboutthelinkbetweendifferentkindsofneuralactivity andthecurrentdipolesthemselves.Detailedbiophysicalmodelinghasplayedanimportantroleinexploringtheneuraloriginofintracranialelectricsignals,like extracellularspikesandlocalfieldpotentials.However,thisapproachhasnotyetbeentakenfulladvantageofinthecontextofexploringtheneuraloriginofthe corticalcurrentdipolesthatarecausingEEG/MEGsignals.

Here,wepresentamethodforreducingarbitrarysimulatedneuralactivitytosinglecurrentdipoles.Wefindthatthemethodisapplicableforcalculatingextracranial signals,butlesssuitedforcalculatingintracranialelectrocorticography(ECoG)signals.Wedemonstratethatthisapproachcanserveasapowerfultoolforinvestigating theneuraloriginofEEG/MEGsignals.Thisisdonethroughexamplestudiesofthesingle-neuronEEGcontribution,theputativeEEGcontributionfromcalcium spikes,andfromcalculatingEEGsignalsfromlarge-scaleneuralnetworksimulations.Wealsodemonstratehowthesimulatedcurrentdipolescanbeuseddirectly incombinationwithdetailedheadmodels,allowingforsimulatedEEGsignalswithanunprecedentedlevelofbiophysicaldetails.

Inconclusion,thispaperpresentsaframeworkforbiophysicallydetailedmodelingofEEGandMEGsignals,whichcanbeusedtobetterourunderstandingof non-inasivelymeasuredneuralactivityinhumans.

1. Introduction

Electroencephalography(EEG) isone of themost importantnon- invasivemethodsforstudyinghumancognitivefunctionanddiagnos- ingbraindiseases(Cohen,2017;Pesaranetal.,2018).Yet,weknow surprisinglylittleabouttheneuraloriginoftheseelectricscalppoten- tials(Cohen,2017):Ontheonehand,wehavearelativelygoodunder- standingofthebiophysicsofEEGs,inknowingthatthesesignalsorigi- natefromcorticalcurrentdipoles,andhavingawell-definedframework forlinkingsuchcorticaldipolestoelectricscalppotentials(Nessetal., 2020;NunezandSrinivasan,2006).Thishasbeentakenadvantageoffor alongtimeinsourcelocalization,byinversemodelingoftheunderlying corticalcurrentdipolesfromEEGdata.Ontheotherhand,eventhough thesecorticaldipolesareassumedtomainlyoriginatefromlargenum- bersofsynapticinputtocorticalpyramidalcellpopulations(Ilmoniemi andSarvas,2019;LopesdaSilva,2013;Nessetal.,2020;Nunezand Srinivasan,2006;Pesaranetal.,2018),thepreciselinkbetweencor- ticaldipolesandtheunderlyingneuralactivityhasremainedunclear.

∗Correspondingauthor.

E-mailaddresses:[email protected](G.T.Einevoll),[email protected](T.V.Ness).

Inotherwords,weknowverylittleaboutexactlywhichtypesofneu- ralactivitythatcauseeventhemostwell-studiedcharacteristicsofthe EEGsignal,suchasdifferenttypesofoscillations(e.g.,alpha,beta,and gammawaves)andstereotypedEEGshapesinresponsetosensorystim- uli(event-relatedpotentials,ERPs)(Cohen,2017).Importantly,these differentEEGcharacteristicsareaffectedinpredictablewaysbyvari- ousbrainconditions,suchassleepandattention(Klimeschetal.,1998;

PalvaandPalva,2011;Siegeletal.,2012),andbybraindisordersin- cludingepilepsy andschizophrenia(Freestoneetal.,2015;Lightand Näätänen,2013; Mäki-Marttunen etal., 2019a;Niedermeyer,2003).

Thismeansthatabetterinsightintohowdifferenttypesofbrainac- tivityisreflectedincorticalcurrentdipolescouldhelpusnotonlyin makingbetterinversemodelsforsourcelocalization,butalsoinprovid- ingabetterunderstandingofthemechanismsofhumancorticalactivity andpossiblycuringbraindiseases(Cohen,2017;Mäki-Marttunenetal., 2019a;Uhlirovaetal.,2016).

ThereasonswhywelackunderstandingoftheneuraloriginofEEG signalsaremany,themainbeing(i)strongethicalconstraintsoninva- sivehumanbrainmeasurementsand(ii)thehighnumberofneurons

https://doi.org/10.1016/j.neuroimage.2020.117467

Received27June2020;Receivedinrevisedform28September2020;Accepted12October2020 Availableonline17October2020

1053-8119/© 2020TheAuthor(s).PublishedbyElsevierInc.ThisisanopenaccessarticleundertheCCBYlicense(http://creativecommons.org/licenses/by/4.0/)

thatcontributetothesignal.However,inrecentyearstherehavebeen majoradvancesinseveralrelevantbranchesofneuroscience,meaning thatabetterunderstandingoftheEEGsignalmaynowbewithinreach (Cohen,2017;Uhlirovaetal.,2016).

Tobypass challenge(i),welooktotherapiddevelopment inthe technologyandmethodsusedtostudyneuralactivityinlabanimals.

Thepossibilitytocontrolandmanipulateneuralactivity,whilesimul- taneouslyrecordingbothintracranialsignalslikethelocalfieldpoten- tial(LFP)(Blomquistetal.,2009;Einevolletal.,2007)andextracranial non-invasivesignalsliketheEEG(Bruyns-Haylettetal.,2017),canbe expectedtomakeimportantcontributionstoourunderstandingofnon- invasivemeasurementsofhumanbrainactivity(Cohen,2017;Lopesda Silva,2013;Pesaranetal.,2018;Uhlirovaetal.,2016).Furthermore,de- tailedbiophysicalmodelingofneuralactivityhasbecomeanimportant toolforunderstandingintracranialLFPmeasurements(Einevolletal., 2013a;Pesaranetal.,2018).GiventhatEEGisexpectedtoreflectthe samebasicprocessasLFP,thatis,largenumbersofsynapticinputto geometricallyalignedpyramidalcells(Buzsákietal.,2012;Nunezand Srinivasan,2006;Pesaranetal.,2018),itseemslikelythatdetailedbio- physicalmodelingcanalsohelpshedlightontheneuraloriginofEEG signals.

Asindicatedin challenge(ii),EEGsignalsareexpectedtoreflect the activityof much larger neuralpopulations than the LFPs,mak- ingthesimulationscomputationallydemanding.Biophysicallydetailed large-scalesimulationsofneuralnetworkshave,however,beengaining substantialmomentumin recentyears,thankstolarge ongoingneu- roscienceinitiatives likeProject MindScopeat theAllenInstitutefor BrainScience,theBlueBrainProjectandtheEUHumanBrainProject (Einevolletal.,2019).ThepossibilitytocalculateEEGsignalsfromsuch existingandfuturelarge-scalebiophysicallydetailedneuralsimulations couldleadtovaluableinsightsintotheneuraloriginoftheEEG.

AnothercomplicatingaspectofEEGmodeling,isthatthesepredic- tionsingeneralrequireaheadmodeltoaccountforthewidelydiffer- entelectricalconductivitiesofthebrain,cerebrospinalfluid(CSF),skull andscalp(IlmoniemiandSarvas,2019;NunezandSrinivasan,2006).

Whilemanysuchheadmodelsexist,theytendtotakecurrentdipoles asinput(NunezandSrinivasan,2006;Pesaranetal.,2018),insteadof thetransmembranecurrentsthatareavailablefrombiophysicalneural simulationsandthatformthebasisformodelingLFPs(Einevolletal., 2013b).

Here,weintroduceanapproachforreducingarbitrarybiophysically detailedsimulatedneuralactivitytocurrentdipoles,whichrepresents anenormousreductionintermofmodelcomplexitywhencomputing brainsignals.Weverifythattheapproachgivesaccurateresultswhen calculatingEEGsignals,butlesssoforintracranialelectrocorticography (ECoG)signals.Next,welookintohowtheapproachcanbe applied forinvestigatingtheoriginofEEGsignals,withaparticularfocuson calciumspikes,beforedemonstratinghowourmethodscanbeapplied forpre-existinglarge-scalenetworkmodels.Finally,weshowhowcur- rentdipolescanbecombinedwithdetailedheadmodels,whichenables simulationofEEGsignalswithunprecedentedbiophysicaldetail.

Notethattheclearseparationbetweencalculationofcurrentdipoles andthecorrespondingEEGisequallyvalidformagnetoencephalogra- phy(MEG)signals.WhileweherefocusmostlyonEEG,thepresented approachforcalculatingcurrentdipolesfromneuralactivityisequally validfor MEGsignals, throughuse ofanappropriateforwardmodel (Hagenetal.,2018;IlmoniemiandSarvas,2019).

2. Methods

Neuralactivitygenerateselectriccurrentsinthebrain,whichinturn createelectromagneticfields.Inthissection,weexplainhowelectric brainsignalscanbemodeledinbothsimpleandmorecomplexvolume conductors.

2.1. Forwardmodelingofelectricpotentials

Weassumenegligiblecapacitiveeffectsinthehead(Micelietal., 2017;PfurtschellerandCooper,1975;Rantaetal.,2017)andthatelec- tricandmagneticsignalseffectivelydecouple.Wecanthenapplythe quasistaticapproximationofMaxwell’sequationsforcalculatingthese signals(Hämäläinenetal.,1993;NunezandSrinivasan,2006).Inother words,forcomputingextracellularelectricpotentials,weenvisionthe headasa3Dvolumeconductor,andcombiningMaxwell’sequations withthecurrentconservationlaw,weobtainthePoissonequationfor computingextracellularpotentials(Griffiths,1999):

𝛁⋅𝐉=𝛁⋅(𝜎𝛁𝜙), (1)

whereJistheelectriccurrentdensityinextracellularspace,𝜎 isthe extracellularconductivityand𝜙istheextracellularelectricpotential.

ThePoissonequationcanbesolvedanalyticallyforsimple,symmetric headmodels,suchasaninfinitelybigspaceandsphericallysymmetric models.Formorecomplexheadmodels,numericalmethodssuchasthe Finite ElementMethod(FEM)can beused (Haufeetal., 2015;Logg etal.,2012;Seoetal.,2016;Vorwerketal.,2014;2019).

2.1.1. Compartment-basedapproach

Extracellularpotentialsgeneratedbytransmembranecurrentscanbe calculatedwithawell-foundedbiophysicaltwo-stepforward-modeling scheme. The first step involves multicompartmental modeling and in- corporates the detailsof reconstructedneuron morphologies forcal- culating transmembrane currents (Sterratt et al., 2011). In the sec- ondstep,Eq.(1)issolvedundertheassumptionthattheextracellular mediumisaninfinitelylarge,linear,ohmic,isotropic,homogeneousand frequency-independentvolumeconductor.Thetransmembranecurrents enteringandescapingtheextracellularmediumcanbeseenascurrent sourcesandsinks,andgivetheextracellularpotential𝜙attheelectrode locationr(Nessetal.,2020),

𝜙(𝐫)= 1 4𝜋𝜎

∑𝑁 𝑛=1

𝐼_𝑛

|𝐫−𝐫_𝑛|, (2)

wherer_nisthelocationoftransmembranecurrentI_n,Nisthenumber oftransmembranecurrentsand𝜎istheextracellularconductivity.This schemeisherereferredtoasthecompartment-basedapproach,andwas appliedusingthePythonpackageLFPy2.0runningNEURONunderthe hood(CarnevaleandHines,2006;Hagenetal.,2018).

2.1.2. Currentdipoleapproximation

Analogoustohowelectricchargescancreatechargemultipoles,a combinationofcurrentsinksandsourcescansetupcurrentmultipoles (NunezandSrinivasan,2006).Fromelectrodynamics,weknowthatex- tracellular potentialsfromavolumeofcurrent sinksandsourcescan bepreciselydescribedbyexpressingEq.(2)asamultipoleexpansion (NunezandSrinivasan,2006):

𝜙(𝑅)= 𝐶monopole

𝑅 +𝐶dipole

𝑅² +𝐶quadrupole

𝑅³ +..., (3)

whenthedistanceRfromthecenterofthevolumetothemeasurement pointislargerthanthedistancefromthevolumecentertothemostpe- ripheralsource(Jackson,1998).Inneuraltissue,thecurrentmonopole contributioniszeroduetocurrentconservation,sincethetransmem- branecurrentssumtozeroatalltimes (Koch,1999;Pettersenetal., 2012).Further,thequadrupole,octopoleandhigher ordertermsare negligiblecomparedtothecurrentdipolecontributionwhenRissuf- ficiently large.Inthiscase,theextracellular potentialfromaneuron modelcanbeestimatedwiththesecondtermofthecurrentmultipole expansion;an approximationknownasthe currentdipoleapproxima- tion(NunezandSrinivasan,2006;Pettersenetal.,2014;Pettersenand Einevoll,2008):

𝜙(𝐫)= 𝐶dipole

𝑅² = 1 4𝜋𝜎

|𝐩|cos𝜃

|𝐫−𝐫p|². (4)

Fig. 1. Illustration of relation between transmembrane currents, axial currents, sourcesandsinks.A:Schematicillustration of a cell model. This toy model only has threecellularcompartments,butnotethatbio- physically detailed neuron models typically have ~ 600–1300 compartments. An exci- tatorysynapticinputinitiatesa currentflow acrossthemembraneandintotheneuron.This currentconsistsof anionic flow of positive ions(e.g.,Na⁺),inadditiontocapacitivecur- rents,andisbyconventionanegativetrans- membranecurrent(I^t),alsoreferredtoasacur- rentsink.Thischangesthemembranepoten- tialatthelocationofthesynapticinput,ini- tiatingaxialcurrents(I^a),thatis,currentsin- sidetheneuron.Theverystrongelectromag- neticattraction of oppositeandrepulsion of equalelectricchargeseffectivelypreventsany chargeaccumulation,ensuringcurrentconser- vation.Thisimpliesthatthesameamountofcurrentthatgoesintoacellularcompartment,mustalsoleavethesamecellularcompartment,enforcingasimple relationshipbetweentransmembranecurrentsandaxialcurrents.Currentconservationalsoensuresthatthesumofalltransmembranecurrentsatanygiventime mustsumtozero,whichimpliesthatanegativetransmembranecurrentcausedbyanexcitatorysynapticinput(currentsink),mustbeexactlybalancedbypositive transmembranecurrentselsewhereonthecell(currentsources).B,C,D:Theextracellularpotentialaroundthecellcanbecalculatedeitherfromthetransmembrane currents(B,Eq.(2)),fromthecurrentdipolemomentsstemmingfromalltheindividualaxialcurrents(C,Eq.(4)),orfromthesinglesummedcurrentdipolemoment (D,Eqs.(6)and(4)).Notethatthesingle-dipoleapproximationisonlyexpectedtobevalidfarawayfromtheneuron,seemaintextfordiscussionofthevalidityof this.

Table1

Radiiandelectricalconductivitiesusedinthefour-spheremodel.Thera- diusofeachsphericalshellinthefour-spheremodel,with𝜎denotingthere- spectiveelectricalconductivities.

Radius (cm) 𝜎(S/m)

Brain 8.9 0.276

CSF 9.0 1.65

Skull 9.5 0.01

Scalp 10.0 0.465

Here,pisthecurrentdipolemomentinamediumwithconductivity 𝜎,𝑅=|𝐑|=|𝐫−𝐫p|isthedistancebetweenthecurrentdipolemoment atr_pandtheelectrodelocationr,and𝜃denotestheanglebetweenpand R.Thecurrentdipolemomentpcanbecalculatedfromanaxialcurrent Iinsideaneuronandthedistancevectordtraveledbytheaxialcur- rent:𝐩=𝐼𝐝,analogoustoachargedipolemoment.Thecurrentdipole approximationisapplicableinthefar-fieldlimit,thatiswhenRismuch largerthanthedipolelength𝑑=|𝐝|(NunezandSrinivasan,2006).

Multi-dipoleapproachFromsomemulticompartmentalneuronsimu- lations(Figs.2–4),wecomputedmultiplecurrentdipolemoments,i.e., oneforeachaxialcurrentflowingbetweenneighboringcompartments intheneuron:

𝐩_𝑘=𝐼_𝑘^axial𝐝_𝑘. (5)

Here, 𝐼_𝑘^axial is an axial current traveling along distance vector d_k, resulting in a current dipole moment p_k. By inserting all the current dipole moments from a neuron simulation into the current dipole approximation (Eq. (4)), we get a good estimate of the ex- tracellular potential at any electrode location where the distance between the electrode and the nearest dipole is sufficiently large (Nunez and Srinivasan, 2006). See Fig. 1 for an illustration of the relation between these different approached for calculating extra- cellular potentials. Note that the length of each (multi-)dipole is equaltohalf thelengthof its corresponding neuronalcompartment.

The calculation of multi-dipoles from simulated neural activity was implemented in LFPy 2.0, and can be used through the function

Cell.get_multi_current_dipole_moments

Single-dipoleapproximationFromeachmulticompartmental neuron simulation,wecomputedonesinglecurrentdipolemoment.Thiscan eitherbedonebysummingupthemultiplecurrentdipolemoments, 𝐩(𝑡)=

∑𝑀 𝑘=1

𝐩_𝑘(𝑡)=

∑𝑀 𝑘=1

𝐼_𝑘^axial(𝑡)𝐝_𝑘, (6)

whereMisthenumberofaxialcurrents,orequivalentlyfromaposition- weightedsumofallthetransmembranecurrents(Hagenetal.,2018;

Lindénetal.,2010):

𝐩(𝑡)=

∑𝑁 𝑘=1

𝐼_𝑘^trans(𝑡)𝐫_𝑘, (7)

whereNisthenumberofcompartmentsinthemulticompartmentalneu- ronmodelandr_kisthepositionoftransmembranecurrent𝐼_𝑘^trans(𝑡).For calculatingEEGsignalsalocationforthecurrentdipolemustbecho- sen,andunlessotherwisespecifiedwepositionedthedipolehalfway betweenthepositionofthesomaandthepositionofthesynapticin- put(formultiplesynapticinputs,weusedtheaveragepositionofthe synapticinputs).Note,however,thatthelargedistancefromtheneu- rontotheEEGelectrode(~ 10 mm)impliesthattheEEGsignalis relativelyinsensitivetosmallchangesinthedipolelocationwithincor- tex. The calculation of current-dipole moments from simulatedneu- ral activitywas implemented in LFPy2.0, andcan be usedthrough

Cell.current_dipole_moment

2.2. Headmodels

Electricpotentialswillbeaffectedbythegeometriesandconductivi- tiesofthevariouspartsofthehead(NunezandSrinivasan,2006),which isespeciallyimportantforelectrodelocationsoutsideofthebrain.This canbeincorporatedintoourextracellularpotentialcalculationsbyap- plyingsimplifiedorcomplexheadmodels.

2.2.1. Four-sphereheadmodel

Thefour-sphereheadmodelisasimpleanalyticalmodelconsisting offourconcentricshellsrepresentingbraintissue,cerebrospinalfluid (CSF),skullandscalp,wheretheconductivitycanbesetindividually foreachshell(NunezandSrinivasan,2006;Srinivasanetal.,1998).The modelsolutionisgiveninNæssetal.(2017)andisfoundbysolvingthe

Poissonequationsubjecttoboundaryconditionsensuringcontinuityof currentandelectricpotentialsovertheboundaries,aswellasnocur- rentescapingtheoutershell.Thismodelisbasedonthecurrentdipole approximation.Theparametersusedinthispaper(Table1)weretaken fromHuangetal.(2013)tobeconsistentwiththeparametersusedin theconstructionofthemorecomplexNewYorkheadmodel(seenext section).

2.2.2. NewYorkHeadmodel

TheNewYorkHeadmodelisadetailedheadmodelbasedonhigh- resolution, anatomical MRI-data from 152 adult heads (Huang and Parra,2015).Themodelwasconstructedbytakingadvantageofthe reciprocitytheorem,statingthatthepositionoftheelectrodeandthe dipolarsourcecanbeswitchedwithoutaffectingthemeasuredpotential (RushandDriscoll,1969).Thismeans,thatvirtuallyinjectingcurrentat thelocationsoftheEEGelectrodesandusingthefiniteelementmethod (Loggetal.,2012)tocomputetheresultingpotentialanywhereinthe brain,givesthelinkbetweencurrentdipolesinthebrainandtheresult- ingEEGsignals(Dmochowskietal.,2017;Huangetal.,2016;Malmivuo andPlonsey,1995;Ziegleretal.,2014).Thislinkwascapturedinama- trixknownastheleadfieldL(NunezandSrinivasan,2006):

𝐋= 𝐄

𝐼 (8)

Here,IistheinjectedcurrentattheelectrodelocationsandEisthe resultingelectricfieldinthebrain.Theleadfieldmatrixgivesusthe preciselinkbetweenacurrentdipolemomentpinthebrainandthe resultingEEGsignals𝚽(NunezandSrinivasan,2006):

𝚽=𝐋⋅𝐩. (9)

Weapplied the NewYork Headmodelby downloadingthe lead fieldLfromhttps://parralab.org/nyhead/. Theunitsincorporatedin the lead field matrix was not immediately obvious. However, from Dmochowskietal.(2017)andHuangetal.(2013)itappearsthatan injectedcurrentIof1mAgivesanelectricpotentialEinV/m,meaning thatacurrentdipolemomentpintheunitofmAmgivesEEGsignalsin theunitofV.

2.3. Simulationofneuralactivity

Allneuronsimulationswereperformedusingthepythonpackage LFPy2.0,runningNEURONunderthehood(Hagenetal.,2018).For investigationsofsingle-cellcontributionstoextracellularpotentials,we appliedthreedifferentmorphologicallyreconstructedcellmodels:The humanlayer-2/3pyramidalcellfrom Eyalet al.(2018), thelayer-5 pyramidalcellfromratcortexconstructed byHay etal.(2011)and a rat layer-5 chandelier cell; an interneuron model developed by Markrametal.(2015).

The pyramidal cell models were downloaded from http://www.senselab.med.yale.edu/modeldb/, with accession num- bers 238347 (2013_03_06_cell03_789_H41_03) and 139653 (cell1) respectively, while we found the interneuron at the Neocortical Microcircuit Collaboration Portal (http://www.bbp.epfl.ch/nmc- portal/microcircuit)underlayer-5,ChandelierCell (ChC),continuous Non-accomodating(cNAC),(rp110201_L_idA_-_Scale_x1.000_y0.975_z1.

000_-_Clone_3).

For all simulations with passive ion channels only (Figs. 2–

4), we used the following cell parameters: membrane resistance of 30000Ωcm²,axialresistanceof150Ωcm(MainenandSejnowski,1996) andamembranecapacitanceof1μF/cm²(Gentetetal.,2000;Sterratt etal.,2011).Whenactivemechanismswereincludedinthesimulations (Fig.5),allcellpropertieswereincorporatedasdescribedinthespecific cell’sdocumentation.

NeuralsimulationsshowninFigs.2–5receivedsynapticinputmod- eled as conductance-based, two-exponential synapses (

Exp2Syn

Table2

Population names and sizes in large-scale neural network modelThenumberofneuronsineachpopulation.E=excitatory, I=inhibitory,andTC=thalamocortical.

Name Population size

L2/3E 20,683

L2/3I 5834

L4E 21,915

L4I 5479

L5E 4850

L5I 1065

L6E 14,395

L6I 2948

TC 902

constantwas3ms,synapticreversalpotentialwas0mVandthesynap- ticweightwassetto0.002μS,unlessotherwisespecified.

2.3.1. Large-scalenetworkmodel

Formodelingofnetworkactivity(Figs. 6and7),weusedtheso- calledhybridschemeproposedbyHagenetal.(2016).Here,theneu- ral network activity is first simulated with point neurons in NEST (Linssenetal.,2018)andtheresultingspikingactivityofallneurons savedtofile.Afterwards,theneuronsaremodeledwithdetailedmulti- compartmentmorphologiesandthespiketimesofthepresynapticneu- ronsareusedasactivationtimesforsynapticinputontotheseneurons inasimulationwheretheextracellularpotentialsarecalculated(Hagen etal.,2016;Senketal.,2018).Thesimulationwasunmodifiedfrom theresultspresentedbyHagenetal.(2016)withtransientthalamocor- ticalinput(theirFigs.1and7),exceptthatallsingle-cellcurrentdipole momentswererecorded,andtheEEGsignalscalculated.Briefly,thenet- workmodelconsistsof8neuralpopulationsacrossfourcorticallayers (L2/3,L4,L5andL6),withoneexcitatoryandoneinhibitorypopula- tionineachofthefourlayers.Thenumberofneuronsineachpopu- lationisgiveninTable2,andtheconnectivitybetween thedifferent populationsisbasedonanatomicaldata(Binzeggeretal.,2004;Potjans andDiesmann,2014),andgiveninHagenetal.(2016)(theirTable5).

Forthefirststep,simulatingthenetworkactivity,thecorticalneurons weremodeledasleakyintegrate-and-fireneurons,connectedwithstatic current-basedexponentialsynapses.Externalinputwassuppliedbothin theformofaconstantcurrentinputwithapopulationspecificstrength, andthalamocorticalinput, whichin thepresentexample correspond tosimultaneousactivationofallthalamocorticalneurons(t=900ms), whichareprojectingtoneuralpopulationsinlayer4andlayer6.Inthe secondstepforcalculatingLFPandEEGsignals,allcellmodelswere passive,withpopulationspecificmorphologies.Theexcitatorypopula- tionswerepyramidalcellsinL2/3,L5andL6,andstellatecellsinL4.All pyramidalcellswereorientedwiththeapicaldendritealongthedepth axisofcortex(z-axis),andrandomlyrotatedaroundthisaxis.Othercell types(stellatecells andinterneurons)wererandomlyrotatedaround allaxes.Toensuresomevariabilityinthemorphologies,the8cortical populationswerefurtherdividedintoatotalof16subpopulationswith differentmorphologies(althoughsomeofthesesubpopulationsusedthe samemorphology).

For a full description of simulation details and parameters used for the large-scale network model, we refer to Hagen et al.

(2016).

2.4. Codeavailability

Simulationcodetoreproduceallfiguresinthispaperisfreelyavail- ablefromhttps://github.com/solveignaess/EEG.git.Foramoregeneral anddetaileddocumentationandexamplesofhowtocalculatecurrent dipolesandEEGsignalsforbiophysicallydetailedcellmodels,werefer thereadertotheLFPydocumentation(https://lfpy.rtfd.io).

Fig. 2. Extracellular potentials become dipolarinthefarfieldlimit.A:Passivelayer- 2/3pyramidalcellfromhuman(Eyaletal., 2016)withanexcitatory,conductance-based, two-exponential synapse placed on apical dendrite (red dot), see Methods (2.3) for parameters. The resulting transmembrane currentsforeachcompartmentareshownas abluearrow(inputcurrent)andredarrows (returncurrents).B: Green arrows represent themultiplecurrentdipolemomentsbetween neighboring neural compartments. C: Gray arrow illustrates the total current dipole moment,thatis,thevectorsumofthedipoles inB.D-F:Extracellularpotentialinimmediate proximityoftheneuron,computedwiththe compartment-based approach, multi-dipole approach and single-dipole approximation, respectively.Notethatthemulti-dipoleresults differ slightly from the compartment-based approach when the distance from the mea- surementpointtothe nearestcurrentdipole momentisshortcomparedtothedipolelength.

G-I:SameasD-F,butatalargerspatialscale (zoomedout).See1mmscalebarinpanelA, DandG.ThecolorbarissharedforpanelsD-I.

3. Results

We introduce an approach for modeling electroencephalography (EEG)andmagnetoencephalography(MEG)signalsfromdetailedbio- physical multicompartment cell models. The approach involves two steps:First,currentdipolemomentsareextractedfromactivityinneu- ronsor networks. Second,the extractedcurrent dipolesareused as sourcesinestablishedforwardmodels.Hereweonlydemonstratethe approachbycomputingEEGsignals,butthecurrentdipolesareequally applicableforcomputingMEGsignalsusingtheappropriatemagnetic- fieldforwardmodels(Hagenetal.,2018;Hämäläinenetal.,1993;Il- moniemiandSarvas,2019).Forillustration,wefirstconsiderEEGsig- nalsstemmingfromsinglesynapticinputontosingleneuronsinanin- finitehomogeneousheadmodel,beforemovingontoasimple,generic headmodel.Finally,westudyEEGsfromlarge-scalesimulatednetwork activity,alsoapplyingadetailedheadmodel.

3.1. Atsufficientlylargedistances,extracellularpotentialsbecomedipolar

Whenmodelingelectricpotentialswithinthebrain,wecanapply thewell-establishedcompartment-basedapproachassumingahomoge- neousvolumeconductor(Section2.1.1)(Einevolletal.,2013a; Holt andKoch,1999).However,thisassumptionisnolongervalidwhenit

comestomodelingEEGsignalsonthescalp,whichcallsforaninhomo- geneousheadmodel(IlmoniemiandSarvas,2019).Suchheadmodels typicallytake currentdipolesasinput,asopposedtoindividualcur- rentsinks/sources, andmustbebasedonthecurrentdipoleapproxi- mation(NunezandSrinivasan,2006).Here,weintroduceanapproach forcomputingcurrentdipolesfromarbitrarysimulatedneuralactivity, andcomparecurrent-basedanddipole-basedmodelingofelectricpo- tentialsgeneratedbyasinglecellreceivingexcitatorysynapticinput.

Excitatorysynapticinputinitiatesanegativecurrentatthesynapselo- cation,sincepositiveionsflowintothecell.Duetocurrentconservation (Koch,1999),thisnegativecurrentisexactlybalancedbyspatiallydis- tributedpositivecurrentsalongthecellularmembrane,asillustratedin Fig.2Aforasingleapicalexcitatorysynapticinputtoapassivehuman corticallayer-2/3pyramidalcellmodel(Eyaletal.,2016).SeeMethods 2.3forsimulationdetails.Inthestandardprocedureformodelingextra- cellularpotentials,herereferredtoasthecompartment-basedapproach, thetransmembranecurrentineachcellularcompartmentcorrespondsto apointcurrentsource/sink.Anotherstrategyistoconsidertheaxialcur- rentofeachcellularcompartmentasasmallcurrentdipole(seeEq.(6)), whichwerefertoasthemulti-dipoleapproach(Fig.2B).Byvectorsum- mationofallthesedipolesintoonesingledipoleataspecificposition, weobtainthesingle-dipoleapproximation(Fig.2C).Forthesakeofcom- paringthesemodelingapproaches,wehaveassumedthatthecellispo-

Fig.3. Single-dipoleapproximationisjus- tifiedforEEGbutnotECoGsignals.A:Illus- trationoffour-sphereheadmodel,wherethe pink,blue,greenandpurplesphericalshells representthebrain,CSF,skullandscalprespec- tively,seeTable1.Thepinkinsetshowsthe humanlayer-2/3neuron(Eyaletal.,2016)lo- catedinthebrain,88.0mmaboveheadcen- ter.41simulationslasting100mswithasingle synapticinputafter20mstocellwithpassive ionchannelsonly,wereperformedforvary- inginput locations,seecolored dots.Thez- componentoftheresultingcurrentdipolemo- mentsfortwosynapticinputlocations(large coloreddots)areshownininsetbelowasfunc- tionsoftime.Theresultspresentedinthisfig- urearecomputedatthesimulationtimepoints producingthelargestcurrentdipolemoment for each synaptic input location. B: Magni- tudeofextracellularpotential|𝜙|asfunction ofdistancefromthetopoftheneuron,shown fortwosimulationswithsynapticinputloca- tionsmarkedbylargecoloreddotsinupperin- setofA.Ineachsimulation,weconsiderthe timepointwiththelargestcurrentdipolemo- ment.Dashedlinesshowextracellularpoten- tialscomputedwithmulti-dipole,andfulllines showsingle-dipolecalculations.C:Relativeer- rorREcomparingthesingle-dipolemodelto themulti-dipolemodel,asfunctionofdistance fromtopofneurontomeasurementpoint.D: RelativeerrorREshowinghowsingle-dipole modeldeviatesfrommulti-dipolemodelEEG calculations,asfunctionofdistancefromsomatosynapselocation.E:MagnitudeofEEGsignal,|EEG|,asfunctionofdistancefromsomatosynapticinputlocation.

F:Relativeerror,RE,showinghowEEGcalculationsperformedwiththesingle-dipoleapproximationdeviatesfrommulti-dipoleapproachasafunctionofamplitude oftheEEGsignal,|EEG|.

Fig.4. EEGsignalsandcurrentdipolemoment from three different cell types with various synapticinput.A:Themorphologiesofahuman L2/3pyramidalcell(blue;Eyaletal.(2016)),arat L5pyramidalcell(red;Hayetal.(2011)),andarat L5interneuron(orange;Markrametal.(2015)).

Theremainingpanelsdisplaydataconnected to eachcelltype,seecell-specificcolors.B:Eachdot representsanexcitatorysynapticinputataspe- cifictime(x-axis)ataspecificheightoftheneu- ron(z-axis,correspondingtopanelA)foraspecific celltype(color).Thebiggerdotswithblackbor- dersmarkinhibitorysynapticinput.Thefourinput bulksrepresent1)100apicalexcitatorysynaptic inputs,2)100basalexcitatorysynapticinputs,3) 400homogeneouslyspread-outexcitatorysynap- ticinputsand4)400homogeneouslyspread-out excitatorysynapticinputsand50inhibitorybasal synapticinputs.Thesynapticweightssumto0.01 μSforallsetsofexcitatory/inhibitorysynapses ineachwave(seeSection2.3fordetails).Forthe interneuron,whichdoesn’thavetypical”apical” or”basal” zones,thesynapseswerespreadoutall overthemorphologyforallinputtypes.C:Thex-, y-andz-componentsofthecurrentdipolemoment pforthethreedifferentcelltypes.D:EEGsignals, 𝜙fromthethreecelltypescomputedwiththefour- spheremodel.

Fig.5. Currentdipolemomentexposeden- dritic calcium spikes. A: Layer-5 cortical pyramidal cell model from rat (Hay et al., 2011), receiving either a single excitatory synapticinputtothe somaevokingasingle somaticactionpotential(bluedot,resultsin B1-4), orin additionanexcitatory synaptic inputtothe apicaldendrite, evokingaden- driticcalciumspikeandtwoadditionalsomatic spikes(orangedot,resultsinC1-4).B1,C1: Membranepotentialatthetwo positionsin- dicatedinA.B2,C2:Extracellularpotential 30μmawayfromthesoma(reddiamondin A),assumingforillustrationaninfinitehomo- geneousextracellularmedium.B3,C3:Single- cellcurrentdipolemoment. B4,C4:Sumof 1000instancesofthesingle-cellcurrentdipole moment(fromB3,C3),thathasbeenrandomly shiftedintimewithanormallydistributedshift withastandarddeviationof10ms.D:Con- tourlinesofextracellularpotentialaroundneu- ronatasnapshotintimeduringthesomatic spikeinB1(t=32.2ms;timemarkedbydashed line).E:Contourlinesofextracellularpotential aroundneuronatasnapshotintimeduringthe calciumspikeinC1(t=36.0ms;timemarked bydashedline).Thesynapticweightwas0.07 and0.15μSforthesomaticandapicalinput location,respectively.

sitionedinaninfinitehomogeneouselectricmedium.Veryclosetothe neuron,theextracellularpotentialwillstronglydependontheexactdis- tributionoftransmembranecurrentsacrossthecellularmorphologyand will,therefore,typicallynottakeapurelydipolarshape(Fig.2D,Ever- susF).However,sincethedipolecontributionwilldominatewhenwe arefurtherawayfromthecurrentsources(seeEq.(3)),theextracellular potentialbecomesmoreandmoredipolarwithincreasingdistancefrom thecell(Lindénetal.,2010).Thisimpliesthatforthepurposeofcalcu- latingextracellularpotentialsfarawayfromthecell,thesingle-dipole approximationmightbewelljustified(Fig.2G-I).Notethattherecan besmalldifferencesbetweentheresultsfromthecompartment-based andthemulti-dipoleapproachesforelectrodelocationsintheimmedi- atevicinityofthecurrentsources,duetotheapproximationsinherent inusingthecurrentdipolemodel(Fig.2DversusFig.2E).

3.2. Single-dipoleapproximationisjustifiedforEEG,butnotECoGsignals

Inordertotesttheapplicabilityofthesingle-dipoleapproximation forcalculatingECoGandEEGsignals,weappliedthefour-spherehead model(Hagenetal., 2018;2019;Næssetal.,2017).Sincethefour- sphereheadmodeltakescurrentdipolesasinput,themulti-dipoleap- proachwasusedasbenchmark:anassumptionthatshouldbewelljus- tifiedforthecell-to-electrodedistancesconsidered,seeSection3.1.

Fordifferentlocationsofasingleexcitatorysynapticinputtoahu- mancorticallayer-2/3pyramidalcellmodel(Eyaletal.,2016)(Fig.3A), wecalculatedtheelectricpotentialatpoint-electrodepositions span- ningfrom100μmabovethetopofthecell,tothesurfaceofthehead, usingboththemulti-dipoleapproachandthesingle-dipoleapproxima- tion(Fig.3B).Inthesimulationsshown,weusedconductance-based synapsesandincludedonlypassivemembraneconductances, butwe confirmedthatusingcurrent-basedsynapsesorafullyactivecellmodel gaveverysimilarresults.

Theelectricpotentialdecreasedsteeplywithdistancewhencrossing thedifferentlayers oftheheadmodel,moststronglyacrossthelow- conductingskull(Fig.3B).Forallsynapticinputlocations,weobserved thattheelectricpotentialcalculatedwiththesingle-dipoleapproxima- tionmarkedlydeviatedfromthemulti-dipoleapproachdirectlyabove theneuron,butthedifferencestronglydecreasedwithdistancefromthe neuron(Fig.3B,fullversusdashedlinesfortwoselectedsynapseloca- tions).Wequantifiedthemodeldissimilaritiesbylookingattherelative erroratthetimepointofthemaximumcurrentdipolemoment,andfora chosendistalsynapticinputtherelativeerrorwas38.9%and0.839%at thepositionoftheECoGandEEGelectrodesrespectively(Fig.3C,green line).Foraspecificproximalsynapticinputweobservedarelativeerror of86.1%attheECoGposition,and13.2%attheEEGposition(Fig.3C, purpleline).Insertingasinglestrongsynapticcurrent(synapticweight 0.05μS)intothesomaofthesamelayer-2/3pyramidalcellwithac- tivemechanisms(Eyaletal.,2018),resultinginasomaticspike,gave relativeerrorsof34.2%and0.813%forthecomputedECoGandEEG signals,respectively(resultsnotshown).WefoundthatcalculatingEEG signalswiththesingle-dipoleapproximationgaverelativeerrorspeak- ingforsynapticlocations ~ 60and400μmabovethesoma,witha sharpdropintherelativeerrorforsynapticinputsfurtherawayfromthe somathan ~ 400μm(Fig.3D).Here,synapticinputsslightlydistalto 400μmawayfromsomaresultedinthemajorityofthereturncurrents escapingthecellbelowthesynapticinput(closertothesoma).Thisgave adistinctlydipole-likesource/sinkdistribution,andtherebylowrela- tiveerrors(~ 0.4%).Synapticinputsslightlyproximalto400μmaway fromsomainsteadresultedinalmostbalancedreturncurrentsabove andbelowthesynapticinput. Thisgaveamultipole-likesource/sink distribution,andtherebylargerrelativeerrors(~ 7%).Note,however, thatthesynapticinputlocationsthatresultedinhigherrelativeerrors, alsogaverelativelyweakEEGsignals(Fig.3E).Thisdemonstratesthat therelativeerrorofthesingle-dipoleapproximationisnegativelycor- relatedwiththeamplitudeofthescalppotential(Fig.3F).Thisisasex-

pected,giventhatthestrongestEEGsignalsareexpectedtobecaused bydipole-likesource/sinkdistributions(section2.1.2).Insummary,the single-dipoleapproximationcanresultinsubstantialerrorsattheposi- tionoftheECoGelectrodes,butgivessmallerrorsatthepositionofthe EEGelectrodesforsynapticlocationsleadingtostrongEEGsignals.

3.3. Single-dipoleapproximationsimplifiesestimateofEEGcontribution

Intheprevioussection,weshowedthatthesingle-dipoleapproxi- mationwasapplicableforcalculationofEEGsignals,andinthissection wedemonstratethatthesingle-dipoleapproximationcansubstantially simplifytheanalysisofthebiophysicaloriginofEEGsignals.

Pyramidalcellshaveapreferredorientationalongthedepthaxisof cortex(herethez-axis),andthedirectionofthecurrentdipolemomentp canbeexpectedtoalignwiththisaxissinceradialsymmetrywilltend tomaketheorthogonalcomponents(p_x,p_y)cancelatthepopulation level(Hagenetal.,2018).Incontrast,interneuronsshowmuchlessof apreferredorientation,andarethereforeexpectedtogiveanegligible contributiontotheEEGsignal,exceptindirectlythroughsynapticinputs ontopyramidalcells(Hagenetal.,2016).Weillustratedthisbyapplying thesingle-dipoleapproximationtothreedifferentcelltypes(Fig.4A), eachreceivingalargenumberofsynapticinputswithtargetregionson thecellssetuptovaryovertime(Fig.4B).

For the previously used human layer-2/3 cell (Fig. 4A, blue;

Eyalet al. (2016)) receivinga volley of 100 excitatorysynaptic in- putsthatwererestrictedtotheuppermost200μmofthecell(t=50ms;

Fig.4B,bluedots),weobservedanegativedeviationofp_z(Fig.4C,blue line).For100basalsynapticinputs(t=100ms;Fig.4B,blueline),the polarityofp_zwasinsteadpositive,butofslightlyloweramplitudethan forapicalinput,ascan beexpectedbecausethelargeareaoftheso- maticregionwillcausestrongreturncurrentsintheimmediatevicinity ofthesynapticinputs,andthereforeanoverallweakercurrent-dipole moment.

Auniformdistributionof400synapticinputsacrossthecellmem- branewitharea-weighted probability(t=150ms;Fig.4B,blueline), onlygaverisetosmallripplesin p_z,duetothesubstantialcancella- tionof currentdipolesof oppositepolarity.Itis sometimesassumed thatexcitatoryinputisrelativelyuniformlydistributedontopyramidal cells,whileinhibitoryinputismoredirectedtotheperisomaticregion (Mazzonietal.,2015;Skaaretal.,2020;Teleńczuketal.,2020;2020).

Asexpected,wefoundthatthiscombinationofuniformlydistributed excitatorysynapticinputandperisomaticinhibitoryinputgaveriseto aclearnegativeresponseinp_z (t=200ms;Fig.4B, blueline),which couldbepartoftheexplanationwhyinhibitorysynapticinputinsome caseshasbeenfoundtodominatetheLFP(Hagenetal.,2016;Teleńczuk etal.,2017).

For a rat cortical layer-5 pyramidal cell model (Fig. 4A, red;

Hayetal.(2011)),theresultingcurrentdipolemomentwasverysim- ilarinshapetothepreviouslydescribedcurrentdipolemomentfrom thehumanlayer-2/3cellmodel. Theamplitudewas however some- what larger, whichwas expectedbecausethelongerapicaldendrite will tend togive larger current dipole moments (Fig. 4C, red line).

Lastly,weusedaratcorticallayer-5interneuronmodel(Fig.4A,or- ange; Markram et al. (2015)), but since the dendrites of interneu- rons are not structured into the same distinctive zones as pyrami- dal cells, the synaptic input caused very small net current dipole moments.

WecalculatedtheEEGsignalswiththefour-sphereheadmodel,us- ingboththemulti-dipole(Fig.4D,dottedlines)andthesingle-dipole (Fig.4D,solidlines)approach.Tocomparetheapproaches,wecom- putedtherelativeerrorasafunctionoftime,thatis,theabsolutediffer- encebetweentheresultsfromthetwoapproaches,normalizedbythe maximumEEGmagnitudecomputed withthemulti-dipoleapproach.

Thesingle-dipoleapproachgaveamaximumerrorof2.14%,3.27%and 0.313%forthehumanlayer-2/3cell,theratlayer-5cellandtheratin- terneuron,respectively.Importantly,theEEGsignalisessentiallyfully

describedbythez-componentofthecurrentdipolemomentp_z,thatis,a singletime-dependentvariable.Thisreductioninsignaldescriptionrep- resentsamassivesimplificationintheunderstandingofthebiophysical originoftheEEGsignal,comparedtoconsideringthetransmembrane currentsandpositionofeachcellularcompartment.

3.4. Currentdipolemomentexposedendriticcalciumspikes

SuzukiandLarkum(2017)recentlydemonstratedthatdendriticcal- ciumspikescanberecordedexperimentallyatthecorticalsurface,and thatthesignalamplitudescanbesimilartocontributionsfromsynaptic inputs.Thisdemonstratesthatactiveconductancesmayplayanimpor- tantroleinshapingECoGandEEGsignals.Furthermore,itsuggeststhat informationaboutcalciumdynamicsmightbepresentinsuchsignals, andthatthisinformationcouldpotentiallybetakenadvantageofwhen studyinglearningmechanismsassociatedwithdendriticcalciumspikes (SuzukiandLarkum,2017).

Thepreviouslyintroducedratlayer-5corticalpyramidalcellmodel from Hay et al. (2011) can exhibit dendritic calcium spikes. When thiscellmodelreceivedasingleexcitatorysynapticinputtothesoma (Fig.5A,bluedot),strongenoughtoelicitasomaticactionpotential (Fig.5B1,blue),asmalldepolarizationwasalsovisibleintheapical dendrite(Fig.5B1,orange).Evenso,thisdidnotinitiateanydendritic calcium spike.However,whencombiningthesamesomaticsynaptic inputwith anadditionalexcitatorysynapticinputtotheapicalden- drite,400μmawayfromthesoma(Fig.5A,orangedot),weobserved adendriticcalciumspike.Thiscalciumspikedid,inturn,inducetwo additionalsomaticspikes(Fig.5C1).Forbothsynapticinputstrategies describedabove,theextracellularpotentialsimulated30μmawayfrom thesomatooktheshapeofstereotypicalextracellularactionpotentials:

thatis,asharpnegativepeakfollowedbyabroaderandweakerposi- tivepeak(Fig.5B2,C2).Further,weobservedthattheslowdendritic calciumspikewasnotreflectedintheextracellularpotentialcloseto thesoma(Fig.5C2).Wefoundthatforthecasewithonlyasomatic spikeandnocalciumspike,thesingle-cellcurrentdipolemomentre- sembledtheinverseoftheextracellularpotential(Fig.5B3),whilefor thecasewithbothsomaticanddendriticspiking,apronouncedslow componentwasalsopresentinthesingle-cellcurrentdipolemoment (Fig.5C3).Somaticactionpotentialsaretypicallynotexpectedtocon- tributesignificantlytoEEGsignals(butseeTeleńczuketal.(2015)), becausetheveryshortdurationofspikeswithbothapositiveandaneg- ativephaseimpliesthatextremesynchronyisneededforspikestosum constructively,andspikesthatareonly partiallyoverlappingtend to sumdestructively.Thesamecannotbeexpectedtoholdforthecalcium spikes,whicharenotonlylonger-lastingbutalsopredominatelycause anegativeresponseinthecurrentcurrentdipolemoment.Tomimica neuralnetworkscenariowithmultiplecellsspikingatslightlydifferent times,wecalculatedthesumof1000instancesofthesingle-cellcurrent dipolemomentthatwasjittered(shifted)intime(normallydistributed, standarddeviation=10ms).Wefoundthatthecasewiththedendritic calciumspikenowhada6.6-foldlargermaximumamplitudethanthe casewithonlythesomaticspike(Fig.5B4versusC4,max|p|=30.8 μAμmand204.2μAμmrespectively).Thisdemonstratesthatdendritic calciumspikesaremuchmorecapableofsummingconstructivelyfora populationofcells,andsubstantiatestheroleofdendriticcalciumspikes inaffectingECoG/EEG/MEGrecordings.

The amplitudeof theslow component of thecurrent dipole mo- mentfromthecalciumspikewasabout0.5μAμm(Fig.5C3),andlater (Sec.3.5)wewillpresentresultsfromasimulatedneuralnetworkwhere theaverageevent-relatedcurrentdipolemomentoflayer5pyramidal cells werefoundtobeabout0.1μAμm(Fig.6D,bottomright).This indicatesthatourresultsarecompatiblewiththeclaimbySuzukiand Larkum(2017)thatsignalamplitudesfromcalciumspikescouldbesim- ilarinamplitudetocontributionsfromsynapticinput.

Wecanmakeaveryroughestimateofthenumberofsimultaneous calcium spikesrequiredtocausea measurableEEG response:Acur-

Fig.6.Large-scaleneural simulationscan beusedtoprobebiophysicaloriginofEEG signals. A: Stimulus-evoked spiking activity fromthalamicinput(time𝑡=900ms,denoted bythinverticalline)inthecorticalmicrocircuit modelfromPotjansandDiesmann(2014).Dots indicatespiketimesofindividualneurons,and populationsare representedin different col- ors(I=inhibitory,E=excitatory).B:Multicom- partmentmodelneuronsusedtoproducethe measurablesignals,withcolorscorresponding topanelA,showingoneexamplemorphology perpopulation.Layerboundariesaremarked atdepthsrelativetocorticalsurface,z=0.A laminarrecordingelectrodewith16contacts separatedby100μm(blackdots)ispositioned inthecenterofthepopulation.C:LFPscalcu- latedatdepthscorrespondingtoblackdotsin B.D:ForthetwoL5populations(L5IandL5E), thethreecomponentsofthecurrentdipolemo- mentisshownforallindividualcells(gray),to- getherwiththepopulationaverage(black).E: Illustrationofthefour-sphereheadmodel,with theredcolumncorrespondingthetheoutline ofthepopulationinpanelB.F:TheEEGsignal fromeachpopulationfoundbysummingthe single-cellEEGcontributionof allindividual cellswithineachpopulation(differentcolors, samecolorschemeasinA,B),togetherwiththe totalsummedEEGsignal(black).Thesimpli- fiedEEGsignalwasfoundbyfirstsummingthez-componentofthecurrentdipolemomentsforallpyramidalcells,thatisL2/3E,L5EandL6E,andcalculatingthe EEGfromthesepopulationdipoles(reddashed).

rentdipolemomentof1μAμmgivesanEEGamplitudeontheorderof 10⁻³𝜇V(seeforexampleFig.4CandD,notedifferentscales).Assum- ingthatanEEGcontributionmustexceed ~ 10μVtobe detectable (Hagenetal.,2018;NunezandSrinivasan,2006)impliesaminimum neededcurrentdipolemomentof ~ 10⁴μAμm.Anumberofperfectly synchronouscalciumspikeswouldeachcontributewith ~ 0.5μAμm (Fig.5C3),suggesting thatabout20,000synchronouscalcium spikes wouldbeneededtocauseameasurableEEGresponse.Further,consid- eringthatthesignalamplitudedecreasesbyabout100-foldfromcorti- calsurfacetoscalp(Fig.3B)andassumingasimilardetectionthreshold, indicatesthatafewhundredsimultaneouscalciumspikeswouldbede- tectablebyECoGelectrodes.

Itmightinitiallyseemsurprisingthatthedendriticcalciumspikeis sostronglyreflectedinthesingle-cellcurrentdipolemoment,giventhat thetransmembranecurrentsassociatedwiththesomaticactionpoten- tialaremuchlargerthanthoseassociatedwiththedendriticcalcium spike:themaximumamplitudeofthetransmembranecurrentsofthe somaticcompartmentwas45.1nA,comparedtojust0.30nAforthe compartmentintheapicaldendrite (Fig.5A, blueandorangedots).

However,thecurrentdipolemomentisgivenastheproductbetween theamplitudeofthecurrentandtheseparationbetweenthesourceand sink(𝐩=𝐼𝐝;Eq.(6)).Whilethecurrentsassociatedwiththesomaticac- tionpotentialwillforthemostpartbecontainedwithinthesomaticre- gion,givingverysmallsink/sourceseparations,thecurrentsassociated withthedendriticcalciumspikewillbedistributedoveramuchlarger partofthecellmembrane.Thiseffectcanbeillustratedbycomparing thespatialprofileoftheextracellularpotentialsaroundtheneuronata snapshotintimeduringasomaticspikeorduringacalciumspike(Fig.5 DversusE).

3.5. EEGfromlarge-scaleneuralnetworksimulations

Sofar,wehaveonlyconsideredEEGcontributionsfromsinglecells, butrealEEG signalsareexpectedtoreflect theactivityof hundreds

ofthousandstomillionsofcells(Cohen,2017;NunezandSrinivasan, 2006).Biophysicallydetailedmodelingoflargepopulationsisstillinits infancy(Einevolletal.,2019)andatpresenttypicallyinclude“only” a fewtens of thousands ofbiophysically detailedcells(Billehetal., 2020;Markrametal.,2015).Networksofpointneurons,ontheother hand,areregularlyusedtosimulatehundredsofthousands(Billehetal., 2020)orevenmillionsofcells(Schmidtetal.,2018;Senketal.,2018), butLFP,ECoG,EEGorMEGsignalscannotbecomputeddirectlyfrom pointneurons(Einevolletal.,2013a;Nessetal.,2020).Toinvestigate EEGsignalsgeneratedbyneuronalnetworks,wethereforeusedahy- bridscheme(Hagenetal.,2016;Senketal.,2018;Skaaretal.,2020), wherethenetworkactivityisfirstsimulatedinahighlycomputationally efficientmannerwithpointneuronsinNEST(Linssenetal.,2018)and theresultingspikingactivityofeachneuronsavedtofile.Afterwards, eachcellismodeledwithbiophysicallydetailedmulticompartmentmor- phologiesandthestoredspikesofallthepresynapticneuronsareusedas activationtimesforsynapticinputontotheseneuronsinasimulation wheretheextracellularpotentialsarecalculated(Hagenetal.,2016;

Senketal.,2018).

We used thelarge-scale point-neuroncorticalmicrocircuit model fromPotjansandDiesmann(2014)andHagenetal.(2016).Thismodel has ~ 80,000neuronsdividedinto8differentcorticalpopulations,that is,oneexcitatoryandoneinhibitorypopulationineachofthefourlay- ersL2/3-L6(seeSection2.3.1).Thismodelcanexhibitadiverseset ofspikingdynamicsincludingdifferentoscillationsandasynchroneous irregular networkstates(Brunel,2000; Hagenetal.,2016).Wehere chosetofocusonthescenariowithtransientthalamocorticalinputfor mimickingevent-relatedpotentials(ERPs).However,focusingondif- ferentneuraloscillations(brainwaves)orspontaneousactivitywould haveservedequallywellforourpurposes.Theonlydifferencefromthe originalsimulationbyHagenetal.(2016)wastheaddedcalculationof currentdipolemomentsandEEGsignals.Wesimulatedtransientthala- micsynapticinputtolayers4and6(Fig.6A),andafterthespikeshad beenmappedontothemulticompartmentcellmodels(Fig.6B),wecal-

Fig. 7. EEGsignals from cortical column networkwithsimpleorcomplexheadmod- els. EEGsignals from thepopulation dipole fromthecorticalmicrocircuitmodel(sameas inFig.6),canbeusedbothwiththecomplex NewYorkHeadmodel(A),orthesimplefour- sphereheadmodel(B).Thedipolewaslocated eitherintheparietallobe(C-G)orintheoc- cipitallobe(H-L).C:Thedipolelocationand orientationin theparietallobe isillustrated withblackarrowonthecorticalsurfaceofthe NewYorkHeadmodel.D,E:EEGsignals(𝜙) onscalpsurfaceelectrodescomputedwiththe NewYorkHeadmodel(D)orthefour-sphere headmodel(E),seenfromabove.Thedatais fromthetime pointofthe strongestcurrent dipolemoment|p|.Dipolelocationismarked byorangestar,withcoordinatesin theNew YorkHeadmodel(55, -49,57) mm.F:EEG tracecomputedwiththeNewYorkHeadmodel (gray)orthefour-sphereheadmodel(black) onthescalpsurfaceelectrodewiththeshort- est distanceto thedipole location. Thedis- tanceswere16.76mm(NewYorkHead)and 16.78mm(four-sphere).G:AbsolutevalueofEEGsignalsfrompanelD,Easafunctionofdistancefromdipoletomeasurementelectrode.H-L:SameasforC-G,but withthedipolelocatedintheoccipitallobe.NotethatpanelI,Jarerotatedtoshowthebackofthehead.ThedipolecoordinatesintheNewYorkHeadmodelwere (-24.3,-105.4,-1.2)mm,andthedistancetotheclosestelectrodeinpanelKwas14.64mm(NewYorkHead)and17.51mm(four-sphere).

culatedtheLFP(Fig.6C)similarlytoHagenetal.(2016)(theirFig.1), inadditiontothecurrentdipolemomentsofeachcell.

Forallcellpopulations,wefoundthatthecurrentdipolemoments fromindividualcellscouldshowlargetransientresponsestothalamic input(Fig.6D; graylinesshow currentdipolemomentfromindivid- ualcellsintwoexamplepopulations:L5inhibitory(L5I)andL5exita- tory(L5E)),butforallinhibitorypopulations,aswellasfortheexci- tatorystellatecellsinL4,thethalamicresponsewasnotvisibleinthe averagecurrentdipolemoment(Fig.6D;blacklines,L5I).Thesame wastrueforthecurrentdipolemomentcomponentsperpendicularto thedepthaxisforexcitatorypopulations(Fig.6D;L5E,p_x,p_y, black lines), butnot for thecomponentalong thedepthaxiswhich hada substantialaverageresponsetothethalamicinput(Fig.6D;L5E, p_z, blackline).Theseobservationsimply,aspreviouslynoted,thatonlythe z-componentofthecurrentdipolemomentfromexcitatorypyramidal cellpopulationscanbeexpectedtocontributesignificantlytotheEEG signal.

OurfindingsinviteasimplifiedapproachtocalculatetheEEGsignal:

Theoriginalapproachinvolvescalculatingallthe ~ 80,000single-cell EEGcontributionsandsummingthem,takingintoaccounttheposition oftheindividualcells,similarlytowhatisdonefortheLFPsignal.A much simpleralternative wouldbe tocomputeasingle summedp_z- componentfromallneuronsineachpyramidalcellpopulation(L2/3E, L5EandL6E),andplaceitinthecenterofthegivenpyramidalcellpop- ulation(withapopulation-specificdepth).Wecan thencalculatethe resultingsimplifiedEEGsignalfromthesepopulationdipoles.Thisap- proximationcanbeexpectedtobeaccuratewhenthepopulationradius issmallcomparedtothedistancefromthepopulationcentertotheEEG electrode.Notethatthedistancefromthetopofcortextothetopofthe headistypically ~ 10mm,whilethediameterofthepresentsimulated populationisonly ~ 1mm(Fig.6;populationoutlineinBisdrawnin redinE).

Totestthissimplifiedapproach, wecombinedthecurrent dipole momentswiththefour-sphereheadmodel(Fig.6E).Wecalculatedthe EEGsignalbythesimplifiedapproach,thatis,fromonetime-dependent array,p_z,foreachpyramidalcellpopulation,locatedinthecenterof thegivenpopulation.Wethencomparedthissimplifiedapproachwith theoriginalapproach,thatis,thesumofEEGcontributions fromall

~ 80,000cells attheirrespective positions.The simplifiedapproach gaveamaximumrelativeerrorof1.1%(maximumabsolutedifference normalizedbymaximumvalueofEEGfromoriginalapproach)(Fig.6F blackversusreddashedline).ThisimpliesthattheEEGsignalfromthe simulatedcorticalactivitycanbenearlyfullyrepresentedbyasingle time-dependentvariablefor eachpyramidalcellpopulation.Further, summingthesepopulationdipolesintoonesingledipole,andlocating itinthecenterofthepopulationcolumnat1mmdepth,insteadgavea maximumrelativeerrorof6.6%.Thishighlightsthattheexactposition of differentpyramidalcellpopulationsarerelativelyunimportantfor shapingtheEEGsignal.

WealsocomparedtherelativeamplitudeoftheEEGsignalfromeach population,andfoundthatforthepresentexample,theexcitatorypop- ulationofL2/3wasthedominantsourceoftheEEGsignal(Fig.6F).

Note,however,thatweexpectthisobservationtobesomewhatmodel- dependent,andthatstronggeneralclaimsaboutthecontributionofdif- ferentpyramidalcellpopulationstotheEEGsignalcannotbemadefrom thisexamplestudyalone.

3.6. Dipoleapproximationincomplexheadmodels

Eventhoughthefour-sphereheadmodelisconvenientforgeneric EEGstudies,manyapplicationssuchasaccurateEEGsourceanalysis, mayrequire moredetailedheadmodels(Daleet al.,1999; Vorwerk et al., 2014). Theconstruction of such complex headmodels is de- pendentonexpensiveequipment,thatismagneticresonanceimaging (MRI),tomaptheelectricalconductivityoftheentireheadatresolu- tionsof ~ 0.5–1.0mm³(HuangandParra,2015;Huangetal.,2016).

Afterwards,numericaltechniquessuchastheFinite ElementMethod (FEM)(Loggetal.,2012)canbeusedtocalculatethesignalattheEEG electrodes forarbitraryarrangementsofcurrentdipolesinthebrain, butatahighcomputationalcost.Thenumberofcomputinghoursis, however,reducedbyapplyingthereciprocityprincipleofHelmholtz.

Thereciprocityprinciplestates,inshort,thatswitchingthelocationof acurrentsourceandarecordingelectrodewillnotaffectthemeasured potential(Dmochowskietal.,2017;Huangetal.,2016;Malmivuoand Plonsey,1995;Ziegleretal.,2014).Thisimpliesthatitsufficestouse

FEMtocalculatetheleadfieldinthebrainfromvirtualcurrentdipoles placedateachoftheEEGelectrodes.Fromtheleadfieldmatrix,wecan inferthepotentialattheEEGelectrodes,givenanarbitraryarrangement ofcurrentdipolesinthebrain.Luckily,severalsuchpre-solvedcom- plexheadmodelsarefreelyavailable,andoneexampleistheNewYork Head(NYH)(Fig.7A),whichwehaveappliedhere(Huangetal.(2016); https://parralab.org/nyhead/).

Toillustratetheuseofpre-solvedcomplexhead-models,weinserted thecurrentdipolemomentobtainedfromthecorticalcolumnmodelin Section3.5intotheNewYorkHeadmodel(Fig.7A),attwomanually chosenpositions:oneintheparietallobe(Fig.7CandD),andoneinthe occipitallobe(Fig.7HandI).Inbothcases,thecurrentdipolemoment wasorientedalongthenormalvectorofthebrainsurface,andtheEEG signalwascalculated.Forcomparisonwitha simplifiedheadmodel, weinsertedthesamecurrentdipolemomentintothefour-spherehead model(Fig.7B)atlocationscomparabletothedipolepositionschosen intheoccipitalandparietallobeintheNYHmodel:thelocationsinthe four-spheremodelwerechosenclosetothebrainsurface,suchthatthe distancefromdipolepositiontotheclosestelectrode(Fig.7EandJ)and thebrainsurfacenormalvectorsweresimilartotherespectivepositions intheNYHmodel.

ThetwoheadmodelsgeneratedEEGsignalsofthesametemporal shape,whichisexpected,giventhatneitherofthemodelsincludedany temporalfiltering(Micelietal.,2017;PfurtschellerandCooper,1975;

Rantaetal.,2017).ThecomputedEEGsignalsfromthetwoheadmodels alsogavecomparableresultsinbothspatialshapeandamplitude(Fig.7 DversusE;IversusJ).TherelativedifferencebetweentheEEGsignals calculatedwiththefour-spheremodelandtheNYHmodelatthetime ofmaximumsignalamplitudewas201%and25.2%forthepositionsin theparietalandoccipitalloberespectively(Fig.7FandK).Notethat whilethefour-sphereheadmodelgaveverysimilarEEGamplitudesfor thetwodifferentdipolelocations(asexpectedfromsymmetry),theEEG amplitudesfromthecomplexheadmodelwasmuchmorevariable,even forsimilardistancestotheclosestelectrode(Fig.7FandK).

Thehighervariabilityofthecomplexheadmodelwasalsoapparent inthedecayofthemaximumEEGamplitudewithdistance,whichwas perfectlysmooth,exponential-like(NunezandSrinivasan,2006),and verysimilarforthetwolocationsinthefour-spheremodel,butvery variableforthecomplexheadmodel,althoughwiththesamegeneral shape(Fig.7GandL).

Notethatdespitethecomplexity,theNYHmodelis substantially fasterthanthefour-spheremodel.InordertosimulatetheEEGsfroma dipolemomentvectorcontaining1200timesteps,theNYHmodelexe- cutiontimeswere~ 0.3s,whilethefour-spheremodelneeded~ 0.9s.

4. Discussion 4.1. Summary

Inthispaper,wehaveintroducedanapproachforreducingarbitrary simulatedneuralactivityfrombiophysicallydetailedneuronmodelsto singlecurrentdipoles(Fig.2).Weverifiedthattheapproachwasap- plicableforcalculatingEEG,butgenerallynotforECoGsignals(Figs.3 and4),andgaveexamplesofhowreducing neuralactivitytoasin- gledipolecanbeapowerfultoolforinvestigatingandunderstanding single-cellEEGcontributions(Figs.4and5).Furthermore,wedemon- stratedthatthepresentedapproachcouldeasilybeintegratedwithex- istinglarge-scalesimulationsofneuralactivity.Moreover,weshowed howsingledipolesareusefulforconstructingcompactrepresentations oftheEEGcontributionsfromentireneuralpopulations,withmethods stillfirmlygroundedintheunderlyingbiophysics(Fig.6).Finally,we demonstratedhowthesimulatedcurrentdipoles,fromsinglecellsor large neuralpopulations, can bedirectly inserted intocomplexhead modelsforcalculatingmorerealisticEEGsignals(Fig.7).

4.2. ApplicationofcurrentdipolesforcomputingEEG,MEGandECoG signals

Wehavehighlightedthatthecalculationofcurrentdipolesfromneu- ralactivityiscleanlyseparatedfromthecalculationoftheensuingEEG signals.SinceMEGsensorslikeEEGelectrodesarepositionedfaraway fromtheneuralsources,thesameistrueforMEGsignals.Thecalculated currentdipolescanthereforealsobeusedincombinationwithsimpli- fiedordetailedframeworksforcalculationofMEGsignals,forexample byfollowingmethodsoutlinedinHagenetal.(2018)andIlmoniemiand Sarvas(2019).

ECoG electrodes are in general positioned closer to the neural sources.ForourexamplesimulationsoftheECoGsignalgeneratedby individualneurons,wefoundthatuseofthesingle-dipoleapproxima- tiongavesubstantialerrors(Fig.3).ThusforcomputationofECoGsig- nals,thestandardcompartment-basedformalismorthemulti-dipoleap- proach(Fig.2)requiringmuchmorecomputationalresources,maybe required. Hereanalternativetousingfullheadmodelsis touse the method ofimages,taking intoaccountthediscontinuityof electrical conductivityatthecorticalsurface(Hagenetal.,2018;Pettersenetal., 2006).AfurthercomplicationofECoGsignalsisthatsincetheelectrodes arebothrelativelylarge andclose totheneuralcurrentsources, the ECoGelectrodesthemselvesmightbeexpectedtosubstantiallyimpact themeasuredsignals(Nessetal.,2015;Rogersetal.,2020;Vermaas etal.,2020).

NotethatwhilewehereusedLFPy2.0(Hagenetal.,2018;2019),a pythoninterfacetoNEURON(CarnevaleandHines,2006),calculation ofcurrentdipolemomentscaneasilybeimplementedintoanyframe- workwherethetransmembranecurrentsareavailable,throughthesim- pleformulagivenineq.(7).

4.3. Generalizationtonon-compartmentalmodels

EEGandMEGrecordingsreflectneuralactivityatthesystems-level (Einevoll etal., 2019; Pesaran etal., 2018). Here,we havefocused on calculatingcurrentdipolesfrom detailedmulti-compartmentneu- ronmodels,butneuralmodelingatthesystems-levelisoftenbasedon higherlevelsofabstraction,likepointneurons(Linssenetal.,2018)or firingratepopulations(Sanz-Leonetal.,2013).Calculationofelectric ormagneticsignalsfromsuchhigher-levelneuralsimulationsmustin generalrelyonsomekindofapproximationtrick,sinceneuronsrequire aspatialstructuretobecapableofproducingelectromagneticsignals (Einevolletal.,2013a).Onesuchtrickthatwetookadvantageofhere isthehybridscheme(Hagenetal.,2016).Thistwo-stepschemeinvolves neuralnetworkactivityfirstbeingsimulatedbypointneurons,before theresultingspiketrainsarereplayedontomulti-compartmentneuron modelsforcalculatingLFPandEEGsignals(Section3.5;Fig.6).

Further,thehybridschemecanbegeneralizedtoalsoallowforcal- culationofEEG/MEGsignalsfromfiring-ratemodelsbyusingtheso- calledkernelmethod,whichhaspreviouslybeensuccessfullyappliedto theLFP(Hagenetal.,2016;Skaaretal.,2020;Teleńczuketal.,2020).

Inpractice,thiscanbedoneintwosteps:First,simultaneouslyactivat- ingalloutgoingsynapsesfromaspecific(presynaptic)simulatedpop- ulation,andrecordingthetotalcurrentdipolemomentoftheresponse (thekernel)(Hagenetal.,2016).Second,computingtheEEG/MEGcon- tributionstemmingfromthis (presynaptic)populationbyconvolving thekernelwiththepopulationfiringrate,andapplyinganappropri- ateforwardmodel.Here,thefiringratewouldbeobtainedseparately inpointneuronnetworkmodelsorfiring-ratemodels.Inthisway,the basicbiophysicsofEEG andMEGsignalsfrom synapticactivationof multi-compartmentneuronmodelsisincluded,avoiding,however,com- putationallyheavymulticompartmentalmodelingofspikingdynamics.

Thecalculatedcurrent-dipolekernelsshouldbeapplicablefordifferent kindsofinputtotheoriginalnetworkmodel,butwouldingeneralhave toberecomputedforchangestocellorsynapticparameters.

4.4. Connectiontootherwork

Calculationofcurrentdipolemomentsfrommorphologicallycom- plex cellmodelshasbeen pursued before,for example tostudythe EEG and MEG contribution of spiking single cells (Murakami and Okada,2006),ortostudyhowthesynapticinputlocationaffectsthecur- rentdipole(AhlforsandWrehII,2015;Lindénetal.,2010).Important workonEEGinterpretationintermsoftheunderlyingneuralactivity hasalsopreviouslybeendonethroughuseof”minimallysufficient” bio- physicalmodels,seeforexampleMurakamietal.(2003,2002),Jones etal. (2009,2007),Sliva etal. (2018)andNeymotin etal. (2020). Here,”minimallysufficient” meansthatthecellmodelsonlyhadmini- mallyneededmulti-compartmentspatialstructure(pointneuronscan- notproducecurrentdipolemoments),onlyconsideredafewcelltypes, andemployedsimplesynapticconnectionrules.Inparticular,theHu- manNeocorticalNeurosolver(HNN)(Neymotinetal., 2020)enables researcherstolinkmeasuredEEGorMEGrecordingstoneuralactivity throughapre-definedcanonicalneocorticalcolumntemplatenetwork.

HNNcomeswithaninteractiveGUI,designedforuserswithlittleorno experienceincomputationalmodeling,andmightthereforebeanap- propriatechoiceforresearchersseekingtogainabetterunderstanding oftheirEEG/MEGdata.However,whiletheuseofsuchminimallysuffi- cientmodelsallowsforquickanddirectcomparisonbetweensimulated andrecordedEEGsignals,itis not(presently)compatiblewithsimu- latingEEGorMEGsfrombiophysicallydetailedsinglecell-ornetwork models,constructedfromdetailedexperimentaldata(Arkhipovetal., 2018;Billehetal.,2020;Eggeretal.,2014;Gratiyetal.,2018;Hagen etal.,2016;Markrametal.,2015;Reimannetal.,2013).

Amorehigh-levelapproachforsimulatingMEG/EEGsignalsfrom theunderlyingneuralactivityhasbeenpursuedthroughneuralfield orneuralmassmodels(Bojaketal.,2010;Coombes,2006;Davidand Friston,2003;Decoetal.,2008;Jirsaetal.,2002;Ritteretal.,2013), whichaimtomodeltheevolutionofcoarse-grainedvariablessuchas themeanmembranepotentialorthefiringrateofneuronpopulations.

Suchcoarse-grainingdrasticallyreducesthenumberofparametersand thecomputationalburdenofthesimulation,andcanbeusedtostudythe interplayamongentirebrainregions,andindeedrunwhole-brainsimu- lations.TheVirtualBrain(TVB)isanexcellentexampleofasoftwarefor whole-brainnetworksimulations(Ritteretal.,2013;Sanz-Leonetal., 2015;2013),wheredetailedandpotentiallypersonalizedheadmod- elscanbecombinedwithtractography-basedmethodsidentifyingthe connectivitybetweenbrainregions(Sanz-Leonetal.,2013).Tocalcu- latemeasurementmodalitieslikeMEGand/orEEGsignalsfromneural fieldorneuralmassmodels,itistypicallyassumedthatthepopulation currentdipolemomentsareroughlyproportionalto,forexample,the averageexcitatorymembranepotential(Bojaketal.,2010;Ritteretal., 2013).Further,EEGscanbecalculatedfromtheresultingcurrentdipole momentsincombinationwithheadmodelsaspresentedinthispaper, orthroughothersoftwaresortechniques(Gramfortetal.,2014).This suggestsanintriguingfuturedevelopment,whereonecouldapplythe above-mentionedkernelmethodbasedonbiophysicallydetailedneu- ronmodelstosubstantiallyincreasetheaccuracyofLFP,EEGandMEG predictionsfromhigh-levellarge-scalesimulationsofneuralactivity.

4.5. Outlook

EEGand,later,MEGsignalshavebeenanimportantpartof neu- roscienceforalongtime,butstillverylittleisknownabouttheneu- raloriginofthesignals(Cohen,2017).Abetterunderstandingofthese signalscouldleadtoimportantdiscoveriesabouthowthebrainworks (Ilmoniemi andSarvas, 2019; Lopes da Silva, 2013; Pesaran et al., 2018;Uhlirovaetal.,2016),andprovidenewinsightsintomentaldis- orders(Mäki-Marttunenetal.,2019a; Sahinetal., 2019).Thiswork layssome of thefoundationfor obtaining abetter understandingof EEG/MEGrecordings,byallowingeasycalculationofthesignalsfrom arbitraryneural activity. Thepresented formalism is well suited for

modeling EEG/MEGcontributions from variouspotential neural ori- gins,includingdifferent celltypes,differention channelsanddiffer- ent synaptic pathways. Forexample, to studythe effect of calcium spikes(SuzukiandLarkum,2017),I_hcurrents(Kalmbachetal.,2018;

Ness et al., 2016; 2018), or geneexpression on EEG signals (Mäki- Marttunenetal.,2019b),oneonlyneedstoknowhowthez-component oftheresultingpopulationcurrentdipoleisaffected.Thisdecouplingof thecurrentdipolemomentandheadmodelallowsforeasierinvestiga- tionandimprovedunderstandingoftheoriginoftheEEG/MEGsignal.

This studyhasonly consideredhow onecan calculateEEG/MEG signalsfromtheunderlyingneuralactivity(so-calledforwardmodel- ing).EEG/MEGmeasurementsare,however,oftenusedforsourcelocal- ization(so-calledinversemodeling),aimingtoidentifytheunderlying corticalcurrentdipoles(Gramfortetal.,2014;IlmoniemiandSarvas, 2019;NunezandSrinivasan,2006).However,suchreconstructedcur- rentdipolesaregenericinthesensethattheyaretypicallynotintended torepresentspecificneuralpopulations.Byallowingforcalculationof currentdipolesfromcorticalpopulations,theworkpresentedheretakes asteptowardsconsolidatingthe,sofar,mostlyseparatescientificdis- ciplinesofneuralmodelingandEEG/MEGdataanalysis(butseealso Neymotinetal.(2020)).

The main challenge in inverse modeling, is that the problem is mathematically ill-posed, because thenumber of possible current sinks/sourcesismuchlargerthanthenumberofrecordingelectrodes.

Thisimpliesthatuniquesolutionscaningeneralnotbefound,without makingadditionalassumptionsoncorticalcurrentdipoledistributions (Pettersenetal., 2012).Inpractice,this meansthattheperformance ofinversemethodswilldependondesignchoicesthatneedtesting.The presentedapproachcanbeusedtocreatesimulatedEEG/MEGdatawith knownunderlyingcurrentdipoles,whichcanbeusedforbenchmark- inginversemethodsforsourcelocalization,similartowhathasprevi- ouslybeendoneforcurrentsourcedensityestimations(Pettersenetal., 2006).Further,Pettersenetal.(2006)demonstratedthatanimproved inversemethodforLFPanalysis,theso-callediCSDmethod,couldbe madebybuildingtheforwardmodelintotheinversemodel.Goingone stepfurther,wecanenvisionincorporatingthepresentedapproachinto anEEGinversemodel,aimingtoidentifysynapticinputregionstodif- ferentpyramidalcellpopulationsinsteadofcorticalcurrentdipoles.

Whiletherearemanyexamplesofdetailedbiophysicalmodelingof neuralactivityimprovinginterpretationofmeasuredintracranialextra- cellularpotentialsinlabanimals(Blomquistetal.,2009;Chatzikalym- niouandSkinner,2018;Einevolletal.,2007;Luoetal.,2018;McColgan etal.,2017;Teleńczuketal.,2020),muchlesshasbeendoneforhu- manEEG/MEGsignals.Thisisnaturalgiventhatstudiesofhealthyhu- manbrainsnecessarilyarelimitedtonon-invasivetechnologies(Cohen, 2017;LopesdaSilva,2013;Uhlirovaetal.,2016).However,givenall the valuableinsightsthatcould be gained fromanincreased under- standingofnon-invasivemeasurementsofneuralactivityinhumans,an importantchallengeinmodernneuroscienceistobuildonthemecha- nisticinsightsfromanimalstudiesandusethemforunderstandingnon- invasivesignalsinhumans(Cohen,2017;Einevolletal.,2019;Lopesda Silva,2013;Mäki-Marttunenetal.,2019a;Uhlirovaetal.,2016).The approachforcalculatingEEG/MEGsignalsinthispapershouldthere- foreideallybeusedincombinationwithanimalstudiessimultaneously measuringmultisitelaminarLFP(andMUA)signalswithincortex,as wellasEEG/MEGsignals(seeforexampleBruyns-Haylettetal.,2017;

Cohen,2017).

Today,wehaveareasonablygoodunderstandingofhowsingleneu- rons operate, thatis, howthey respond tosynapticinput, and how multitudes of synaptic inputs combine to produce action potentials (Einevoll etal., 2019). Similarly, we can,toa highdegree, explain themeasurementphysicsofEEG/MEG,thatis,howneuralcurrentsaf- fectelectromagneticbrainsignalsrecordedoutsideofthehead(Cohen, 2017;IlmoniemiandSarvas,2019;NunezandSrinivasan,2006).The challengeofunderstandingEEG/MEGsignalsisthereforecloselyrelated tothegreatestchallengeinmodernneuroscience:understandingneural