Olfactory bulb microcircuits model with dual-layer inhibition (Gilra & Bhalla 2015)

 Download zip file 
Help downloading and running models
Accession:153574
A detailed network model of the dual-layer dendro-dendritic inhibitory microcircuits in the rat olfactory bulb comprising compartmental mitral, granule and PG cells developed by Aditya Gilra, Upinder S. Bhalla (2015). All cell morphologies and network connections are in NeuroML v1.8.0. PG and granule cell channels and synapses are also in NeuroML v1.8.0. Mitral cell channels and synapses are in native python.
Reference:
1 . Gilra A, Bhalla US (2015) Bulbar microcircuit model predicts connectivity and roles of interneurons in odor coding. PLoS One 10:e0098045 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Olfactory bulb;
Cell Type(s): Olfactory bulb main mitral GLU cell; Olfactory bulb main interneuron periglomerular GABA cell; Olfactory bulb main interneuron granule MC GABA cell;
Channel(s): I A; I h; I K,Ca; I Sodium; I Calcium; I Potassium;
Gap Junctions:
Receptor(s): AMPA; NMDA; Gaba;
Gene(s):
Transmitter(s): Gaba; Glutamate;
Simulation Environment: Python; MOOSE/PyMOOSE;
Model Concept(s): Sensory processing; Sensory coding; Markov-type model; Olfaction;
Implementer(s): Bhalla, Upinder S [bhalla at ncbs.res.in]; Gilra, Aditya [aditya_gilra -at- yahoo -period- com];
Search NeuronDB for information about:  Olfactory bulb main mitral GLU cell; Olfactory bulb main interneuron periglomerular GABA cell; Olfactory bulb main interneuron granule MC GABA cell; AMPA; NMDA; Gaba; I A; I h; I K,Ca; I Sodium; I Calcium; I Potassium; Gaba; Glutamate;
# -*- coding: utf-8 -*-

import math, os
import pickle

from scipy import optimize
from pylab import *

## USAGE: python2.6 calc_differential_kernels.py

sys.path.extend(["..","../networks","../generators","../simulations"])

from networkConstants import * # has central_glom
from data_utils import *
from stimuliConstants import *

## study the main mitral that mit2 of glom1 (sideglom) connects to
## typically this is 0 i.e. mit0 of glom0 (mainglom)
fitted_mitral = DIRECTED_CONNS[2]

paramsnum = 12
STAtime = 1.0 # s

if paramsnum == 1:
    ## filelist with 5ms NOISEDT
    NOISEDT = 5e-3 # s
    NUMWHITETRAINS = 6
    rate_seednum = 191.0
    filelist = [
    '2011_06_22_20_00_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype-1.pickle',
    None,
    None,
    '2011_06_23_16_21_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype2.pickle',
    '2011_06_23_15_17_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype3.pickle',
    '2011_06_23_14_20_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype4.pickle',
    '2011_06_23_13_27_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype5.pickle'
    ]
elif paramsnum == 2:
    ## filelist with 25ms NOISEDT
    NOISEDT = 25e-3 # s
    NUMWHITETRAINS = 6
    rate_seednum = 191.0
    filelist = [
    '2011_06_23_17_49_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype-1.pickle',
    None,
    None,
    None,
    None,
    None,
    '2011_06_23_18_16_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype5.pickle'
    ]
elif paramsnum == 3:
    ## filelist with 5ms NOISEDT
    NOISEDT = 5e-3 # s
    NUMWHITETRAINS = 240
    rate_seednum = 191.0
    filelist = [
    '2011_06_24_12_02_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype-1.pickle',
    '2011_06_24_02_19_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype0.pickle',
    None,
    None,
    None,
    None,
    '2011_06_23_20_52_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype5.pickle'
    ]
elif paramsnum == 4:
    ## filelist with 25ms NOISEDT
    NOISEDT = 25e-3 # s
    NUMWHITETRAINS = 250
    rate_seednum = 191.0
    filelist = [
    None,
    None,
    None,
    [(191.0,'2011_06_24_15_36_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype2.pickle'),
    (181.0,'2011_06_25_13_44_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype2.pickle'),
    (171.0,'2011_06_27_20_15_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype2.pickle'),
    (161.0,'2011_06_28_00_31_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype2.pickle'),
    (151.0,'2011_06_28_12_13_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype2.pickle')],
    [(191.0,'2011_06_24_17_22_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype3.pickle'),
    (181.0,'2011_06_25_16_02_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype3.pickle')],
    None,
    [(191.0,'2011_06_24_13_35_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype5.pickle'),
    (181.0,'2011_06_25_19_13_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype5.pickle')]
    ]
elif paramsnum == 5:
    ## filelist with 50ms NOISEDT
    NOISEDT = 50e-3 # s
    NUMWHITETRAINS = 250
    filelist = [
    '2011_06_24_20_20_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype-1.pickle',
    None,
    None,
    '2011_06_24_18_50_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype2.pickle',
    None,
    None,
    None
    ]
elif paramsnum == 6:
    ## with extra noise sd=meanA, not sqrt(meanA)
    ## you'll need to rename the rate file
    ## filelist with 25ms NOISEDT
    NOISEDT = 25e-3 # s
    NUMWHITETRAINS = 250
    rate_seednum = 191.0
    filelist = [
    '2011_06_25_02_32_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype-1.pickle',
    None,
    '2011_06_25_00_46_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype1.pickle',
    '2011_06_24_23_25_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype2.pickle',
    None,
    None,
    None
    ]
elif paramsnum == 7:
    ## usual sqrt(meanA) noise
    ## filelist with 25ms NOISEDT
    NOISEDT = 25e-3 # s
    NUMWHITETRAINS = 250
    rate_seednum = 171.0

    ## without any interneurons
    #'2011_06_29_13_27_whitenoise_NOSINGLES_NOJOINTS_NOPGS_numgloms2_simtype-1.pickle'
    ## without PGs
    #'2011_06_29_12_51_whitenoise_SINGLES_JOINTS_NOPGS_numgloms2_simtype-1.pickle'
    ## all present
    #'2011_06_29_11_37_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype-1.pickle'
    ## only joints, 0.025 directed granules, mit2 to mit0
    #'2011_06_29_18_04_whitenoise_NOSINGLES_JOINTS_NOPGS_numgloms2_simtype2.pickle'
    ## only joints, 0.05 directed granules, mit2 and mit3 to mit0
    #'2011_06_29_19_46_whitenoise_NOSINGLES_JOINTS_NOPGS_numgloms2_simtype2.pickle'
    ## all present for averaging:
    #[(151.0,'2011_06_28_16_30_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype2.pickle'),
    #(161.0,'2011_06_28_18_43_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype2.pickle'),
    #(171.0,'2011_06_28_20_20_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype2.pickle')],

    filelist = [
    '2011_06_29_13_27_whitenoise_NOSINGLES_NOJOINTS_NOPGS_numgloms2_simtype-1.pickle',
    None,
    None,
    [(151.0,'2011_06_28_16_30_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype2.pickle'),
    (161.0,'2011_06_28_18_43_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype2.pickle'),
    (171.0,'2011_06_28_20_20_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype2.pickle')],
    '2011_06_28_22_44_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype3.pickle',
    '2011_06_29_00_36_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype4.pickle',
    None
    ]
elif paramsnum == 8:
    ## with extra noise sd=meanA, not sqrt(meanA)
    ## filelist with 100ms NOISEDT
    ## only joints, 0.05 directed granules, mit2 and mit3 to mit0, dt = 100ms
    #[(171.0,'2011_06_30_00_38_whitenoise_NOSINGLES_JOINTS_NOPGS_numgloms2_simtype2.pickle'),
    #(181.0,'2011_06_30_12_56_whitenoise_NOSINGLES_JOINTS_NOPGS_numgloms2_simtype2.pickle'),
    #(191.0,'2011_06_30_14_24_whitenoise_NOSINGLES_JOINTS_NOPGS_numgloms2_simtype2.pickle')],
    ## without interneurons, 0.05 directed granules, mit2 and mit3 to mit0, dt = 100ms
    #'2011_06_30_10_22_whitenoise_NOSINGLES_NOJOINTS_NOPGS_numgloms2_simtype2.pickle'
    NOISEDT = 100e-3 # s
    NUMWHITETRAINS = 250
    rate_seednum = 181.0
    filelist = [
    None,
    None,
    None,
    [(171.0,'2011_06_30_00_38_whitenoise_NOSINGLES_JOINTS_NOPGS_numgloms2_simtype2.pickle'),
    (181.0,'2011_06_30_12_56_whitenoise_NOSINGLES_JOINTS_NOPGS_numgloms2_simtype2.pickle'),
    (191.0,'2011_06_30_14_24_whitenoise_NOSINGLES_JOINTS_NOPGS_numgloms2_simtype2.pickle')],
    None,
    None,
    None
    ]
elif paramsnum == 9:
    ## with extra noise sd=meanA, not sqrt(meanA)
    ## filelist with 50ms NOISEDT
    NOISEDT = 100e-3 # s
    NUMWHITETRAINS = 250
    rate_seednum = 171.0
    filelist = [
    None,
    None,
    None,
    #[(191.0,'2011_06_30_23_58_whitenoise_NOSINGLES_JOINTS_NOPGS_numgloms2_simtype2.pickle'),
    #(171.0,'2011_07_01_11_42_whitenoise_NOSINGLES_JOINTS_NOPGS_numgloms2_simtype2.pickle')],
    #'2011_07_01_14_27_whitenoise_NOSINGLES_JOINTS_NOPGS_numgloms2_simtype2.pickle',
    #'2011_07_01_17_17_whitenoise_NOSINGLES_JOINTS_NOPGS_numgloms2_simtype2.pickle',
    [(191.0,'2011_07_02_01_03_whitenoise_NOSINGLES_JOINTS_NOPGS_numgloms2_simtype2.pickle')],
    None,
    ## 50ms
    #[(191.0,'2011_07_01_23_14_whitenoise_NOSINGLES_JOINTS_NOPGS_numgloms2_simtype4.pickle')],
    [(191.0,'2011_07_02_00_49_whitenoise_NOSINGLES_JOINTS_NOPGS_numgloms2_simtype4.pickle')],
    None
    ]
elif paramsnum == 10:
    ## filelist with 5ms NOISEDT
    NOISEDT = 5e-3 # s
    NUMWHITETRAINS = 250
    rate_seednum = 191.0
    noisemean = 5.0
    ## with extra noise sd=meanA, not sqrt(meanA)
    noisesd = sqrt(noisemean)
    
    ## Each element of filelist could be just a filenamestr,
    ## or [(rateseed,filenamestr),...] or
    ## [(rateseed,filenamestr,noisemean,noisesd),...]
    filelist = [
    #[(191.0,'2011_07_26_18_53_whitenoise_NOSINGLES_NOJOINTS_NOPGS_numgloms2_simtype-1.pickle',2.0,sqrt(2.0)),
    #(191.0,'2011_07_26_19_25_whitenoise_NOSINGLES_NOJOINTS_NOPGS_numgloms2_simtype-1.pickle',5.0,sqrt(5.0)),
    #(191.0,'2011_07_26_19_50_whitenoise_NOSINGLES_NOJOINTS_NOPGS_numgloms2_simtype-1.pickle',8.0,sqrt(8.0))],
    [(191.0,'2011_07_27_17_23_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype-1.pickle',2.0,sqrt(2.0)),
    (191.0,'2011_07_27_18_41_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype-1.pickle',5.0,sqrt(5.0)),
    (191.0,'2011_07_28_00_08_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype-1.pickle',8.0,sqrt(8.0))],
    None,
    #[(191.0,'2011_07_27_12_13_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype0.pickle',2.0,sqrt(2.0)),
    #(191.0,'2011_07_27_10_50_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype0.pickle',5.0,sqrt(5.0)),
    #(191.0,'2011_07_27_09_28_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype0.pickle',8.0,sqrt(8.0))],
    None,
    #[(191.0,'2011_07_27_13_19_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype2.pickle',2.0,sqrt(2.0))],
    None,
    None,
    None,
    None
    ]
elif paramsnum == 11:
    ## filelist with 25ms NOISEDT
    NOISEDT = 25e-3 # s
    NUMWHITETRAINS = 250
    rate_seednum = 191.0
    noisemean = 5.0
    ## with extra noise sd=meanA, not sqrt(meanA)
    noisesd = sqrt(noisemean)
    
    ## Each element of filelist could be just a filenamestr,
    ## or [(rateseed,filenamestr),...] or
    ## [(rateseed,filenamestr,noisemean,noisesd),...]
    filelist = [
    [(191.0,'2011_07_28_19_56_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype-1.pickle',2.0,sqrt(2.0)),
    (191.0,'2011_07_28_16_58_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype-1.pickle',5.0,sqrt(5.0)),
    (191.0,'2011_07_28_12_51_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype-1.pickle',8.0,sqrt(8.0))],
    None,
    None,
    None,
    None,
    None,
    None
    ]
elif paramsnum == 12:
    ## filelist with 5ms NOISEDT
    NOISEDT = 5e-3 # s
    NUMWHITETRAINS = 250
    rate_seednum = 191.0
    noisemean = 5.0
    ## with extra noise sd=meanA, not sqrt(meanA)
    noisesd = sqrt(0.1*noisemean/NOISEDT)
    
    ## Each element of filelist could be just a filenamestr,
    ## or [(rateseed,filenamestr),...] or
    ## [(rateseed,filenamestr,noisemean,noisesd),...]
    filelist = [
    [(191.0,'2011_08_01_20_27_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype-1.pickle',2.0,sqrt(0.1*2.0/NOISEDT)),
    (191.0,'2011_08_02_10_36_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype-1.pickle',5.0,sqrt(0.1*5.0/NOISEDT)),
    (191.0,'/2011_08_02_19_50_whitenoise_SINGLES_JOINTS_PGS_numgloms2_simtype-1.pickle',8.0,sqrt(0.1*8.0/NOISEDT))],
    None,
    None,
    None,
    None,
    None,
    None
    ]

STAtimelist = arange(-STAtime+NOISEDT,0.0+1e-10,NOISEDT)
kerneltimelist = arange(0.0,STAtime,NOISEDT)
lenSTA = (STAtime/NOISEDT)

def get_avg_STA(responseList,frateList):
    STAtotal = zeros(lenSTA)
    totalspikes = 0
    numtrials = 0
    ## same mean and sd of noise, but different firing rate realization
    for trainnum,frate in enumerate(frateList):
        ## same white noise firing rate realization, but different ORN spike times
        for trialnum,trialresponse in enumerate(responseList[:,trainnum]):
            numspikes,STAsum = calc_STA(frate,trialresponse,NOISEDT,STAtime)
            STAtotal += STAsum
            totalspikes += numspikes
            numtrials += 1
    ## only those response spikes are counted which occur after STAtime
    responserate = totalspikes/float(numtrials)/(PULSE_RUNTIME-STAtime)
    if totalspikes == 0:
        print "No target spikes"
        STA = STAtotal # essentially zeros if no target spikes
    else:
        STA = STAtotal/totalspikes
    return STA, responserate

def load_files(seednum, noisedt, numtrains, fn, fileidx):
    fratefile = '../generators/firerates/firerates_whitenoise_seed'\
        +str(seednum)+'_dt'+str(noisedt)+'_trains'+str(numtrains)+'.pickle'
    f = open(fratefile,'r')
    frateResponseList = pickle.load(f)
    f.close()
    print "Loaded the stimulus", fratefile
    ## frateResponseList[glomnum][trainnum]

    filenamefull = '../results/odor_whitenoise/'+fn
    f = open(filenamefull,'r')
    responseList = pickle.load(f)
    f.close()
    print "Loaded responses file",filenamefull
    ## responseList[trainnum][trialnum][mitralnum][spikenum]
    
    ## since last dimension has variable length, I need dtype=object below
    ## Now responseList[trialnum][trainnum][mitralnum] is a numpy array() of lists
    ## index as responseList[trialnum,trainnum,mitralnum][spikenum]
    responseList = array(responseList,dtype=object)

    if fileidx == 0:
        ## exc STA - input is to mainglom, study effect on mainmit
        frateList = frateResponseList[central_glom]
    else:
        ## inh STA - input is to sideglom, study effect on mainmit
        frateList = frateResponseList[central_glom+1]

    return frateList, responseList

def chisqfunc(params, STAdata_rev, ratedata):
    chisq = []
    kernel = params[0:lenSTA]
    kinteg = sum(kernel)*NOISEDT
    for i,(ratemean,ratesd,responserate) in enumerate(ratedata):
        diffs = (kernel+kinteg*ratemean**2)*(ratesd**2)/responserate - STAdata_rev[i]
        chisq.extend( [diff**2 for diff in diffs] )
    numdatapoints = len(STAdata_rev)*len(STAdata_rev[0])
    return array(chisq)/float(numdatapoints-lenSTA) # normalize to number of dof

if __name__ == "__main__":
    kernelslist = []
    fig = figure(facecolor='w') # 'none' is transparent
    ax = fig.add_subplot(111)
    fig2 = figure(facecolor='w') # 'none' is transparent
    ax2 = fig2.add_subplot(111)
    
    numfiles = len(filelist)-1
    for fileidx,filename in enumerate(filelist):
        if filename is None: continue
        if type(filename).__name__ == 'list':
            ## compute average STA for different seeds, but same mean, sd of noise
            if len(filename[0])==2:
                STAtotal = zeros(lenSTA)
                responserate_total = 0.0
                for seednum,fn in filename:
                    frateList,responseList = \
                        load_files(seednum, NOISEDT, NUMWHITETRAINS, fn, fileidx)
                    STA,responserate = get_avg_STA(responseList[:,:,fitted_mitral],frateList)
                    STAtotal += STA
                    responserate_total += responserate
                numSTAs = float(len(filename))
                STA = STAtotal/numSTAs
                responserate = responserate_total/numSTAs
                STAdata_rev = [ STA[::-1] ] # reversed STA
                ratedata = [(noisemean,noisesd,responserate)]
            ## compute average STA for diff seeds, diff means, diff sd-s of noise
            else:
                STAlist = []
                STAdata_rev = []
                ratedata = []
                for seednum,fn,meanHz,sdHz in filename:
                    frateList,responseList = \
                        load_files(seednum, NOISEDT, NUMWHITETRAINS, fn, fileidx)
                    STA,responserate = get_avg_STA(responseList[:,:,fitted_mitral],frateList)
                    STAlist.append( STA )
                    STAdata_rev.append( STA[::-1] ) # reversed STA
                    ratedata.append((meanHz,sdHz,responserate))

        else:
            frateList,responseList = \
                load_files(rate_seednum, NOISEDT, NUMWHITETRAINS, filename, fileidx)
            STA,responserate = get_avg_STA(responseList[:,:,fitted_mitral],frateList)
            STAdata_rev = [ STA[::-1] ] # reversed STA
            ratedata = [(noisemean,noisesd,responserate)]

        #### avg kernel from the reversed, calibrated STA
        kernel = zeros(lenSTA)
        for i,(noisemean,noisesd,responserate) in enumerate(ratedata):
            print i,noisemean, noisesd, responserate
            STArev = array(STAdata_rev[i])
            STArev = STArev - STArev.mean()
            thiskernel = STArev*responserate/noisesd**2
            kernel += thiskernel
            ax2.plot(kerneltimelist,thiskernel,label=str(i))
        kernel /= float(len(ratedata))

        #### fit THE KERNEL
        ## but I'm not using the fitted kernel, as it doesn't remove the offset :(
        ## initial params are an estimate of the kernel
        params0 = kernel
        ## You need more data points that params to be fitted!!
        ## That means at least two STA for different stim firing rates should be there.
        #if len(ratedata)>1:
        #    ## args is a tuple! if only one element write (elem, )
        #    params = optimize.leastsq( chisqfunc, params0,
        #        args=(STAdata_rev, ratedata), full_output=1 )
        #    print params[3]
        #    params = params[0] # leastsq returns a whole tuple of stuff - errmsg etc.
        #    kernel = params[0:lenSTA]
        kernelslist.append((fileidx,kernel))

        if fileidx == 0:
            ## exc STA - input is to mainglom, study effect on mainmit
            color = (0,0,0)
            legend = 'exc'
        else:
            ## inh STA - input is to sideglom, study effect on mainmit
            color = ( fileidx/float(numfiles), 1-fileidx/float(numfiles), 0 )
            legend = 'inh @ '+str(varied_mainrate[fileidx-1])+'Hz'
        ax.plot(STAtimelist,STA,color=color,label=legend)
        ax2.plot(kerneltimelist,kernel,color=color,label=legend)
 
    biglegend("upper left",ax)
    biglegend("upper right",ax2)
    axes_labels(ax,'time (s)','mean rate of 400 ORNs (Hz)')
    axes_labels(ax2,'time (s)','mitral kernel')
    ax.set_title('exc and inh STAs', fontsize=24)
    ax2.set_title('exc and inh kernels', fontsize=24)
       
    filename = '../results/odor_whitenoise/kernels.pickle'
    kernelsfile = open(filename,'w')
    pickle.dump(kernelslist, kernelsfile)
    kernelsfile.close()
    print "wrote",filename
    
    show()

Loading data, please wait...