Layer V pyramidal cell model with reduced morphology (Mäki-Marttunen et al 2018)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:187474
" ... In this work, we develop and apply an automated, stepwise method for fitting a neuron model to data with fine spatial resolution, such as that achievable with voltage sensitive dyes (VSDs) and Ca2+ imaging. ... We apply our method to simulated data from layer 5 pyramidal cells (L5PCs) and construct a model with reduced neuronal morphology. We connect the reduced-morphology neurons into a network and validate against simulated data from a high-resolution L5PC network model. ..."
References:
1 . Hay E, Hill S, Schürmann F, Markram H, Segev I (2011) Models of neocortical layer 5b pyramidal cells capturing a wide range of dendritic and perisomatic active properties. PLoS Comput Biol 7:e1002107 [PubMed]
2 . Hay E, Segev I (2015) Dendritic Excitability and Gain Control in Recurrent Cortical Microcircuits. Cereb Cortex 25:3561-71 [PubMed]
3 . Mäki-Marttunen T, Halnes G, Devor A, Metzner C, Dale AM, Andreassen OA, Einevoll GT (2018) A stepwise neuron model fitting procedure designed for recordings with high spatial resolution: Application to layer 5 pyramidal cells. J Neurosci Methods 293:264-283 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism: Neocortex;
Cell Type(s): Neocortex L5/6 pyramidal GLU cell;
Channel(s): I Na,p; I Na,t; I L high threshold; I T low threshold; I A; I M; I h; I K,Ca; I Calcium; I A, slow;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON; NEURON (web link to model); Python; NeuroML;
Model Concept(s):
Implementer(s): Maki-Marttunen, Tuomo [tuomo.maki-marttunen at tut.fi]; Metzner, Christoph [c.metzner at herts.ac.uk];
Search NeuronDB for information about:  Neocortex L5/6 pyramidal GLU cell; I Na,p; I Na,t; I L high threshold; I T low threshold; I A; I M; I h; I K,Ca; I Calcium; I A, slow;
/
reducedhaymodel
snmf
models
morphologies
README.html
Ca_HVA.mod *
Ca_LVAst.mod *
CaDynamics_E2.mod *
epsp.mod *
Ih.mod *
Im.mod *
K_Pst.mod *
K_Tst.mod *
Nap_Et2.mod *
NaTa_t.mod *
SK_E2.mod *
SKv3_1.mod *
drawfits_withmids_combfs.py
drawfits_withmids_combfs_nseg5.py
drawmorphs.py
drawobjective_evolution.py
emoo.py
mytools.py
originalrun.sav
pars_withmids_combfs_1a.sav
pars_withmids_combfs_1a_0a.sav
pars_withmids_combfs_1a_0a_5a.sav
pars_withmids_combfs_1a_0a_5a_5a.sav
pars_withmids_combfs_final.sav *
pars_withmids_combfs_fixed_final.sav
screenshot.png
snmf_protocols.py
snmf_target.py
snmf_withmids_combfs.py
                            
#Target data generation for the simpilified neuron model fitter.
#Loads the full Hay neuron model (Hay et al. 2011 Models of Neocortical Layer 5b Pyramidal Cells Capturing a Wide Range of 
#  Dendritic and Perisomatic Active Properties), simulates the response to the stimuli defined in snmf_protocols.py,
#  and saves the quantities defined in snmf_protocols.py
#hoc-code based on the NEURON implementation of the Hay model, publicly available in https://senselab.med.yale.edu/modeldb/showModel.cshtml?model=139653
#Tuomo Maki-Marttunen, 2015-2016
#
#Running the script:
#  python snmf_target.py
#
#Files needed:
#  snmf_protocols.py:                                                        
#    Fitting protocols (which quantities are fitted at each step; what stimuli are for each objective function; etc.)
#  mytools.py:                                                                                                       
#    General tools for e.g. spike detection                                                                          
#                                                                                                                      
#Output files:      
#  originalrun.sav  
#    Target data (generated from the full model based on snmf_protocols.py)  

from neuron import h
import matplotlib
matplotlib.use('Agg')
import numpy
from pylab import *
import mytools
import pickle
import time
import sys
import random
import snmf_protocols

dists_apical = []
dists_basal = []
BACdt = 5.0

def initialize():
  global dists_apical, dists_basal, BACdt
  v0 = -80
  ca0 = 0.0001
  BACdt = 5.0
  fs = 8
  tstop = 15000.0
  icell = 0
  nsegs = 10
  nrecsperseg = 5

  if nrecsperseg == 1:
    xsperseg = [0.5]
  else:
    xsperseg = [1.0*x/(nrecsperseg+1) for x in range(1,nrecsperseg+1)]

  h("""
load_file("stdlib.hoc")
load_file("stdrun.hoc")
objref cvode
cvode = new CVode()
cvode.active(1)
cvode.atol(0.00005)
load_file("import3d.hoc")
objref L5PC              
load_file(\"models/L5PCbiophys3.hoc\")
load_file(\"models/L5PCtemplate.hoc\")
L5PC = new L5PCtemplate(\"morphologies/cell1.asc\")
access L5PC.soma
objref st1
objref vsoma, vdend, vdend2
vsoma = new Vector()
vdend = new Vector()
vdend2 = new Vector()
objref syn1, tvec, sl
tvec = new Vector()
sl = new List()
double siteVec[2]
sl = L5PC.locateSites("apic",620)       
maxdiam = 0                             
for(i=0;i<sl.count();i+=1){             
  dd1 = sl.o[i].x[1]                    
  dd = L5PC.apic[sl.o[i].x[0]].diam(dd1)
  if (dd > maxdiam) {                   
    j = i                               
    maxdiam = dd                        
  }                                     
}                                       
siteVec[0] = sl.o[j].x[0]               
siteVec[1] = sl.o[j].x[1]
L5PC.apic[siteVec[0]] cvode.record(&v(siteVec[1]),vdend,tvec)
L5PC.apic[siteVec[0]] syn1 = new epsp(siteVec[1])
sl = L5PC.locateSites("apic",800)       
maxdiam = 0                             
for(i=0;i<sl.count();i+=1){             
  dd1 = sl.o[i].x[1]                    
  dd = L5PC.apic[sl.o[i].x[0]].diam(dd1)
  if (dd > maxdiam) {                   
    j = i                               
    maxdiam = dd                        
  }                                     
}                                       
siteVec[0] = sl.o[j].x[0]               
siteVec[1] = sl.o[j].x[1]
L5PC.apic[siteVec[0]] cvode.record(&v(siteVec[1]),vdend2,tvec)
syn1.imax = 0
L5PC.soma cvode.record(&v(0.5),vsoma,tvec)
forall nseg = """+str(nsegs)+"""
tstop = """+str(tstop)+"""

objref vrecs_apical["""+str(nrecsperseg*109)+"""], vrecs_basal["""+str(nrecsperseg*84)+"""], carecs_apical["""+str(nrecsperseg*109)+"""]

cai0_ca_ion = """+str(ca0)+"""
v_init = """+str(v0)+"""
""")

  dists_apical = []
  dists_basal = []
  for i in range(0,109):
    h("access L5PC.apic["+str(i)+"]")
    for j in range(0,nrecsperseg):
      dists_apical.append(h.distance(xsperseg[j]))
      h("L5PC.apic["+str(i)+"] vrecs_apical["+str(i*nrecsperseg + j)+"] = new Vector()")
      h("L5PC.apic["+str(i)+"] carecs_apical["+str(i*nrecsperseg + j)+"] = new Vector()")
      h("L5PC.apic["+str(i)+"] cvode.record(&v("+str(xsperseg[j])+"),vrecs_apical["+str(i*nrecsperseg + j)+"],tvec)")
      h("L5PC.apic["+str(i)+"] cvode.record(&cai("+str(xsperseg[j])+"),carecs_apical["+str(i*nrecsperseg + j)+"],tvec)")
  for i in range(0,84):
    h("access L5PC.dend["+str(i)+"]")
    for j in range(0,nrecsperseg):
      dists_basal.append(h.distance(xsperseg[j]))
      h("L5PC.dend["+str(i)+"] vrecs_basal["+str(i*nrecsperseg + j)+"] = new Vector()")
      h("L5PC.dend["+str(i)+"] cvode.record(&v("+str(xsperseg[j])+"),vrecs_basal["+str(i*nrecsperseg + j)+"],tvec)")
  h("access L5PC.soma")
  maxLenApic = 1300
  maxLenBasal = 280


  STIMTYPES = snmf_protocols.get_stimulus_types()
  for istim_type in range(0,len(STIMTYPES)):
    ST = STIMTYPES[istim_type]
    if ST[3][0] == "fixed":
      segname = "L5PC."+ST[2]
      segx = ST[3][1]
    elif ST[3][0] == "distance":
      h("""
sl = L5PC.locateSites(\""""+ST[2]+"""\","""+str(ST[3][1])+""")
maxdiam = 0                             
for(i=0;i<sl.count();i+=1){             
  dd1 = sl.o[i].x[1]                    
  dd = L5PC."""+ST[2]+"""[sl.o[i].x[0]].diam(dd1)
  if (dd > maxdiam) {                   
    j = i                               
    maxdiam = dd                        
  }                                     
}                                       
siteVec[0] = sl.o[j].x[0]               
siteVec[1] = sl.o[j].x[1]
""")
      segname = "L5PC."+ST[2]+"["+str(int(h.siteVec[0]))+"]"
      segx = h.siteVec[1]
    else:
      print "Error: use \"fixed\" or \"distance\" for ST[3][0]!"
    print "istim_type = "+str(istim_type)
    print("""
objref """+ST[0]+"""
""" +segname + """ """ + ST[0] + """ = new """+ST[1]+"""("""+str(segx)+""")
""") 
    h("""
objref """+ST[0]+"""
""" +segname + """ """ + ST[0] + """ = new """+ST[1]+"""("""+str(segx)+""")
""")

  return [dists_apical, dists_basal]



def run(saveFig=False):

  STIMULI = snmf_protocols.get_stimuli()
  SETUP = snmf_protocols.get_setup()

  approx_to_real_seclist = {"soma": "somatic", "dend": "basal", "apic[0]": "apical", "apic[1]": "apical"}

  STIMTYPES = snmf_protocols.get_stimulus_types()
  DATASTORTYPES = snmf_protocols.get_data_storage_types()

  myReturn = [[]]*len(SETUP)
  for istep in range(len(SETUP)-1,-1,-1):
    vars_zero = snmf_protocols.get_variable_params(istep+1)
    for iparam in range(0,len(vars_zero)):
      if vars_zero[iparam][0]=='g' and vars_zero[iparam].find('bar') > -1:
        print vars_zero[iparam]
        underscoreind = vars_zero[iparam].rfind('_')
        section = vars_zero[iparam][underscoreind+1:len(vars_zero[iparam])]
        print "forsec L5PC."+approx_to_real_seclist[section]+" "+vars_zero[iparam][0:underscoreind]+" = 0" 
        h("forsec L5PC."+approx_to_real_seclist[section]+" "+vars_zero[iparam][0:underscoreind]+" = 0")

    myValsAllAll = []
    for istims in range(0,len(SETUP[istep])):
      stims = SETUP[istep][istims]
      myValsAll = []
      for iamp in range(0,len(stims[1])):
        for istim in range(0,len(stims[0])):
          stimulus = STIMULI[stims[0][istim]]
          st = STIMTYPES[stimulus[0]]          
          if type(stims[1][iamp]) is list:
            myamp = stims[1][iamp][istim]
          else:
            myamp = stims[1][iamp]
          print(st[0]+"."+st[4]+" = "+str(myamp))
          h(st[0]+"."+st[4]+" = "+str(myamp))
          for iprop in range(0,len(stimulus[1])):
            thisPropVal = stimulus[1][iprop][1]
            if stimulus[1][iprop][0] == "onset" or stimulus[1][iprop][0] == "del":
              thisPropVal = thisPropVal + istim*BACdt
            print st[0]+"."+stimulus[1][iprop][0]+" = "+str(thisPropVal)
            h(st[0]+"."+stimulus[1][iprop][0]+" = "+str(thisPropVal))
        h.init()
        timenow = time.time()
        h.run()
        print "Step "+str(istep)+", istims = "+str(istims)+", iamp = "+str(iamp)+" completed in "+str(time.time()-timenow)+" seconds"
        times=np.array(h.tvec)
        Vsoma=np.array(h.vsoma)
        Vdend=np.concatenate((np.array(h.vrecs_apical),np.array(h.vrecs_basal)))
        Cadend=np.array(h.carecs_apical)
        spikes = mytools.spike_times(times,Vsoma,-35,-45)
        if saveFig:
          close("all")
          f,axarr = subplots(2,2)
          axarr[0,0].plot(times,Vsoma)
          axarr[0,0].set_xlim([9990,11000])
          axarr[0,0].set_title("nSpikes="+str(len(spikes)))
          axarr[0,1].plot(dists_apical+[-x for x in dists_basal], [max([x[i] for i,t in enumerate(times) if t >= 9000]) for x in Vdend], 'bx')
          axarr[1,0].plot(dists_apical, [max([x[i] for i,t in enumerate(times) if t >= 9000]) for x in Cadend], 'b.')
          f.savefig("original_step"+str(istep)+"_istims"+str(istims)+"_iamp"+str(iamp)+".eps")
        #times = range(0,15000); Vsoma = range(0,15000); Vdend = [range(0,15000), range(0,15000)]; Cadend = Vdend; spikes = []
        print "istep="+str(istep)+", istims="+str(istims)+", iamp="+str(iamp)+": nSpikes="+str(len(spikes))
        irecs = stims[2]
        if type(irecs) is not list:
          irecs = [irecs]
        myVals = []
        for iirec in range(0,len(irecs)):
          irec = irecs[iirec]
          if irec == 0:
            myData = [Vsoma]
          elif irec == 1:
            myData = Vdend
          elif irec == 2:
            myData = Cadend
          else:
            print "Unknown recording type: "+str(irec)
            continue
          print "istep="+str(istep)+": "
          print DATASTORTYPES[stims[3]]
          try:
            print "min(myData)="+str(min(myData))+", max(myData)="+str(max(myData))
          except:
            try:
              print "min(myData[0])="+str(min(myData[0]))+", max(myData[0])="+str(max(myData[0]))
            except:
              print myData
          if DATASTORTYPES[stims[3]][0] == "fixed":
            print str(next(i for i,t in enumerate(times) if t >= DATASTORTYPES[stims[3]][1]))
            myVals.append([x[next(i for i,t in enumerate(times) if t >= DATASTORTYPES[stims[3]][1])-1] for x in myData])
          elif DATASTORTYPES[stims[3]][0] == "max":
            myVals.append([max(x[next(i for i,t in enumerate(times) if t >= DATASTORTYPES[stims[3]][1][0]):next(i for i,t in enumerate(times) if t >= DATASTORTYPES[stims[3]][1][1])]) for x in myData])
          elif DATASTORTYPES[stims[3]][0] == "trace" or DATASTORTYPES[stims[3]][0] == "highrestrace":
            myVals.append([mytools.interpolate(times,x,DATASTORTYPES[stims[3]][1]) for x in myData])
          elif DATASTORTYPES[stims[3]][0] == "highrestraceandspikes":
            myVals.append([[mytools.interpolate(times,x,DATASTORTYPES[stims[3]][1]) for x in myData],spikes])
          elif DATASTORTYPES[stims[3]][0] == "nspikes":
            myVals.append(sum([1 for x in spikes if x >= DATASTORTYPES[stims[3]][1][0] and x < DATASTORTYPES[stims[3]][1][1]]))
          else:
            print "Unknown data storage type: "+DATASTORTYPES[stims[3]][0]
            continue
        myValsAll.append(myVals[:])
      myValsAllAll.append(myValsAll[:])
          
      for istim in range(0,len(stims[0])):
        stimulus = STIMULI[stims[0][istim]]
        st = STIMTYPES[stimulus[0]]          
        print(st[0]+"."+st[4]+" = 0")
        h(st[0]+"."+st[4]+" = 0")

    myReturn[istep] = myValsAllAll[:]
  picklelist = [myReturn, dists_apical, dists_basal]
  file = open('originalrun.sav', 'w')
  pickle.dump(picklelist,file)
  file.close()
  #print "run():"
  #mytools.printlistlen(myReturn)
  return myReturn


#calculate the error from mean experimental values (measured in terms of SDs, except #APs, which are absolute)
def howfarfrommeans():
  A = run()
  errors = [mean([abs(A[0][0]-9.0)/0.88,      abs(A[0][1]-14.5)/0.56,     abs(A[0][2]-22.5)/2.22]),
            mean([abs(A[1][0]-0.0036)/0.0091, abs(A[1][1]-0.0023)/0.0056, abs(A[1][2]-0.0046)/0.0026]),
            mean([abs(A[2][0]-0.1204)/0.0321, abs(A[2][1]-0.1086)/0.0368, abs(A[2][2]-0.0954)/0.0140]),
            mean([abs(A[3][0]-57.75)/33.48,   abs(A[3][1]-6.625)/8.65,    abs(A[3][2]-5.38)/0.83]),
            mean([abs(A[4][0]-43.25)/7.32,    abs(A[4][1]-19.13)/7.31,    abs(A[4][2]-7.25)/1.0]),
            mean([abs(A[5][0]-26.23)/4.97,    abs(A[5][1]-16.52)/6.11,    abs(A[5][2]-16.44)/6.93]),
            mean([abs(A[6][0]+58.04)/4.58,    abs(A[6][1]+60.51)/4.67,    abs(A[6][2]+59.99)/3.92]),
            mean([abs(A[7][0]-0.238)/0.03,    abs(A[7][1]-0.279)/0.027,   abs(A[7][2]-0.213)/0.037]),
            mean([abs(A[8][0]-1.31)/0.17,     abs(A[8][1]-1.38)/0.28,     abs(A[8][2]-1.86)/0.41]),
            abs(A[9]-6.73)/2.54,
            abs(A[10]-37.43)/1.27,
            abs(A[11]-3.0),
            abs(A[12]-9.9)/0.85,
            abs(A[13]+65.0)/4.0,
            abs(A[14]-25.0)/5.0,
            abs(A[15]-2.0)/0.5,
            abs(A[16]-1.0)]
  return errors
            
  
initialize()
run(True)

Loading data, please wait...