# ---------------------------------------------------------- # Library of cell classes and functions # # Naoki Hiratani, N.Hiratani@gmail.com # # Based on the model by Dr. Tiago Branco # ---------------------------------------------------------- import numpy as np import neuron from neuron import h from neuron import load_mechanisms from neuron import gui from math import * import random as pyrnd import matplotlib.pyplot as plt import scipy.special as scisp h('objref nil') # ---------------------------------------------------------- # MODELS class L23(object): def __init__(self): h('xopen("./L23.hoc")') props(self) self._geom() self._topol() self._biophys() self._dist_from_soma() def _geom(self): self.axon = h.Section() self.axon.L = 300 self.axon.diam = 1 def _topol(self): self.soma = h.soma self.axon.connect(self.soma,1,0) #self.axon = h.axon self.dends = [] #list of dendritic section for sec in h.allsec(): self.dends.append(sec) self.dends.remove(self.soma) # Remove soma from the list self.dends.remove(self.axon) # and the Axon dnames = ['dend1_1','dend1_11','dend1_111','dend1_1111','dend1_1112','dend1_112','dend1_1121','dend1_1122','dend1_12','dend1_121','dend1_1211','dend1_1212','dend1_122','dend1_1221','dend1_1222','dend1_12221','dend1_12222','dend2_1','dend2_11','dend2_111','dend2_1111','dend2_1112','dend2_112','dend2_1121','dend2_1122','dend2_12','dend2_121','dend2_1211','dend2_12111','dend2_121111','dend2_121112','dend2_12112','dend2_121121','dend2_1211211','dend2_1211212','dend2_12112121','dend2_12112122','dend2_121122','dend2_1212','dend2_12121','dend2_121211','dend2_121212','dend2_12122','dend2_122','dend2_1221','dend2_12211','dend2_12212','dend2_1222','dend2_12221','dend2_12222','dend3_1','dend3_11','dend3_111','dend3_1111','dend3_1112','dend3_11121','dend3_11122','dend3_112','dend3_1121','dend3_1122','dend3_12','dend3_121','dend3_1211','dend3_1212','dend3_12121','dend3_12122','dend3_121221','dend3_121222','dend3_1212221','dend3_1212222','dend3_12122221','dend3_12122222','dend3_121222221','dend3_1212222211','dend3_12122222111','dend3_12122222112','dend3_1212222212','dend3_121222222','dend3_1212222221','dend3_1212222222','dend3_12122222221','dend3_12122222222','dend3_122','dend3_1221','dend3_12211','dend3_12212','dend3_1222','dend4_1','dend4_11','dend4_111','dend4_1111','dend4_1112','dend4_11121','dend4_11122','dend4_112','dend4_1121','dend4_1122','dend4_12','dend4_121','dend4_122','dend4_1221','dend4_1222','dend4_12221','dend4_12222'] def _biophys(self): for sec in h.allsec(): sec.cm = self.CM sec.insert('pas') sec.e_pas = self.E_PAS sec.g_pas = 1.0/self.RM sec.Ra = self.RA sec.nseg = int(ceil(sec.L/5.0)) #5um segmentation def _dist_from_soma(self): self.dists = [] #distance from the soma to dendrites h('soma distance(0.0,0.5)') lidx = 0 for sec in h.allsec(): if lidx < len(self.dends): h('''proc calc_distance(){ disttmp = distance(1.0,0.0)}''') h.calc_distance() self.dists.append(h.disttmp) lidx += 1 self.branchdists = [] #dendritic distance between two branches l1idx = 0 for sec1 in h.allsec(): if l1idx < len(self.dends): self.branchdists.append([]) h('distance()') l2idx = 0 for sec2 in h.allsec(): if l2idx < len(self.dends): h('''proc calc_distance(){ disttmp = distance(1.0,0.0)}''') h.calc_distance() self.branchdists[l1idx].append( h.disttmp ) l2idx += 1 l1idx += 1 # ---------------------------------------------------------- # INSTRUMENTATION FUNCTIONS def props(model): # Passive properties model.CM = 1.0 model.RM = 7000.0 model.RA = 100 model.E_PAS = -75 model.CELSIUS = 35 # Active properties model.Ek = -90 model.Ena = 60 model.Eca = 140 model.gna_axon = 1000 model.gkv_axon = 100 model.gna_soma = 1000 model.gkv_soma = 100 model.gkm_soma = 2.2 model.gkca_soma = 3 model.gca_soma = 0.5 model.git_soma = 0.0003 model.gna_dend = 27#80 model.gna_dend_hotSpot = 600 model.gkv_dend = 1#3 model.gkm_dend = 0.3#1 model.gkca_dend = 3 model.gca_dend = 0.5 model.git_dend = 0.00015 model.gh_dend = 0 def init_active(model, axon=False, soma=False, dend=True, dendNa=False, dendCa=False): if axon: model.axon.insert('na'); model.axon.gbar_na = model.gna_axon model.axon.insert('kv'); model.axon.gbar_kv = model.gkv_axon model.axon.ena = model.Ena model.axon.ek = model.Ek if soma: model.soma.insert('na'); model.soma.gbar_na = model.gna_soma model.soma.insert('kv'); model.soma.gbar_kv = model.gkv_soma model.soma.insert('km'); model.soma.gbar_km = model.gkm_soma model.soma.insert('kca'); model.soma.gbar_kca = model.gkca_soma model.soma.insert('ca'); model.soma.gbar_ca = model.gca_soma model.soma.insert('it'); model.soma.gbar_it = model.git_soma #model.soma.insert('cad'); model.soma.ena = model.Ena model.soma.ek = model.Ek model.soma.eca = model.Eca if dend: for d in model.dends: d.insert('na'); d.gbar_na = model.gna_dend*dendNa d.insert('kv'); d.gbar_kv = model.gkv_dend d.insert('km'); d.gbar_km = model.gkm_dend d.insert('kca'); d.gbar_kca = model.gkca_dend d.insert('ca'); d.gbar_ca = model.gca_dend*dendCa d.insert('it'); d.gbar_it = model.git_dend*dendCa #d.insert('cad') d.ena = model.Ena d.ek = model.Ek d.eca = model.Eca def init_params(model, Kin, gmax, gI, uk_min, sd_bias, release_prob): #simulation control model.trials = 1001# #trial number #synapses model.Min = 200 #number of input neurons model.Kin = Kin #redanduncy in synaptic connections model.Nin = model.Kin*model.Min #number of synaptic inputs model.gmax = gmax#0.0015 #standard conductance [muS] model.sd_bias = sd_bias #bias in synaptic distribution model.uk_min = uk_min #lower boundary condition of uks model.uk_max = 1.0 #upper boundary condition of uks #rewiring parameters model.ukthreshold = model.uk_min #spine cutoff threshold model.ukinit = 1.0/float(Kin) #model.uk_min #size of a newly created spine model.taumr = 10.0 #the time constant for mean firing rate model.rewiring_freq = 0.2 #relative frequency of rewiring model.pre_rates_th = 0.05 #threshold for the presynaptic firing rate in supervised trials #input structure model.spn_rate = 1.5*pi #spontaneous firing rate model.epsilon = 0.01 #the minimum relative rate model.min_rate = model.epsilon*model.spn_rate model.release_prob = release_prob# presynaptic release probability model.p_theta = 0.0 #preferred orientation model.n_theta = model.p_theta + pi/2.0 #non-preferred orientation model.rfdist_zero = 1.0 # standard distance between somatic and synaptic receptive fields model.rfdist_min = 0.01 # minimum distance between somatic and synaptic receptive fields model.rfdist_max = 3.0 # maximum distance between somatic and synaptic receptive fields model.kappa_zero = 2.0 # orientation selectivity model.kappa_phi = 4.0 # association field selectivity model.krfdist_const = exp(model.kappa_phi)/100.0 #constant for rfdist in calculation of k_i (to avoid numerical instability) #inhibitory synapses model.inhN = 200 #Number of inhibitory synapses model.gI = gI #inhibitory conductance model.Inh_thr = -90.0 #mV #stimulation protocol model.sstart = 20.0 #starting timing of the stimulation (from the initiation of the simulation) model.sduration = 20.0 #stimulation duration[ms] model.snoise = 0.0 #noise in spike train def plot_morphology(model,locs): h('access soma') h('objref sh') h('sh = new PlotShape(0)') h('sh.size(-300,300,-299.522,299.522)') h('''proc plot_synapse(){ xdtmp = x3d($1) ydtmp = y3d($1) ctmp = 7 sh.mark(xdtmp,ydtmp,"o",6,ctmp,4) }''') lidx = 0 for sec in h.allsec(): if lidx < len(model.dends): for i in range(len(locs)): loc = locs[i] if loc[0] == lidx: h('n3dtmp = n3d()') ndtmp = int(floor(h.n3dtmp*loc[1])) h.plot_synapse(ndtmp,model.preidx[i]) lidx += 1 h('sh.view(-300, -400.522, 600, 900.043, 265, 450, 200.64, 400.32)') def add_Estims(model,locs,sstart=20.0,sinterval=20,snumber=1,snoise=0): model.Estimlist = [] lidx = 0 for loc in locs: Estim = h.NetStim() Estim.interval = sinterval Estim.number = snumber Estim.start = sstart Estim.noise = snoise model.Estimlist.append(Estim) lidx += 1 def add_Istims(model,inh_locs,sstart=20.0,sinterval=20,snumber=1,snoise=0): model.Istimlist = [] lidx = 0 for loc in inh_locs: Istim = h.NetStim() Istim.interval = sinterval Istim.number = snumber Istim.start = sstart Istim.noise = snoise model.Istimlist.append(Istim) lidx += 1 #select dendritic locations for calculation of the unit EPSPs def calc_locs_kappa(model): dL = 5.0 locs = [] lidx = 0 nidx = 0 for sec in model.dends: nsec = int(ceil(sec.L/dL)) for sidx in range( nsec ): locs.append([]) locs[nidx].append( lidx ) locs[nidx].append( (sidx+0.5)/float(nsec) ) if locs[nidx][1] < 0.0 or locs[nidx][1] > 1.0: print lidx,nidx,locs[nidx][0],locs[nidx][1] nidx += 1 lidx += 1 #print '#kappa_locs',nidx return locs #Uniform selection of dendritic locations of the N synaptic contacts def calc_locs_uniform(model): locs = [] Ltot = 0.0 for sec in model.dends: Ltot += sec.L dL = Ltot/(model.Nin+1.0) Ltmp = 0.0 nidx = 0 lidx = 0 for sec in model.dends: Ltmp += sec.L while nidx*dL < Ltmp and Ltmp <= (nidx+1)*Ltmp: if nidx > 0 and nidx <= model.Nin: locs.append([]) locs[nidx-1].append(lidx) locs[nidx-1].append( 1.0 - (Ltmp-nidx*dL)/(sec.L) ) nidx += 1 lidx += 1 #print Ltot, dL, len(locs) return locs #Random selection of dendritic locations of the N synaptic contacts from a Beta distribution characterized by (bsd, 2-bsd) def calc_locs_random(model): locs = [] Lcums = [] Lcums.append(0.0) for sec in model.dends: Lcums.append(Lcums[-1] + sec.L) Ltot = Lcums[-1] dmax = 500.0#max(model.dists) bsda = model.sd_bias; bsdb = 2.0 - model.sd_bias #parameters for Beta distribution Zbeta = 10.0*scisp.beta(bsda,bsdb) #print dmax, bsda, bsdb, Zbeta nidx = 0 while(nidx < model.Nin): Ltmp = Ltot*pyrnd.random() for j in range(len(model.dends)): if Lcums[j] <= Ltmp and Ltmp < Lcums[j+1]: rdtmp = (model.dists[j] + (Ltmp-Lcums[j]))/dmax if pyrnd.random() < pow(rdtmp,bsda-1.0)*pow(1.0-rdtmp,bsdb-1.0)/Zbeta: locs.append([]) locs[nidx].append(j) locs[nidx].append( (Ltmp-Lcums[j])/(Lcums[j+1]-Lcums[j]) ) nidx += 1 #print Ltot, len(locs) return locs #Uniform selection of dendritic locations of the inhN synaptic contacts def calc_inh_locs_uniform(model): inh_locs = [] Ltot = 0.0 for sec in model.dends: Ltot += sec.L dL = Ltot/(model.inhN+1.0) Ltmp = 0.0 nidx = 0 lidx = 0 for sec in model.dends: Ltmp += sec.L while nidx*dL < Ltmp and Ltmp <= (nidx+1)*Ltmp: if nidx > 0 and nidx <= model.inhN: inh_locs.append([]) inh_locs[nidx-1].append(lidx) inh_locs[nidx-1].append( 1.0 - (Ltmp-nidx*dL)/(sec.L) ) nidx += 1 lidx += 1 print Ltot, dL, len(inh_locs) return inh_locs #restricted resampling: the posisiton of the new synapse is restricted to the set of dendritic branches where original connections were made def restricted_resampling(model, locs, i): loctmp = [] Lcums = [] Lcums.append(0.0) for j in range(model.Nin): if model.preidx[j] == model.preidx[i]: Lcums.append(Lcums[-1] + model.dends[locs[j][0]].L) Ltot = Lcums[-1] Ltmp = Ltot*pyrnd.random() lidx = 0 for j in range(model.Nin): if model.preidx[j] == model.preidx[i]: if Lcums[lidx] <= Ltmp and Ltmp < Lcums[lidx+1]: loctmp.append( locs[j][0]) loctmp.append( (Ltmp-Lcums[lidx])/(Lcums[lidx+1]-Lcums[lidx]) ) lidx += 1 return loctmp def calc_dist_bt_syns(model,locs):#calculate the distance between synapses disttmps = [] for i1 in range(model.Nin): jidx = model.preidx[i1] d1 = locs[i1][0] for i2 in range(i1+1,model.Nin): if model.preidx[i2] == jidx: d2 = locs[i2][0] disttmp = model.branchdists[d1][d2] + model.dends[d1].L*locs[i1][1] + model.dends[d2].L*locs[i2][1] disttmps.append( disttmp ) return disttmps #presynaptic index allocation def allocate_syns(model): rndidx = range(model.Nin) pyrnd.shuffle(rndidx) model.preidx = [] for i in range(model.Nin): model.preidx.append(0) for iidx in range(model.Nin): model.preidx[rndidx[iidx]] = iidx/model.Kin #add AMPA synapses def add_AMPAsyns(model, locs=[[0, 0.5]], gmax=0.5, tau1=0.5, tau2=2.5): model.AMPAlist = [] model.ncAMPAlist = [] for lidx in range(len(locs)): loc = locs[lidx] AMPA = h.Exp2Syn(float(loc[1]), sec=model.dends[int(loc[0])]) AMPA.tau1 = tau1 AMPA.tau2 = tau2 #NC = h.NetCon(h.nil, AMPA, 0, 0, gmax) NC = h.NetCon(model.Estimlist[lidx], AMPA, 0, 0, gmax) model.AMPAlist.append(AMPA) model.ncAMPAlist.append(NC) def add_GABAsyns(model, inh_locs=[[0, 0.5]], gI=0.5, tau1=0.5, tau2=2.5): model.GABAlist = [] model.ncGABAlist = [] for lidx in range(len(inh_locs)): inh_loc = inh_locs[lidx] GABA = h.Exp2Syn(float(inh_loc[1]), sec=model.dends[int(inh_loc[0])]) GABA.tau1 = tau1 GABA.tau2 = tau2 GABA.e = model.Inh_thr NC = h.NetCon(model.Istimlist[lidx], GABA, 0, 0, model.gI) model.GABAlist.append(GABA) model.ncGABAlist.append(NC) def rewire_synapse(model,loc,lidx,tau1=0.5,tau2=2.5): AMPAtmp = h.Exp2Syn(float(loc[1]), sec=model.dends[int(loc[0])]) AMPAtmp.tau1 = tau1 AMPAtmp.tau2 = tau2 NCtmp = h.NetCon(model.Estimlist[lidx], AMPAtmp, 0, 0, model.gmax) model.AMPAlist[lidx] = AMPAtmp model.ncAMPAlist[lidx] = NCtmp def generate_pre_tuning(model): model.phis = [] #polar direction of RF model.thetas = [] #preferred orientation model.rfdists = [] #Receptive field distance rfdist_zero = model.rfdist_zero rfdist_min = model.rfdist_min rfdist_max = model.rfdist_max for j in range(model.Min): model.phis.append( 2.0*pi*np.random.random() ) model.thetas.append( pi*np.random.random() ) model.rfdists.append( rfdist_min + (rfdist_max-rfdist_min)*np.random.random() ) def generate_log_rates(model, thetao): rfdist_zero = model.rfdist_zero krfdist_const = model.krfdist_const kappa_zero = model.kappa_zero kappa_phi = model.kappa_phi min_rate = model.min_rate spn_rate = model.spn_rate rates = []; log_rates = [] for j in range(model.Min): phi = model.phis[j]; rfdist = model.rfdists[j]; theta = model.thetas[j] kappa_r = rfdist_zero*exp(kappa_phi*cos(2.0*(phi-thetao)))/(rfdist + krfdist_const) ktmp = sqrt( kappa_zero*kappa_zero + kappa_r*kappa_r + 2*kappa_zero*kappa_r*cos( 2.0*(theta-thetao) ) ) rates.append( spn_rate*exp(-rfdist/rfdist_zero)*scisp.i0(ktmp)/(2*pi*scisp.i0(kappa_zero)*scisp.i0(kappa_r)) ) if rates[-1] > min_rate: log_rates.append( log(rates[-1]/min_rate) ) else: log_rates.append( 0.0 ) return rates, log_rates def generate_input_rates(model): # spontaneous firing rates model.s_rates = [] for j in range(model.Min): model.s_rates.append( model.spn_rate ) # firing rate for the preferred orientation p_rates, log_p_rates = generate_log_rates(model, model.p_theta) model.p_rates = p_rates model.log_p_rates = log_p_rates # firing rate for the non-preferred orientation n_rates, log_n_rates = generate_log_rates(model, model.n_theta) model.n_rates = n_rates model.log_n_rates = log_n_rates #print "p_rate, n_rate: ", np.average(p_rates), np.average(n_rates) def generate_prespikes(model, pre_rates): release_prob = model.release_prob Min = model.Min; Nin = model.Nin true_pre_spikes = []; pre_spikes = [] for j in range(Min): #Poisson pre-spikes true_pre_spikes.append( np.random.poisson( pre_rates[j] ) ) for i in range(Nin): jidx = model.preidx[i] pre_spikes.append( np.random.binomial( true_pre_spikes[jidx], release_prob ) ) return true_pre_spikes, pre_spikes def generate_inputs(model, uks, pre_spikes): release_prob = model.release_prob for i in range(model.Nin): if pre_spikes[i] > 0: model.Estimlist[i].interval = model.sduration/float(pre_spikes[i]) else: model.Estimlist[i].interval = model.sduration model.Estimlist[i].number = pre_spikes[i] model.Estimlist[i].start = model.sstart + (model.Estimlist[i].interval)*np.random.random() model.Estimlist[i].noise = model.snoise weight_tmp = model.gmax*uks[i]/release_prob #weight_tmp = model.gmax*uks[i] NC = h.NetCon(model.Estimlist[i], model.AMPAlist[i], 0, 0, weight_tmp) model.ncAMPAlist[i] = NC def generate_inh_inputs(model, true_pre_spikes): inhN = model.inhN Nin = model.Nin; Min = model.Min Erate = 0.0 for i in range(Nin): Erate += true_pre_spikes[ model.preidx[i] ]/float(Nin) Irate = ( float(Nin)/float(inhN) )*Erate Ispikes = np.random.poisson( Irate, (inhN) ) #print Ispikes for i in range(inhN): if Ispikes[i] > 0: model.Istimlist[i].interval = model.sduration/float(Ispikes[i]) else: model.Istimlist[i].interval = model.sduration model.Istimlist[i].number = Ispikes[i] model.Istimlist[i].start = model.sstart + (model.Istimlist[i].interval)*np.random.random() model.Istimlist[i].noise = model.snoise NC = h.NetCon(model.Istimlist[i], model.GABAlist[i], 0, 0, model.gI) model.ncGABAlist[i] = NC # ---------------------------------------------------------- # SIMULATION RUN def simulate(model): trec, vrec = h.Vector(), h.Vector() trec.record(h._ref_t) vrec.record(model.soma(0.5)._ref_v) h.celsius = model.CELSIUS #h.FInitializeHandler(1, initSpikes2) h.finitialize(model.E_PAS) neuron.run(model.tstop) return np.array(trec), np.array(vrec)