from neuron import h, gui
from math import * #exp,log, sqrt
from axonMorph import *
h.load_file('import3d.hoc')
####################################################################
# Setting up params
####################################################################
# h-current
Vrest = -90.220509483
h.v_init = -70.0432010302
h.celsius = 34.0 # same as experiment
h.tstop = 3000
curr = 0.300 # stimulus amplitude
# passive properties (optimization proxies)
cap = 0.700445928608
rall = 137.494564931
rm = 38457.4393085
spinecapfactor = 1.48057846279
#Na, K reversal potentials calculated from BenS internal and external solutions via Nernst equation
p_ek = -104
p_ena = 42
# h-current
h_erev = -37.0
h_gbar = 6.6129403774e-05 # mho/cm^2
h_gbar_tuft = 0.00565 # mho/cm^2 (based on Harnett 2015 J Neurosci)
# d-current
kdmc_gbar = 0.00110907315064
kdmc_gbar_axonm = 20
# spiking currents
nax_gbar = 0.0153130368342
nax_gbar_axonm = 5
kdr_gbar = 0.0084715576279
kdr_gbar_axonm = 5
# A few kinetic params changed vis-a-vis kdr.mod defaults:
# These param values match the default values in kdr.mod (except where indicated).
kdr_vhalfn = 13
kdr_a0n = 0.0075 # def 0.02
kdr_nmax = 20 # def 2
kap_gbar = 0.0614003391814
kap_gbar_axonm = 5
# These param values match the default values in kap.mod (except where indicated).
kap_vhalfn = 35 # def 11
kap_vhalfl = -56
kap_nmin = 0.4 # def 0.1
kap_lmin = 5 # def 2
kap_tq = -45 # def -40
# new ion channel parameters
cal_gcalbar = 5.73945708921e-06
can_gcanbar = 4.74427101753e-06
calginc = 1.0
kBK_gpeak = 7.25128017201e-05 # original value of 268e-4 too high for this model
kBK_caVhminShift = 46.9679440782 # shift upwards to get lower effect on subthreshold
cadad_depth = 0.102468419281 # original 1; reduced for tighter coupling of ica and cai
cadad_taur = 16.0181691392
h.erev_ih = h_erev # global
# K delayed rectifier current
h.a0n_kdr = kdr_a0n # global
h.nmax_kdr = kdr_nmax # global
# K-A current
h.nmin_kap = kap_nmin # global
h.lmin_kap = kap_lmin # global
###############################################################################
# Morph Cell
###############################################################################
class Cell:
def __init__(self,morph):
morph = str(morph)
self.axon = []
cell = h.Import3d_Neurolucida3()
cell.input(morph)
i3d = h.Import3d_GUI(cell, 0)
i3d.instantiate(self)
self.add_axon()
self.init_once()
def add_axon(self):
self.axon.append(h.Section(name="axon[0]"))
self.all.append(self.axon[0])
self.axon[0].connect(self.soma[0], 0.0, 0.0)
# clears 3d points
h.pt3dclear()
# define a logical connection point relative to the first 3-d point
h.pt3dstyle(axonPts[0][0], axonPts[0][1], axonPts[0][2], axonPts[0][3], sec=self.axon[0])
# add axon points after first logical connection point
for x, y, z, d in axonPts[1:]:
h.pt3dadd(x, y, z, d, sec=self.axon[0])
def init_once(self):
for sec in self.all: #forall within an object just accesses the sections belonging to the object
sec.insert('pas') # passive
sec.insert('savedist') # mechanism to keep track of distance from soma even after multisplit
sec.insert('ih') # h-current in Ih_kole.mod
sec.insert('nax') # Na current
sec.insert('kdr') # K delayed rectifier current
sec.insert('kap') # K-A current
sec.insert('k_ion')
sec.insert('na_ion')
if not (h.ismembrane('ih', sec= sec)):
print(sec)
self.setallprop()
self.addapicchan()
self.apicchanprop()
self.addbasalchan()
self.basalchanprop()
# spiny:
# for i=2, 102 apic[i] spiny.append()
# for i=0, 68 dend[i] spiny.append()
# Skips first two apicals
for i in range(len(self.apic)):
if (i != 0) and (i != 1):
self.apic[i].cm = spinecapfactor * self.apic[i].cm # spinecapfactor * cm
self.apic[i].g_pas = spinecapfactor / rm # spinecapfactor * (1.0/rm)
for sec in self.dend:
sec.cm = spinecapfactor * sec.cm
sec.g_pas = spinecapfactor / rm
self.optimize_nseg()
self.setgbarnaxd()
self.setgbarkapd() # kap in apic,basal dends
self.setgbarkdrd() # kdr in apic,basal dends
self.sethgbar() # distributes HCN conductance
for sec in self.soma:
sec.insert('kdmc')
sec.insert('ca_ion')
sec.insert('cadad')
sec.insert('cal')
sec.insert('can')
sec.insert('kBK')
self.setsomag()
# axon has I_KD, and more I_Na, I_KA, I_KDR, and no I_h
# K-D current in soma and axon only
for sec in self.axon:
sec.insert('kdmc')
self.setaxong()
def addapicchan(self):
for sec in self.apic:
sec.insert('ca_ion')
sec.insert('cadad') # cadad.mod
sec.insert('cal') # cal.mod
sec.insert('can') # can.mod
sec.insert('kBK') #kBK.mod
def apicchanprop(self):
for sec in self.apic:
sec.gcalbar_cal = cal_gcalbar
sec.gcanbar_can = can_gcanbar
sec.gpeak_kBK = kBK_gpeak
sec.caVhmin_kBK = -46.08 + kBK_caVhminShift
sec.depth_cadad = cadad_depth
sec.taur_cadad = cadad_taur
def addbasalchan(self):
# basal == dend
for sec in self.dend:
sec.insert('ca_ion')
sec.insert('cadad') # cadad.mod
sec.insert('cal') # cal.mod
sec.insert('can') # can.mod
sec.insert('kBK') #kBK.mod
def basalchanprop(self):
# basal == dend
for sec in self.dend:
sec.gcalbar_cal = cal_gcalbar
sec.gcanbar_can = can_gcanbar
sec.gpeak_kBK = kBK_gpeak
sec.caVhmin_kBK = -46.08 + kBK_caVhminShift
sec.depth_cadad = cadad_depth
sec.taur_cadad = cadad_taur
def geom_nseg(self):
# local freq, d_lambda, before, after, tmp
# these are reasonable values for most models
freq = 100 # Hz, frequency at which AC length constant will be computed
d_lambda = 0.1
before = 0
after = 0
for sec in self.all:
before += sec.nseg
#soma area(0.5) # make sure diam reflects 3d points
for sec in self.all:
# creates the number of segments per section
# lambda_f takes in the current section
sec.nseg = int((sec.L/(d_lambda*self.lambda_f(sec))+0.9)/2)*2 + 1
for sec in self.all:
after += sec.nseg
print("geom_nseg: changed from ", before, " to ", after, " total segments")
def lambda_f(self, section):
# these are reasonable values for most models
freq = 100 # Hz, frequency at which AC length constant will be computed
d_lambda = 0.1
# The lowest number of n3d() is 2
if (section.n3d() < 2):
return 1e5*sqrt(section.diam/(4*pi*freq*section.Ra*section.cm))
# above was too inaccurate with large variation in 3d diameter
# so now we use all 3-d points to get a better approximate lambda
x1 = section.arc3d(0)
d1 = section.diam3d(0)
self.lam = 0
#print section, " n3d:", section.n3d(), " diam3d:", section.diam3d(0)
for i in range(section.n3d()): #h.n3d()-1
x2 = section.arc3d(i)
d2 = section.diam3d(i)
self.lam += (x2 - x1)/sqrt(d1 + d2)
x1 = x2
d1 = d2
# length of the section in units of lambda
self.lam *= sqrt(2) * 1e-5*sqrt(4*pi*freq*section.Ra*section.cm)
return section.L/self.lam
def optimize_nseg(self):
# Ra, cm, spinecapfactor
# use worst case Ra, cm values
for sec in self.all:
sec.Ra = rall
sec.cm = cap
# spiny:
# for i=2, 102 apic[i] spiny.append()
# for i=0, 68 dend[i] spiny.append()
# Skips the first two apicals
for i in range(len(self.apic)):
if (i != 0) and (i != 1):
self.apic[i].cm = spinecapfactor * self.apic[i].cm # spinecapfactor * cm
self.apic[i].g_pas = spinecapfactor / rm # spinecapfactor * (1.0/rm)
for sec in self.dend:
sec.cm = spinecapfactor * sec.cm
sec.g_pas = spinecapfactor / rm
# optimize # of segments per section (do this afer setting Ra, cm)
self.geom_nseg()
# need to reassign distances after changing nseg
self.recalculate_x_dist()
def setallprop(self):
for sec in self.all:
# passive
sec.g_pas = 1 / rm
sec.Ra = rall
sec.cm = cap
sec.e_pas = Vrest
# Na current
sec.gbar_nax = nax_gbar
# K delayed rectifier current
sec.gbar_kdr = kdr_gbar
sec.vhalfn_kdr = kdr_vhalfn
# K-A current
sec.gbar_kap = kap_gbar
sec.vhalfn_kap = kap_vhalfn
sec.vhalfl_kap = kap_vhalfl
sec.tq_kap = kap_tq
# reversal potentials
sec.ena = p_ena
sec.ek = p_ek
def setaxong(self):
# axon has I_KD, and more I_Na, I_KA, I_KDR, and no I_h
for sec in self.axon:
# axon has no Ih
sec.gbar_ih = 0
# increase the I_Na, I_Ka, I_KD and I_KDR
sec.gbar_nax = nax_gbar_axonm * nax_gbar
sec.gbar_kap = kap_gbar_axonm * kap_gbar
sec.gbar_kdmc = kdmc_gbar_axonm * kdmc_gbar
sec.gbar_kdr = kdr_gbar_axonm * kdr_gbar
def setgbarnaxd(self):
# alldend:
# for i=0, 68 dend[i] alldend.append()
# for i=0, 102 apic[i] alldend.append()
for sec in self.dend:
sec.gbar_nax = nax_gbar
for sec in self.apic:
sec.gbar_nax = nax_gbar
def setgbarkapd(self):
# alldend:
# for i=0, 68 dend[i] alldend.append()
# for i=0, 102 apic[i] alldend.append()
for sec in self.dend:
sec.gbar_kap = kap_gbar
for sec in self.apic:
sec.gbar_kap = kap_gbar
def setgbarkdrd(self):
# alldend:
# for i=0, 68 dend[i] alldend.append()
# for i=0, 102 apic[i] alldend.append()
for sec in self.dend:
sec.gbar_kdr = kdr_gbar
for sec in self.apic:
sec.gbar_kdr = kdr_gbar
def setsomag(self):
for sec in self.soma:
sec.gbar_kdmc = kdmc_gbar
sec.gcalbar_cal = cal_gcalbar
sec.gcanbar_can = can_gcanbar
sec.gpeak_kBK = kBK_gpeak
sec.caVhmin_kBK = -46.08 + kBK_caVhminShift
sec.depth_cadad = cadad_depth
sec.taur_cadad = cadad_taur
def recalculate_x_dist(self):
# set the center of the soma as origin
h.distance(0,0.5,sec=self.soma[0])
for sec in self.all:
for seg in sec.allseg():
seg.x_savedist = h.distance(seg.x,sec=sec)
# stays nested to compare soma to apicals
# list of apicals that are Upper Trunk
self.nexusdist = nexusdist = 1e9 # calculate distance to tuft here
apicUpperTrunk = [23, 24, 25, 32, 33, 46, 51, 52, 53, 78]
for apicIndx in apicUpperTrunk:
# excludes any apicals not in the list
for seg in self.apic[apicIndx].allseg():
nd = h.distance(seg.x,sec=self.apic[apicIndx])
if (nd < self.nexusdist):
self.nexusdist = nd
def reconfig (self):
self.setallprop() # set initial properties, including g_pas,e_pas,cm,Ra,etc.
self.optimize_nseg() # set nseg based on passive properties
self.apicchanprop() # initial apic dend properties
self.basalchanprop() # initial basal dend properties
self.setgbarnaxd() # nax in apic,basal dends
self.setgbarkapd() # kap in apic,basal dends
self.setgbarkdrd() # kdr in apic,basal dends
self.sethgbar() # h in all locations
self.setsomag() # soma-specific conductance values
self.setaxong() # axon-specific conductance values
def sethgbar(self):
h.distance(0,0.5,sec=self.soma[0]) # middle of soma is origin for distance
self.h_gbar_tuftm = h_gbar_tuftm = h_gbar_tuft / h_gbar
self.h_lambda = h_lambda = self.nexusdist / log(h_gbar_tuftm)
for sec in self.soma:
sec.gbar_ih = h_gbar
# basal:
# basal == all dends
for sec in self.dend:
sec.gbar_ih = h_gbar
# apical_oblique:
# for i=79, 102 apic[i]
for i in range(79,103):
self.apic[i].gbar_ih = h_gbar
#apical_maintrunk:
# for i=0, 22 apic[i]
for i in range (0,23):
for seg in self.apic[i]:
xd = h.distance(seg.x,sec=sec)
if (xd <= self.nexusdist):
seg.gbar_ih = h_gbar * exp(seg.x_savedist/self.h_lambda)
else:
seg.gbar_ih = h_gbar_tuft
if (seg.gbar_ih < 0):
seg.gbar_ih = 0
# apical upper trunk list
apicUpperTrunk = [23,24,25,32,33,46,51,52,53,78]
for i in apicUpperTrunk:
self.apic[i].gbar_ih = h_gbar * h_gbar_tuftm
# apical tuft:
# i=26, 31 apic[i]
# i=34, 45 apic[i]
# i=47, 50 apic[i]
# i=54, 77 apic[i]
for i in range(26,32):
self.apic[i].gbar_ih = h_gbar * h_gbar_tuftm
for i in range(34,46):
self.apic[i].gbar_ih = h_gbar * h_gbar_tuftm
for i in range(47,50):
self.apic[i].gbar_ih = h_gbar * h_gbar_tuftm
for i in range(54,77):
self.apic[i].gbar_ih = h_gbar * h_gbar_tuftm
# Initializing cells
#c1 = Cell('BS0284')
#c2 = Cell('BS0409')
# Showing 3D Plot of Cells
#shape_window = h.PlotShape()