A stepwise neuron model fitting procedure designed for recordings with high spatial resolution: Application to layer 5 pyramidal cells

ContentslistsavailableatScienceDirect

Journal of Neuroscience Methods

j ou rn a l h o m epa g e :w w w . e l s e v i e r . c o m / l o c a t e / j n e u m e t h

A stepwise neuron model fitting procedure designed for recordings with high spatial resolution: Application to layer 5 pyramidal cells

Tuomo Mäki-Marttunen

, Geir Halnes

, Anna Devor

, Christoph Metzner

, Anders M. Dale

, Ole A. Andreassen

, Gaute T. Einevoll

aNORMENT,KGJebsenCentreforPsychosisResearch,InstituteofClinicalMedicine,UniversityofOslo,Oslo,Norway

bSimulaResearchLaboratory,Lysaker,Norway

cFacultyofScienceandTechnology,NorwegianUniversityofLifeSciences, ˚As,Norway

dDepartmentofNeurosciences,UniversityofCaliforniaSanDiego,LaJolla,CA,USA

eDepartmentofRadiology,UniversityofCalifornia,SanDiego,LaJolla,CA,USA

fMartinosCenterforBiomedicalImaging,MassachusettsGeneralHospital,HarvardMedicalSchool,Charlestown,MA,USA

gBiocomputationResearchGroup,UniversityofHertfordshire,Hatfield,UK

hDivisionofMentalHealthandAddiction,OsloUniversityHospital,Oslo,Norway

iDepartmentofPhysics,UniversityofOslo,Oslo,Norway

h i g h l i g h t s

•NewVSDandCa²⁺-imagingtechniquesallowhigh-resolutionimagingofneurons.

•Wepresentastepwisemodel-fittingschemewithpossibilitiestoapplytosuchdata.

•Weapplyourmethodtosimulateddatatoconstructareduced-morphologyL5PCmodel.

•Ourmodeliscost-efficientandreproducesthemainfeaturesoftheoriginalmodel.

•OurmodelpredictsthatinterconnectedL5PCscanamplifylow-frequencyinputs.

a r t i c l e i n f o

Articlehistory:

Received8September2016

Receivedinrevisedform7September2017 Accepted5October2017

Availableonline7October2017 Keywords:

Multi-compartmentalneuronmodels Biophysicallydetailedmodeling Modelfittingusingimagingdata Automatedfittingmethods Parameterpeeling

a b s t r a c t

Background:Recentprogressinelectrophysiologicalandopticalmethodsforneuronalrecordingspro- videsvastamountsofhigh-resolutiondata.Inparallel,thedevelopmentofcomputertechnologyhas allowedsimulationofever-largerneuronalcircuits.Achallengeintakingadvantageofthesedevelop- mentsistheconstructionofsingle-cellandnetworkmodelsinawaythatfaithfullyreproducesneuronal biophysicswithsubcellularlevelofdetailswhilekeepingthesimulationcostsatanacceptablelevel.

Newmethod:Inthiswork,wedevelopandapplyanautomated,stepwisemethodforfittinganeuron modeltodatawithfinespatialresolution,suchasthatachievablewithvoltagesensitivedyes(VSDs)and Ca²⁺imaging.

Result:Weapplyourmethodtosimulateddatafromlayer5pyramidalcells(L5PCs)andconstructamodel withreducedneuronalmorphology.Weconnectthereduced-morphologyneuronsintoanetworkand validateagainstsimulateddatafromahigh-resolutionL5PCnetworkmodel.

Comparisonwith existingmethods:Ourapproachcombinesfeaturesfromseveralpreviouslyapplied model-fittingstrategies.Thereduced-morphologyneuronmodelobtainedusingourapproachreliably reproducesthemembrane-potentialdynamicsacrossthedendritesaspredictedbythefull-morphology model.

Conclusions:Thenetworkmodelsproducedusingourmethodarecost-efficientandpredictthatinter- connectedL5PCsareabletoamplifydelta-rangeoscillatoryinputsacrossalargerangeofnetworksizes andtopologies,largelyduetothemediumafterhyperpolarizationmediatedbytheCa²⁺-activatedSK current.

©2017TheAuthors.PublishedbyElsevierB.V.ThisisanopenaccessarticleundertheCCBYlicense (http://creativecommons.org/licenses/by/4.0/).

∗Correspondingauthorat:SimulaResearchLaboratory,Lysaker,Norway.

E-mailaddress:tuomo@simula.no(T.Mäki-Marttunen).

1. Introduction

Automatedmethods for neuron model fitting have replaced theneedformanualtuningofmodelparameters(VanGeitetal., https://doi.org/10.1016/j.jneumeth.2017.10.007

0165-0270/©2017TheAuthors.PublishedbyElsevierB.V.ThisisanopenaccessarticleundertheCCBYlicense(http://creativecommons.org/licenses/by/4.0/).

2008).Duetotheease oftheiruse, theycouldprovidea solu- tion for exploiting computational properties of single neurons andneuralcircuits(Markrametal.,2015).Novelalgorithmsand strategiesforautomatedneuronmodelfittinghavebeenproposed (Kerenet al.,2005;Druckmann etal.,2007,2011; Hendrickson et al., 2011; Bahl et al., 2012; Brookings et al., 2014; Rumbell etal.,2016).Thesemethodsspanawiderangeoftypesofneu- ronsandtheirelectrophysiologicalcharacteristics.Manyofthese strategiesonlyusevoltagetracesrecordedfromthesoma,while others rely on electrophysiological recordings at one or more additionaldendritic locationsin orderto reproducethecorrect membrane-potentialdynamicsdistributedacrossthesub-cellular compartments. However, recording with multiple intracellular electrodesisanexperimentallydemandingprocedureultimately limitedbythenumberofmicromanipulatorsthatcanfitinasetup, theexperimenter’sskills,andintegrityofthecellinthepresence ofmultiplerecordingelectrodes(Wangetal.,2015).Bycontrast, recent developments of opticalimaging technologiesand engi- neeringofnovelvoltage-sensitivedyes(VSDs)andCa²⁺indicators haveenabledhigh-resolutionsamplingoftransmembranevoltage andintracellularCa²⁺concentrationinsingleneuronswithsub- cellularresolution(belandHelmchen,2007,1;Peterkaetal.,2011;

Hochbaumetal.,2014;Grienbergeretal.,2015;Anticetal.,2016;

Louetal.,2016).Inthiswork,wedevelopanautomated,stepwise procedureforfittingamulticompartmentalneuronmodeltodata fromsomaticpatch-clamprecordingsincombinationwithVSDand Ca²⁺-imagingdata.

Interactionsbetweensynapticinputstothedendritesandfir- ing of the somaare a hallmark of neural computation (Smith etal.,2013).ThisisespeciallytrueforL5PCs,whicharecharac- terizedbyalongapicaldendritethatspansacrosscorticallayers andreceivesinputsfromvariousneuronpopulationsindifferent partsofthedendritictree(Larkum,2013).Theapicaldendriteis rich involtage-gated Ca²⁺channels thatcontribute tothegen- erationofa dendriticCa²⁺spike(Schilleretal.,1997).ThisCa²⁺

spikeplaysanimportantroleinintegrationofsynapticinputsto theapicaltuft,communicationofthesesignalstothesoma,and coincidencedetectionintheformoftheback-propagatingaction potential-activatedCa²⁺spike(BAC)firing(Larkum,2013).L5PCs expressmanytypesofvoltage-gatedionchannels(Korngreenand Sakmann,2000;Christopheetal.,2005;Almogetal.,2009),and anumberofcomputationalmodelshavebeendevelopedaccount- ingforthesebiophysicalproperties(Kerenetal.,2005,2009;Bahl etal.,2012;MainenandSejnowski,1996;Durstewitzetal.,2000;

Hayetal.,2011;AlmogandKorngreen,2014).Themultitudeof typesofvoltage-gatedionchannels,however,representsachal- lengeformodelingofthemembrane-potentialdynamics:unless specificcare istaken, therole ofa specific ion-channelspecies maybeassignedtoanotherion-channelspecieswhenbothcon- ductancesarefittedsimultaneously,i.e.,constrainedbythesame objectivefunctions(Achardetal.,2006).Totacklethisproblem, inKerenetal.(2009),aparameterpeelingexperimentalproce- durewasintroduced,inwhich specifictypesofionchannelsin L5PCareblockedsequentiallyusingdrugs,andtheneuronresponse todifferentstimuliarerecorded ateachstage.Theion-channel conductancesarethenfittedstep-by-steptothesedata.Another strategywasexploredinBahletal.(2012),where experimental datafromL5PCswithandwithoutapicaldendrite(occludedusinga

“pinching”method(BekkersandHäusser,2007))wereusedduring oneofthethreestagesoffitting.Inthiscase,thedatafordifferent stagesoffittingwereobtainedfromseparateexperiments.Both techniquesfacilitatetheoptimizationprocedurebyreducingthe numberoffreeparametersthatarefittedsimultaneously.

Reduced-morphologymodelsmaybecrucialinsimulationsof largenetworksduetothelightercomputationalloadtheyimpose.

Whilethelevelofdetail inthemorphologiesobtainedfrom3D

reconstructionsishigh,theelectrophysiologicalpropertiesofthe distal dendritic segments,as well as the heterogeneity of ion- channelpopulationsbetweendifferentdendriticbranches,remain elusive (Häusser et al., 2000) and are generally not taken into account in the models. However, in many neuron types, den- driticelectrophysiologicalpropertiesvarymonotonicallywiththe distance from the soma(Migliore and Shepherd, 2002), which favorstheuseofsimplified(yetmulti-compartmental)morpholo- gies.Thesesimplifiedmodelsshouldreproducetheexperimentally observedpropertiesofcommunicationbetweenperisomaticand (proximaltomid-distal)dendriticsectionsoftheconsideredneu- ron while reducing the computational load in comparison to full-morphologymodels.

Inthiswork,weusetheexperimentallyvalidatedmodelintro- ducedinHayetal.(2011)(“Hay model”)togeneratesimulated VSDand Ca²⁺-imagingdataaswellassimulated electrophysio- logicalrecordingsinaL5PC.Wesimulatetheparameterpeeling procedurebysequentiallysettingchannelconductancesofdiffer- ention-channelspeciestozerointheHaymodelandmeasuringthe neuronresponsestodifferentstimuliundertheseblockades.We thenfitthechannelconductancesinafour-compartmentmodelto reproducethesedata.Weproposeandapplyafour-stepscheme, wherethethreefirststepsutilizeinformationonvoltageandintra- cellularCa²⁺concentrationsalongthedendriteswithhighspatial resolution.Thefirststepfitstheparametersofreducedmorphol- ogy,inasimilarfashionasinBahletal.(2012),andparameters controllingpassivemembraneproperties.Thesecondstepfitsthe non-specific ion-channelconductances.Theseare importantfor correctlydescribingthedistaldendriticexcitability.Thethirdstep fitstheCa²⁺channelconductancesandSKchannelconductances andtheparametersdescribingCa²⁺dynamics.Thefourthstepfits therestoftheactiveconductances, includingtheconductances responsibleforthespikingbehavior.Weshowthattheobtained reduced-morphologyL5PC model iscost-efficient and faithfully reproducesthemembranedynamicsandspikingbehavior,includ- ingtheBACfiring.Furthermore,wetestourmethodforfittinga neuronmodelwithafull,reconstructedmorphology,andfindthat acceptablefittingresultsareobtainedalsowhenusingthiscomplex dendriticmorphology.

Theobtainedreduced-morphologymodelisespecially useful in network simulationsdue toits lighter requirementsof ran- dom access memory and computation time. We validated our modelbyintroducingit ina biophysicallydetailed L5PCmicro- circuitmodel (Hay andSegev, 2015), which originallyincluded thefull-morphology Haymodel neurons,and showingthat the twomodelsyieldedsimilarnetworkdynamics.Ourcircuitmodel ofreduced-morphologyL5PCspredictsthatinterconnectedL5PCs amplifycertaindelta-range frequenciesdue tothelargecontri- butionsoftheCa²⁺-activatedK⁺currents(SKcurrents)tothecell electrophysiology.

2. Materialsandmethods 2.1. TheL5PCmodel

TheHaymodelofanL5PC,aswellasthereduced-morphology modeldevelopedhere,includethefollowingioniccurrents:fast inactivatingNa⁺current(INat),persistentNa⁺current(INap),non- specific cation current (I_h), muscarinic K⁺ current (I_m), slow inactivatingK⁺current(I_Kp),fastinactivatingK⁺current(I_Kt),fast non-inactivating K⁺ current (IKv3.1), high-voltage-activated Ca²⁺

current(I_CaHVA),low-voltage-activatedCa²⁺current(I_CaLVA),small- conductanceCa²⁺-activatedK⁺current(I_SK),andfinally,thepassive leakcurrent(I_leak).Thecurrentbalanceequationofeachsegment

oftheneuronalmembranecanthusbewrittenas Cm∂V

∂t =I_Nat+I_Nap+I_h+Im+I_Kp+I_Kt+I_Kv3.1 +I_CaHVA+I_CaLVA+I_SK+I_leak+I_axial,

whereeachcurrentspecies,exceptfortheaxialcurrent,isaproduct ofactivationandinactivationvariablesas

I=gm¯ ^Nmh^Nh(E−V).

Here, ¯gisthemaximalconductanceoftheionchannels,mandhare theactivationandinactivationvariables,NmandN_hareconstants describingthegatingmechanismsofthechannel,andEistherever- salpotentialcorrespondingtotheionicspecies.Reversalpotentials ofNa⁺andK⁺areconstants(E_Na=50mV,E_K=−85mV),whilethe reversalpotentialofCa²⁺isdeterminedthroughtheNernstequa- tionbytheintracellular[Ca²⁺](typically,valuesofECaintheHay modelvarybetween96and120mV(Mäki-Marttunenetal.,2016)).

TheaxialcurrentI_axialisdeterminedbytheaxialresistanceandthe voltagedifferencebetweentheconsideredmembranesegmentand itsneighbors.Fordetailsonthemodelequations,seeHayetal.

(2011).

Theintracellular[Ca²⁺]obeysthefollowingdynamics:

d[Ca²⁺]

dt =ICaHVA+ICaLVA

2Fd −[Ca²⁺]−c_{min decay} ,

where I_CaHVA and I_CaLVA arethe highand low-voltage activated Ca²⁺currentsenteringtheconsideredcellsegment,represents thefractionofCa²⁺ionsentering thecellthatcontributetothe intracellular[Ca²⁺],FtheFaradayconstant,disthedepthofthe sub-membranelayerconsideredforcalculationofconcentration, c_min therestingintracellular[Ca²⁺],and_decayisthedecaytime constantoftheintracellular[Ca²⁺].

Inthiswork,werunparameter-fittingtaskswhereweassume theion-channeldynamicsfixed,andonlyvarytheparametersgov- erningthemaximalconductances,namely,g_Nat,g_Nap,g_h,gm,g_Kp, g_Kt,g_Kv3.1,g_CaHVA,g_CaLVA,g_SK,andg_l,andparameters and_decay thatcontroltheCa²⁺dynamics.

2.1.1. Synapsemodel

Themodelforsynapticcurrents,whichareusedinsimulations includinginvivo-likebackgroundsynapticfiring,isadoptedfrom HayandSegev(2015).EachexcitatorysynapseconductsAMPA- andNMDA-mediatedcurrents,andeachinhibitorysynapsecon- ductsGABA_A-mediated currents.These aremodeledas follows:

I_AMPA=w_AMPAg_Glu(B_AMPA−A_AMPA)(E_Glu−V) I_NMDA=w_NMDAg_Glu(B_NMDA−A_NMDA)(E_Glu−V)c_Mg2+

I_GABA=w_GABAg_GABA(B_GABA−A_GABA)(E_GABA−V).

Here, w_AMPA, w_NMDA, and w_GABA are synaptic weights, g_Glu=0.0004␮S is the baseline conductance of a single gluta- matergic synapse and g_GABA=0.001␮S that of an inhibitory synapse.VariablesB_AMPA,A_AMPA,B_NMDA,A_NMDA,B_GABA,andA_GABA areincreasedbyapositive constant(chosensuchthatthepeak conductance of the synaptic current after a long silent period wouldbeg_Gluorg_GABA)ateach timeofsynapticactivation,and otherwisedecaytowardszero,AfasterthanB,asfollows:

dA dt =−A

A

and dB dt =−B

B

. (1)

Timeconstants for the differentsynaptic species are _A,AMPA = 0.2ms,_B,AMPA=1.7ms,_A,NMDA=0.29ms,_B,NMDA=43ms,_A,GABA

=0.2ms,andB,GABA=1.7ms.Variablec_Mg2+representstheMg²⁺

blockoftheNMDA-activatedchannel,anditsvalueisdetermined bythemembranepotential(JahrandStevens,1990).Parameters E_Glu=0mVandEGABA=−80mVarethereversalpotentialsofthe glutamatergicandinhibitorysynapse,respectively.

ThebackgroundsynapsesareactivatedatPoisson-distributed activationtimes,whiletheintra-networksynapsesthatareexclu- sively glutamatergic (when present) are activated immediately aftera spike inthepre-synaptic neuron. Theactivations ofthe synapsesfollowaprobabilisticrule(Ramaswamyetal.,2012)that allowsthemodelingofbothshort-termdepressionandfacilitation.

Fordetails,seeHayandSegev(2015).

2.1.2. Modelimplementation

Throughout the work, the NEURON software (Hines and Carnevale,1997)isusedforsimulatingtheL5PC.Forsingle-cellsim- ulations,theadaptivetime-stepintegrationmethodisused,while for themulti-neuron simulations,the fixedtime step 0.025ms is used. Our NEURON and Python scriptsboth for runningthe parameterfittingsandsimulatingthefittedmodelsarepublicly availableat:https://senselab.med.yale.edu/ModelDB/showModel.

cshtml?model=187474 – in addition, theprincipal model for a singlereduced-morphologyL5PC isdescribed using NeuroML-2 (Cannonetal.,2014).TheNEURONscriptsarebasedonthepublicly availablemodelspublishedinHayetal.(2011)andHayandSegev (2015).WhenimplementingthemodelofHayandSegev(2015) usingourreduced-morphologyL5PCs,weupdatedthebackground synapsemodelasfollows.Wegroupedthebackgroundsynapses thatwerelocatedin thesamesegmentintoone synapsegroup, wheretheindividualsynapsessharedthevariablesAandB(see Eq.(1))–onlytheactivationtimesofeachindividualsynapsewere savedintothememory.Thisradicallydecreasedthecomputational loadimposedbysolvingthedifferentialequations(Eq.(1)).

2.2. Stepwisefittingprocedure

Wepresentaflexiblestepwiseneuronmodel-fittingframework andapplyittofitafour-compartmentmodeltosimulatedelectro- physiologyandimagingdatafromanL5PC.Thefittedparameters arelistedinTable1.Theparametervaluerangesaretakenfrom Table2inHayetal.(2011),withcertainexceptions(seeSupple- mentarymaterial,SectionS1).Duringeachstep,weuseaPython implementation(publishedbytheauthorsofBahletal.,2012)ofthe non-dominatedsortinggeneticalgorithmII(NSGA-II)(Debetal., 2002)fortheparameteroptimization.Thecrossoverandmutation parametersandtheprobabilityofmutationperparameterarekept fixedasc=20,m=20and pm=0.5.Thethreefirstfittingsteps are performedusing Nsamp=1000samples and Ngen=20 genera- tions,thefourthfittingstepusingNsamp=2000samplesandNgen=20 generations–thesevalueswerefoundadequateforobtainingan acceptablefittothedata,althoughnoconvergencetoalocalopti- mumisguaranteed.Thecapacityofthegeneticalgorithmisdouble thepopulationsize.

TheobjectivesaregiveninTable2–seeSection2.2.5fordetails oftheusedobjectivefunctions.Duringeachstep,theparameters thatwereoptimizedduringthepreviousstepsarekeptfixed,while theparametersthatwillbeoptimizedinthefollowingstepsare settozero(simulatingaperfectblockadeofthecorrespondingion channels).

2.2.1. Firststep:Morphology

Wefollowtheprocedureofthefirstfittingstepaspresentedin Bahletal.(2012),withcertainchanges.InBahletal.(2012),theleak conductancesweresettofixedvaluesbothinthetargetmodeland inthereduced-morphologymodelunderoptimization.Here,we donotchangetheleakconductanceinthetargetmodelbutkeepit

Table1

Thetableofparametersusedineachstepandtheirboundariesandreferencevalues.

Thefirstcolumnshowsthestepinwhichtheparameterisfitted.Thesecondcolumn namestheparameters,andthethirdandfourthcolumnshowthevaluesusedby theoptimizationalgorithmasthelowerandupperlimits,respectively.Thefifth columnshowsthecorrespondingvaluesinthefullmodel–notethatinthefull model,theapicaldendritictreeisdividedintotwoequallylongsections,andincase theparametersvaryspatially,thevaluesshownaretheminimumandmaximumof thevaluesacrossthesesections.ConductancesaregiveninS/cm²,lengthsin␮m, axialresistancesincm,andcapacitancesin␮F/cm².

Step Parameter Lower Upper Fullmodel

1 L^soma 11.58 46.34 23.17

1 L^basal 141.06 564.26 282.12

1 L^apic 325.0 1300.0 650

1 L^tuft 325.26 1301.06 650.53

1 R^soma_a 20 500 100

1 R^basal_a 10 1000 100

1 R^apic_a 10 1000 100

1 R^tuft_a 10 1000 100

1 c^somam 0.5 2.0 1.0

1 c^basal_m 0.5 4.0 2.0

1 c^apicm 0.5 4.0 2.0

1 c^tuft_m 0.5 4.0 2.0

1 g^soma_l 2×10⁻⁵ 0.0001 3.38×10⁻⁵

1 g^basal_l 1.5×10⁻⁵ 0.0001 4.67×10⁻⁵

1 g^apic_l 1.5×10⁻⁵ 0.0001 5.89×10⁻⁵

1 g^tuft_l 1.5×10⁻⁵ 0.0001 5.89×10⁻⁵

2 Eh −55 −35 −45

2 g^soma_h 0 0.0008 1.29×10⁻⁶

2 g^basal_h 0 0.0008 1.30×10⁻⁶–1.71×10⁻⁶

2 g^apic_h 0 0.008 1.78×10⁻⁶–0.000127

2 g^soma_l 2×10⁻⁵ 0.0001 3.38×10⁻⁵

2 g^basal_l 1.5×10⁻⁵ 0.0001 4.67×10⁻⁵

2 g^apic_l 1.5×10⁻⁵ 0.0001 5.89×10⁻⁵

2 g^tuft_l 1.5×10⁻⁵ 0.0001 5.89×10⁻⁵

3 g^soma_CaHVA 0 0.001 3.66×10⁻¹⁰

3 g^apic_CaHVA 0 0.0025 5.55×10⁻⁵

3 g^tuft_CaHVA 0 0.025 5.55×10⁻⁵–0.000555

3 g^soma_CaLVA 0 0.01 3.86×10⁻⁸

3 g^apic_CaLVA 0 0.1 0.000187

3 g^tuft_CaLVA 0 1.0 0.000187–0.0187

3 ^soma 0.0005 0.05 0.000501

3 ^apic 0.0005 0.05 0.000509

3 ^tuft 0.0005 0.05 0.000509

3 decay^soma 20.0 1000.0 460.0

3 decay^apic 20.0 200.0 122.0

3 decay^tuft 20.0 200.0 122.0

3 g^soma_SK 0 0.1 4.18×10⁻⁵

3 g^apic_SK 0 0.005 0.0012

3 g^tuft_SK 0 0.005 0.0012

4 g^soma_Nat 0 4.0 2.04

4 g^soma_Nap 0 0.01 0.00172

4 g^soma_Kt 0 0.1 0.0812

4 g^soma_Kp 0 1.0 0.00223

4 g^soma_Kv3.1 0 2.0 0.693

4 g^apic_m 0 0.0005 6.75×10⁻⁵

4 g^apic_Nat 0 0.02 0.0213

4 g^apic_Kv3.1 0 0.02 0.000261

4 g^tuft_m 0 0.0005 6.75×10⁻⁵

4 g^tuft_Nat 0 0.02 0.0213

4 g^tuft_Kv3.1 0 0.02 0.000261

fixed.However,similarlyasinBahletal.(2012),wefittheleakcon- ductancesagaininthesecondstep.Thus,theoutputparametersof ourfirststepcanbeinterpretedasoptimallengths,axialresistances andmembranecapacitancesforwhichtheconsideredobjectivescan bewellmetwithsomevaluestheofleakconductances.IntheSup- plementarymaterial,SectionS7,weshowthatthemethodused inBahletal.(2012)produces validfittingresultsin ourframe- workaswell,andinSectionS9,weshowthatthefirststepcanbe combinedwiththesecondonewhenfittinganeuronmodelwith reconstructedmorphology.

FollowingBahletal.(2012),weresetthediametersofeachcom- partmentduringeachiterationoftheoptimizationalgorithmsuch thatthemembraneareaofthecompartmentisequaltothetotal membraneareaofthecorrespondingsegmentsofthefullmodel.

The objectivefunctionsare chosen suchthattheresponse of a reduced-morphologymodelneuronwithoptimalparameters to somaticand apicalinputs isassimilaraspossibletothecorre- spondingresponseinthefull-morphologyneuron.Theaccuracy ofthemodelneuronresponsetosomaticinputs ismeasuredin termsofbothsomaticmembrane-potentialtimeseries(objective 1.3)anddendriticsteady-statevoltagedistribution(objective1.1), whiletheaccuracyoftheresponsetoapicalinputsismeasuredonly inthelatterterms(objective1.2).

2.2.2. Secondstep:Passiveand“nearlypassive”properties

After fixing the morphological parameters (compartment lengths,axialresistancesandmembranecapacitances),weopti- mize therestof theparameters that are majorcontributors to membrane response propertiesat rest.These arethe leakcon- ductances, non-specific ion-channel conductances (g_h), and the reversalpotentialofthenon-specificioncurrent(E_h).Inthesame wayasinthefirststep,wesettheobjectivessuchthataneuron modelwithoptimalparametervalueswouldrespondaccurately tosomaticinputbothintermsofmembrane-potentialtimeseries (objective2.1)anditssteady-state distributionacrossdendrites (objective 2.2). We assume the leak reversal potential known (−90mV),butifitisunknown,itshouldbefittedhereaswell.

2.2.3. Thirdstep:Ca²⁺dynamicsandSKcurrents

Oncethepassivepropertiesandnon-specificionchannelprop- ertieshavebeenoptimized,wesearchfortheoptimalparameters thatgovernCa²⁺currents,Ca²⁺dynamics,andCa²⁺-dependentK⁺ currents.TheseparametersareintheHaymodelthefollowing:

high-voltageactivated(HVA)andlow-voltageactivated(LVA)Ca²⁺

channelconductances(g_CaHVAandg_CaLVA),percentageofCa²⁺cur- rentinclusionintosub-membranespace(),timeconstantofdecay offreeCa²⁺(_decay),andtheSKchannelconductance(gSK).Itmay beimportantthattheseparametersbefittedsimultaneously,as theCa²⁺currentshaveastrongfeedbackeffectonthemembrane potentialthroughtheSKchannels.Forthisreason,weuseobjec- tivefunctionsthatpromote both a correctmembrane-potential response and an accurate [Ca²⁺] response to different stimuli.

A goodfit postulatesthat thesteady-state membrane potential (objective3.2)andCa²⁺concentration(objective3.3)distributions acrossthedendrites,asaresponsetosomaticDC,aresimilarto thoseinthetargetneuron,andthatthemaximalmembranepoten- tial(objective3.4)andCa²⁺concentration(objective3.5)responses toanEPSP-likecurrentinjectionareaccurateaswell.Thedistribu- tionofmembranepotentialismeasuredacrossthewholeneuron, whiletheCa²⁺concentrationonlyneedstobemeasuredattheapi- caldendrite:asthesignalpropagationinthebasaldendritesof anL5PCwasnearlypassive(Nevianetal.,2007),theHaymodel doesnotincludeCa²⁺channelsinthebasaldendriticcompartments (Hayetal.,2011).Inaddition,themembrane-potentialtimeseries atsomaduringa100-msDCpulseshouldbeascloseaspossible totheoneinthetargetneuron(objective3.1).Thecorrecttempo- ralactivationanddeactivationoftheCa²⁺andSKcurrentsatthe somaisaprerequisiteforanaccuratespikingbehavior,especially regardingthemediumafterhyperpolarization(AHP)period.

2.2.4. Fourthstep:Correctspikingbehavior

AfterfixingtheparametersgoverningCa²⁺dynamicsandSKcur- rents,weoptimizetherestoftheparameterstomakethemodel produceacceptablespikingbehavior.Werequirethatthef–Icurve beclosetothatinthetargetneuron(objective4.4),andthatthe somaticmembrane-potentialtimeseriesbesimilartothatinthe

Table2

Objectivefunctionsusedineachstep.Thefirstcolumnnamestheobjectivessuchthatthenumberleftofthedotindicatesthestepandthenumberrightofthedotindicates theobjectivenumber.Thesecondcolumndescribesthetypeofobjectivefunction,thethirdcolumnindicatesthesiteofrecordingintheneuron,andthefourthcolumn describesthetypeofstimulusused.Thefifthcolumnreferstotheequationcorrespondingtothetypeofobjectivefunction(seeSection2.2.5).Thesixthcolumnliststhe amplitudesofthestimuli:eachvaluerepresentsanindividualrunwhereonlytheshownstimuluswiththeshownamplitudeisappliedtotheneuron,exceptforobjective 4.3,whereacombinationoftwostimuliwitha5msinter-stimulusintervalisapplied.

Objectivenumber Objectivefunction Wheremeasured Stimulustype Eq. Stimulus

amplitude 1.1 Differenceindistributionofsteady-state

membranepotential

Dendrites 3000-msDCpulseatsoma (5) 0.5nA

1.2 Differenceindistributionofpeak membranepotential

Dendrites EPSP-likecurrentinjection attheapicaldendrite 620␮mfromsoma

(5) 0.5nA

1.3 Differenceinmembranepotentialtime series

Soma 100-msDCpulseatsoma (2) −1.0nA

0.0nA 1.0nA 2.1 Differenceindistributionofsteady-state

membranepotential

Dendrites 3000-msDCpulseatsoma (5) 0nA

0.5nA 1.0nA 2.2 Differenceinmembranepotentialtime

series

Soma 100-msDCpulseatsoma (2) 0.5nA

1.0nA 3.1 Differenceinmembranepotentialtime

series

Soma 100-msDCpulseatsoma (2) 1.0nA

2.0nA 3.2 Differenceindistributionofsteady-state

membranepotentials

Dendrites 3000-msDCpulseatsoma (5) 0.5nA

1.0nA 3.3 Differenceindistributionofsteady-state

intracellular[Ca²⁺]

Apicaldendrite 3000-msDCpulseatsoma (5) 0.5nA 1.0nA 3.4 Differenceindistributionofpeak

membranepotential

Dendrites EPSP-likecurrentinjection attheapicaldendrite620

␮mfromsoma

(5) 0.5nA 1.0nA 3.5 Differenceindistributionofpeak

intracellular[Ca²⁺]

Apicaldendrite EPSP-likecurrentinjection attheapicaldendrite620

␮mfromsoma

(5) 0.5nA 1.0nA

4.1 Differenceinmembranepotentialtime seriesandnumbersofspikes

Soma 100-msDCpulseatsoma (4) 0.25nA

0.5nA 4.2 Differenceinmembranepotentialtime

seriesandnumbersofspikes

Soma 5-msDCatsoma (4) 1.9nA

4.3 Differenceinmembranepotentialtime seriesandnumbersofspikes

Soma 5-msDCatsomaand

EPSP-likecurrentatapical dendrite

(4) 1.9nA(soma)+ 0.5nA(apical)

4.4 Differenceinnumbersofspikes Soma 3000-msDCatsoma (3) 0.78nA

1.0nA 1.9nA

targetneuron,bothwhengivensomaticsub-thresholdandsupra- thresholdDC(objective4.1and4.2)andacombinationofsomatic andapicalstimulithatinducesBACfiring(objective4.3).Thisstepis computationallythemostchallengingone,bothbecauseofthecost ofevaluatingthef–Icurves(albeitheredoneforonlythreevalues ofcurrentamplitude),andbecauseofthelargegeneticpopulation thatisneededtosecurethattheobjectivefunctionsbemetclosely enough.Itmightbeadvisabletodividethisstepfurtherintotwo steps,where,e.g.,firsttheslowerion-channelconductances(Im, IKp,INap)areoptimized,andfinallytheconductancesofthefaster ones(I_Kt,I_Nat,I_Kv3.1).Thisisleftforfuturestudies.

2.2.5. Distancemetricsoftheobjectivefunctions

TheobjectivefunctionsofTable2aredesignedtocapturecor- rectmembrane-potentialbehavioracrossthespatialextentofthe neuron.Certainobjective functionsonly considerthequantities measuredat soma(objectives 1.3, 2.2, 3.1, and 4.1–4.4). These objectivesarefurthercategorizedtothosethataimatcapturing thecorrecttimeseries(objectives1.3,2.2,3.1,4.1),thosethatonly aimatcapturingthecorrectnumbersofspikes(objective4.4),and thosethataimatcapturingboth(objectives4.2and4.3).Thediffer-

enceintimeseriesisquantifiedusingtheL¹norm(meanabsolute difference)betweenthetargetandcandidatemembranepotential (afunctionofparametersp)acrossatimewindowrangingfrom 50msbeforethestartofthestimulusto200msafterthestartof thestimulus:

f(p)= 1 250ms·1mV

−50 ms

|V(p,t)−Vtarget(t)|dt. (2)

TheL¹normisrelativelylessstrictagainstdifferencesinspiketim- ingsbetweenthetargetandcandidatedata(butelaboratesmore thebreadthofthetimewindowinwhichdifferencesinmembrane potentialsoccur)thantheL²norm.Thedifferenceinf–Icurvesis quantifiedbythe2-norm,i.e.,

f(p)=

|N_spikes,Ii(p)−Nspikes,target,I_i|² (3)

whennoothermeasuresareconsidered(objective4.4).However, objectives4.2and4.3useacombinationofspikenumberandtim- ingdata and thetime coursedata– in theseobjectivesallthe differencesarequantifiedusing1-norms:

f(p)=a1

−50 ms

|V(p,t)−Vtarget(t)|dt+a2

min(N_spikes(p),N

|t_spike⁽ⁱ⁾ (p)−tspike,target⁽ⁱ⁾ |+|N_spikes(p)−Nspikes,target|. (4)

We chose the coefficients as a₁= _250ms¹_·12mV and a₂= ₂₀¹_ms, whichmeansthatanaveragedifferenceof12mVispenalizedas muchasthesummeddistanceof20msbetweenthespiketim- ings,andfurthermoreasmuchasadifferenceofonespikeinspike counts.Thesevaluesgavetheerrorfunctionadesiredempiricalbal- ance,suchthatvoltagetracesthatlooked(byeye)moredifferent fromthetargetdatathanothersalsoreceivedlargererrorvalues.

Therestoftheobjectives,namely,objectives1.1–1.2,2.1,and 3.2–3.5,concernquantitiesmeasuredalongthedendrites–both neartoand farfromsoma. Fortheseobjectives,themembrane potentialsV(eitherthesteady-statemembranepotentialfollow- ingalongstimulusasinobjectives1.1,2.1,and3.2,ormaximal membranepotentialduringapulsestimulusasinobjectives1.2 and 3.4) or Ca²⁺ concentrations c (objectives 3.3 and 3.5) are quantified across the spatial extent. This is done using either 20 (in simulations of reduced-morphologyL5PCs when synap- tic backgroundinputsare notmodeled)or 5(in simulationsof reconstructed-morphologyL5PCsandreduced-morphologyL5PCs withspontaneoussynaptic inputs)recordingsitespercompart- ment.Withfourcompartmentsinthereduced-morphologyneuron and 193 in the full-morphology neuron, this implies that the total numbers of recording sites are N_recs^(reduced)=80 (for objec- tives1.1–3.5)andN^(full)_recs =965.Foreachrecordingsite,aspatial coordinatedis determinedbyitsdistancefromsoma;negative valuesof dareassigned torecording sitesalongthebasal den- driteandpositivevaluesalongtheapicaldendrite.Thedifference ofthese2-dimensionaldata,(d⁽ⁱ⁾(p),V⁽ⁱ⁾(p))i=1,...,N_recs^(reduced)and (d⁽ⁱ⁾_target,V_target⁽ⁱ⁾ )i=1,...,N^(full)_recs is quantified asfollows. First, we disregardthedatacorrespondingtothoserecordingsitesatthe reducedmorphologythat arefurtheraway fromthesomathan thefarthestsitesinthereconstructed morphology.Inourwork, thisisdoneforrecordingsitesforwhichd<−282␮mord>1301

␮m (i.e. d outside the spatial extent of cell #1 in Hay et al.

(2011)).Second,wenormalizetheremainingdatabythemaximal range(max_i(d⁽ⁱ⁾_target)−min_i(d⁽ⁱ⁾_target)ormax_i(V_target⁽ⁱ⁾ )−min_i(V_target⁽ⁱ⁾ )), andhencetheresultingdataareinR² plane,wheretheaverage distanceofthereduced-morphology-neurondatafromtheirnearest neighborinthereconstructed-morphology-neurondataisusedasthe distancemetric.Inmathematicalterms,wecanthuswrite

f(d, V, dtarget, Vtarget)=

1 Nacc

N^(reduced)_recs

i=1 d⁽ⁱ⁾≥−282 and

d⁽ⁱ⁾≤−1301 min^N

(full) recs j=1

⎛

⎜ ⎝

d⁽ⁱ⁾−d^(j)_target max_k(d^(k)_target)−min_k(d^(k)_target)

+

V⁽ⁱ⁾−V_target^(j) max_k(V_target^(k) )−min_k(V_target^(k) )

⎞

⎟ ⎠

(5)

whereNaccisthenumberofaccepteddatapoints,i.e.,datapointsfor whichd⁽ⁱ⁾≥−282␮mandd⁽ⁱ⁾≤1301␮m.Notethatthelimitation oftheaccepteddatapointsallowsneuronsthataresignificantly longertobecreated.However,asthediametersofthesegments arerestrictedbytherulethatconservesthetotalmembranearea, neuronswithtoolongdendritesenduphavingtoothindiameters, whichislikelytopreventagoodfittothedata.

2.2.6. Combiningthesteps

Asweuseageneticmulti-objectiveoptimizationalgorithm,at theendofeachstepwehaveapopulationofPareto-efficient(see, e.g.,VanGeitetal.,2008)parametersets.Weapplyanexploratory scheme, where the best parameter sets of each objective and parametersetsthatperformwellintwoseparateobjectivesare handedontothenextstep.Intheinterestofreducingthenum- beroftheparametersetsthatarefixedduringtheearlysteps,we

reducethenumberofobjectivefunctionsbygroupingsomeofthem together.

Firstly,mostoftheobjectivesofTable2consistofstimuliof differentamplitudes(objectives1.3–4.1and4.4).Inthesecases, theobjectivefunctionsaredefinedassumsofthesub-objective functions.Asanexample,forobjective1.3,wehave

f_1.3(p)=f_{1.3,−1.0 nA}(p)+f_{1.3,0.0 nA}(p)+f_{1.3,1.0 nA}(p), wherefunctionsf_1.3,I(p)areoftheformofEq.(2).

Secondly, we further combine the objectives for cor- rect membrane-potential distribution across the dendrites as a response to somatic DC and that as a response to EPSP-like stimulus. Namely, we group together objectives 1.1 and 1.2 as f1.1+1.2(p)=f1.1(p)+5f1.2(p), objectives 3.2 and 3.4 as f_3.2+3.4(p)=f_3.2(p)+f_3.4(p), objectives 3.3 and 3.5 as f_3.3+3.5(p)=f_3.3(p)+f_3.5(p), and objectives 4.2 and 4.3 as f4.2+4.3(p)=f4.2(p)+f4.3(p). The factor 5 is chosen for f1.2(p) due to the small effect of EPSP-like stimulus on dendritic peak membrane potentialswhen allactiveconductancesare blocked comparedtothatofthesomaticDC;inotherobjectivestheeffects are approximately of the same order of magnitude (data not shown).

Dueto this grouping,themulti-objective optimizationalgo- rithmisgiventwoobjectivefunctionsforthefirstandsecondsteps, andthreeforthethirdandfourthsteps.Inpractice,thismeansthat thefirststepoptimizationisperformedonce,andthreecandidate parametersetsareobtained–onethat performsbestinf_1.1+1.2, onethatperformsbestinf_1.3,andonethatperformswellinboth.

Whenpickingtheparametersetthatproducesagoodfittotwo objectivefunctions,wefirstnormalizetheerrorfunctionvalues ofbothobjectivesbytheirmedians(acrossthewholepopulation ofparametersattheendoftheoptimization)andthensumthem together:theparametersetpithatproducethesmallestvalueof thesum

f_1.1+1.2(p_i)

median_j(f_1.1+1.2(p_j))+ f_1.3(p_i)

median_j(f_1.3(p_j)) (6) ischosen.Thesecondstepoptimizationisthenperformedthree times(onceusingeachofthesethreeparametersets)andninecan- didateparametersetsareobtained.Thethirdstepoptimizationis

thusperformedninetimes,andeachoptimizationgivessixcandi- dateparametersets(onethatperformsbestinf_3.1,onethatisbest inf3.2+3.4,onethatisbestinf3.3+3.5,andthreeintermediateonesthat performwellintwoofthethreeobjectivefunctions),andtherefore, thefinalstepisperformedforamaximumnumberof54param- etersets.Notethatsomeoftheoptimizationsperformedduring thefirstthreestepsmaygiveparametersetsthatdonotproducea goodfittothedata.Insuchcases,onlythefeasibleparametersets arehandedovertothefinalfittingstep.

2.3. Powerspectra

We illustrate and quantify the possible oscillations in the neuronalnetworkdynamicsusingthepowerspectraofthepop- ulationspiketrains.Thepowerspectrum ofa spiketrains(t)=

j=1 ıt_j(t),wherevariablest_jrepresentthespiketimes,isdeter- minedas

Ps(f)=|Fs(f)|²,

whereFs(f)istheFouriercomponentforthefrequencyf.Thiscom- ponentcanbedeterminedas

Fs(f)=

−∞

s(t)e^−2itfdt=

N

j=1

e^−2itjf,

whereiistheimaginaryunit.

3. Results

3.1. Morphologyparametersandion-channelconductancesfor thereducedmodelcanbefittedtothefullmodeldata

Weappliedthestepwisemodel-fittingprocedureusingsim- ulated data(obtained fromsimulations of theHaymodel with reconstructed morphology)as the targetdata for the objective functions.We used the multi-objective optimizationalgorithm developedinDebetal.(2002)tofindtheoptimalparameterval- ues,inasimilarfashionasdoneinBahletal.(2012).Figs.1–4show theperformanceofthefittedmodelinfulfillingtheobjectives,and Table3liststheobtainedparametervalues.

Thereducedmodelshowsagoodfittomostofthemeasured quantities.However,thethirdsteprevealsthatallaspectsoftheion channeldistributionscannotbeaccuratelypreservedinthefour- compartmentmodel:theresponsesofthefullmodelwerebest reproducedbylettingtheCa²⁺channelconductancesgotozeroin theapicaltrunk,whilekeepingthecorrespondingvaluesattheapi- caltuftrelativelylarge.Bycontrastintheoriginalmodel,theCa²⁺

channelconductanceswerenon-zeroallalongtheapicaldendrite, buthad extremely largevalues at segments685–885␮maway fromthesoma(inthe“hotzoneofCa²⁺channels”).Fig.S1(Sup- plementarymaterialSectionS2)showsthattheerrorfunctionsfor thespatialdistributionsofmembranepotentialandCa²⁺concen- trationinthethirdstepfittingcouldbedecreasedbyincludingan extracompartmentthatrepresentedthehotzoneinthereduced- morphologyneuron.Nevertheless,thiswouldmakethemodela six-compartmentmodel,asthefurthestcompartmentwouldhave tobedividedtothreecompartments.Furthermore,inarecentstudy (AlmogandKorngreen,2014),theCa²⁺channelswerebetterfitted byusinglinearincreasesinCa²⁺channeldensitiesthanaGaussian- shapeddistributionwithlargestconductances aroundthemain bifurcationpointalongtheapicaldendrite.Therefore,weapplied theparametersobtainedfromthefittingofthefour-compartment model(Fig.3)intherestofthiswork.

Fig.5showstheevolutionoftheobjectivefunctionsacrossthe generationsinallfourstepsoffitting:inmostcasesthevaluesof theerrorfunctiondramaticallydecreasedduringthefirstfivetoten generations,afterwhichonlymodestimprovementwasachieved.

FortheresultsshowninFigs.1–4,onlythefinalparametersetwas used,buttomakesurethefittingprocedureisrobust,werepeated thefinalfittingstep10times.Allofthese10samplesshowedthe correctnumbersofspikesinresponsetothestimuliofFig.4B–D (datanotshown).

Thestepwiseapproachgainsadvantagefromthefactthatthe parametersearchspaceissmallerthanwhenallparametersare fittedatonce.In Fig.6,we showthatfittingallmodelparame- terssimultaneously(withoutsimulationofionchannelblockades) didnotproduceasgoodfittingresultsasthestepwisemethod:

outoftentrials,onlyoneoftheobtainedparametersetsproduced thecorrectnumberofspikesasresponsetothestimuliofobjec- tives4.1–4.3.Thisparametersetwascharacterizedbyrelatively

Table3

Parametervaluesobtainedfromthemulti-objectiveoptimizationsofFigs.1–4.Note thatthepassiveleakconductancesareinitiallyfittedatthefirststep,butrefittedat thesecondstep.ConductancesaregiveninS/cm²,lengthsin␮m,axialresistances incm,andcapacitancesin␮F/cm².Thefirststepparametersetwastheonethat minimizedf1.3,whilethesecondstepparametersetminimizedf2.1.Thethirdstep parametersminimizedthesumoff3.2+3.4andf3.3+3.5andthefourthstepparameters thesumoff4.2+4.3andf4.4,wherethefunctionvalueswerenormalizedbythemedians ofthecorrespondingfunctionsacrossthegeneticalgorithmpopulation(e.g.,seeEq.

(6)).

STEP1

Variable Value

L^soma 24.5

L^basal 426

L^apic 400

L^tuft 702

R^soma_a 380

R^basal_a 197

R^apic_a 958

R^tufta 224

c_m^soma 1.22

c_m^basal 1.94

c_m^apic 1.45

c_m^tuft 2.6

g^soma_l 7.8·10⁻⁵

g^basal_l 2.56·10⁻⁵

g^apic_l 5.92·10⁻⁵

g^tuft_l 6.75·10⁻⁵

STEP2

Variable Value

Eh −40.7

g^soma_h 0.000279

g^basal_h 0.000294

g^apic_h 0

g^tuft_h 0.00493

g^soma_l 4.37·10⁻⁵

g^basal_l 3.79·10⁻⁵

g^apic_l 5.29·10⁻⁵

g^tuft_l 6.83·10⁻⁵

STEP3

Variable Value

g^soma_CaHVA 0.000838

g^apic_CaHVA 0

g^tuft_CaHVA 0.000977

g^soma_CaLVA 0.00311

g^apic_CaLVA 0

g^tuft_CaLVA 0.000487

^soma 0.0005

^apic 0.0347

^tuft 0.0005

decaysoma 488

decayapic 142

decaytuft 95.4

g^soma_SK 0.0479

g^apic_SK 0.000231

g^tuft_SK 0.00365

STEP4

Variable Value

g^soma_Nat 2.41

g^soma_Nap 0.00206

g^soma_Kt 0.0239

g^soma_Kp 0.000176

g^soma_Kv3.1 0.701

g^apic_m 0.000143

g^apic_Nat 0.0135

g^apic_Kv3.1 0.00121

g^tuft_m 0.000113

g^tuft_Nat 0.0131

g^tuft_Kv3.1 0