A model of optimal learning with redundant synaptic connections (Hiratani & Fukai 2018)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:225075
This is a detailed neuron model of non-parametric near-optimal latent model acquisition using multisynaptic connections between pre- and postsynaptic neurons.
Reference:
1 . Hiratani N, Fukai T (2018) Redundancy in synaptic connections enables neurons to learn optimally. Proc Natl Acad Sci U S A 115:E6871-E6879 [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:
Cell Type(s): Neocortex V1 L2/6 pyramidal intratelencephalic cell;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON; Python;
Model Concept(s): Synaptic Plasticity;
Implementer(s): Hiratani,Naoki [N.Hiratani at gmail.com];
Search NeuronDB for information about:  Neocortex V1 L2/6 pyramidal intratelencephalic cell;
/
HirataniFukai2018
data
README.html
ca.mod *
cad.mod *
caL3d.mod *
CaT.mod
exp2synNMDA.mod
h.mod *
HH2.mod *
kca.mod *
kir.mod *
km.mod *
kv.mod *
na.mod *
L23.hoc
libcell.py
md_readout.py
neuron_simulation.py
screenshot.png
                            
# ----------------------------------------------------------
# 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)


Loading data, please wait...