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]
Citations  Citation Browser
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;
import sys, pickle

sys.path.extend(["..","../networks","../simulations"])
from networkConstants import *
from stimuliConstants import *
from simset_odor import *

from moose_utils import *
from poisson_utils import *
from data_utils import * # has mpi import and variables also

from pylab import * # part of matplotlib that depends on numpy but not scipy

### generate firefiles i.e. list of spike times from the firerates computed by generate_firerates.py
############## For generating random pulses -- takes maximum time, hence parallelized
## from node000 (not from gj, else master process on node000 cannot connect to display for graphs):
## mpiexec -machinefile ~/hostfile -n <numpulses*NUM_GLOMS+1> ~/Python-2.6.4/bin/python2.6 generate_firefiles_odors.py
## typical value for numpulses = RANDOM_PULSE_NUMS*2
## typically for 10 gloms,
## mpiexec -machinefile ~/hostfile -n 61 ~/Python-2.6.4/bin/python2.6 generate_firefiles_odors.py
############## For generating odor morphs, etc. not parallelized
## python2.6 generate_firefiles_odors.py

## REALRUNTIME is from simset_odor.py
RUNTIME = REALRUNTIME + SETTLETIME

## binning to plot odor responses
bindt = (RUNTIME-SETTLETIME)/respbins
## to make bins causal, keep bin_t at right edge of bin
tlist = arange(SETTLETIME+bindt,RUNTIME+bindt/2.0,bindt)
pulsetlist = arange(pulsebindt,PULSE_RUNTIME+pulsebindt/2.0,pulsebindt)

# time points for the firing rate which is read from a pickled file
firingtsteps = arange(0,RUNTIME+1e-10,FIRINGFILLDT)# include the last RUNTIME point also.
numt = len(firingtsteps)
extratime = arange(0,2*RESPIRATION,FIRINGFILLDT)
pulsetsteps = arange(0,PULSE_RUNTIME+1e-10,FIRINGFILLDT)
numtpulse = len(pulsetsteps)

stim_firefiles_dirname = '../firefiles/firefiles'+str(stim_rate_seednum)
if NONLINEAR_ORNS: stim_firefiles_dirname += '_NL'+NONLINEAR_TYPE

def logisticfn(x,offset,width):
    return 1.0/(1.0+exp(-(x-offset)/width))

def thresholded_sigmoid(response):
    ## take a numpy array and threshold and sigmoid() it
    response[ where(response<=0)[0] ] = 0 # clip negative values to 0

    ## I want a sigmoid for only non-negative Reals
    ## Rather than a truncated sigmoid,
    ## I use the cumulative distribution function of the gamma distribution (non-negative domain).
    ## Gamma distribution with scale = 9.0 and shape = 0.5 is close to Gaussian.
    ## Its integral i.e. CDF is a 'sigmoid' over the non-negative Reals.
    #### No I can't use it as it requires scipy.stats.gamma and it is not available on cluster nodes
    #response = scipy.stats.gamma.cdf(response,[9.0]*scipy.stats.gamma.numargs,loc=0,scale=0.5)
    
    ## erf(0.5)~=0.5; set NONLINEAR_CENTER as the center of the sigmoid.
    ## Width should be normal NONLINEAR_WIDTH (spans 2-3 orders of mag in concentration).
    #### No I can't use it as it requires scipy.stats.gamma and it is not available on cluster nodes
    #response = scipy.special.erf((response-NONLINEAR_CENTER)/NONLINEAR_WIDTH)
    ## set the response to zero where negative
    #response[ where(response<=0)[0] ] = 0

    ## finally since gammaCDF and erf can't be used, we go for logistic function.
    ## 1.5*FIRINGMEANA is the center of the sigmoid, FIRINGMEANA/2.0 is width
    response = logisticfn(response,NONLINEAR_CENTER,NONLINEAR_WIDTH)

    ## peak ORN firing set as 3*FIRINGMEANA
    return 3*FIRINGMEANA*response

def NLcondition(frate,glomnum):
    if not NONLINEAR_ORNS: return frate
    if NONLINEAR_TYPE == 'P':
        if glomnum == 0: return thresholded_sigmoid(frate)
        else: return frate
    else:
        if glomnum == 0: return frate
        else: return thresholded_sigmoid(frate)

def odor_and_air_stimuli():
    ## mitral and PG odor ORNs firing files
    seed([stim_rate_seednum]) # Seed numpy's random number generator.
    for glomnum in range(NUM_GLOMS):
        if SHOWFIGS:
            figure()
            title(str(glomnum))
            ylabel('Hz')
            xlabel('time (s)')
        frateperglomList = []
        for odoridx,(odorA, odorB) in enumerate(inputList):
            frate = frateOdorList[glomnum][odoridx]
            frate = NLcondition(frate,glomnum) # output non-linearity if switch is set
            for avgnum in range(MAXNUMAVG):
                mitralfirefilename = stim_firefiles_dirname+\
                    '/firetimes_2sgm_glom_'+str(glomnum)+\
                    '_odor_'+str(odorA)+'_'+str(odorB)+'_avgnum'+str(avgnum)
                ornstimvector_merged = write_odor_files(NUM_ORN_FILES_PER_GLOM, frate,\
                    mitralfirefilename, RUNTIME, firingtsteps)
            if SHOWFIGS:
                ## plotBins here returns firing rate of NUM_ORN_FILES_PER_GLOM combined.
                ## So divide by NUM_ORN_FILES_PER_GLOM.
                ## We just plot ornstimvector_merged for the last avgnum 
                ratebins = [rate/NUM_ORN_FILES_PER_GLOM \
                    for rate in plotBins(ornstimvector_merged, respbins, RUNTIME, SETTLETIME)]
                plot(tlist, ratebins, color=(odorA,odorB,0), marker=',')
            ## extra unmodelled sister cell excitation to granules.
            ## No need of separate files, above files will do for sister excitation
            #mitralfirefilename = '../firefiles/firetimes_2exp_sisters_glom_'+str(glomnum)+\
            #    '_odor_'+str(odorA)+'_'+str(odorB)+'_avgnum'+str(avgnum)
            #write_odor_files(NUM_ORN_FILES_PER_GLOM, frate, mitralfirefilename)

def random_pulse_stimuli(glomnum, pulse_i):
    numpulses = 2*RANDOM_PULSE_NUMS
    ## mitral and PG odor ORNs firing files
    frate_air = randomResponseList[glomnum][0]
    ## the first 'frate' in randomPulseList[glomnum] is for air (rectangle pulse),
    ## the second is a random pulsed air pulse.
    ## then alternately odorA and odorB.
    ## the last odorA and odorB frates are combined and given simultaneously.
    frate = randomResponseList[glomnum][pulse_i+1]
    ## randomPulseList[glomnum][[pulses,frate],...]
    ## pulse_i goes from 0 to numpulses-1 as the
    ## mpi slave allocator skips the first (rectangle air)
    ## and last (combined with 2nd last) pulse
    if pulse_i == 0:
        ## this a just a random air pulse sequence,
        ## so rectangular air pulse 'frate_air' is not added
        frate = frate + ORN_BGND_MEAN
    elif pulse_i < (numpulses-1):
        ## the ORN background is added here, and not used in
        ## fitting the kernel in generate_firerates_odors.py,
        ## as dF/F is not sensistive to background.
        frate = frate + frate_air + ORN_BGND_MEAN
    else:
        ## the last odorA and odorB frates are combined and given simultaneously.
        ## the ORN background is added here, and not used to fit the kernel,
        ## as dF/F is not sensistive to background.
        frate = frate + randomResponseList[glomnum][pulse_i+2] + frate_air + ORN_BGND_MEAN
    frate = NLcondition(frate,glomnum) # output non-linearity if switch is set
    for avgnum in range(MAXNUMAVG):
        mitralfirefilename = stim_firefiles_dirname+\
            '/firetimes_rndpulse_glom_'+str(glomnum)+\
            '_pulse_'+str(pulse_i)+'_avgnum'+str(avgnum)
        ornstimvector_merged = write_odor_files(NUM_ORN_FILES_PER_GLOM, frate,\
            mitralfirefilename, PULSE_RUNTIME, pulsetsteps)
    ## plotBins here returns firing rate of NUM_ORN_FILES_PER_GLOM combined.
    ## So divide by NUM_ORN_FILES_PER_GLOM.
    ## We just plot ornstimvector_merged for the last avgnum 
    ratebins = [rate/NUM_ORN_FILES_PER_GLOM \
        for rate in plotBins(ornstimvector_merged, pulsebins, PULSE_RUNTIME, 0.0)]
    ## extra unmodelled sister cell excitation to granules.
    ## No need of separate files, above files will do for sister excitation
    #mitralfirefilename = '../firefiles/firetimes_2exp_sisters_glom_'+str(glomnum)+\
    #    '_odor_'+str(odorA)+'_'+str(odorB)+'_avgnum'+str(avgnum)
    #write_odor_files(NUM_ORN_FILES_PER_GLOM, frate, mitralfirefilename)
    return ratebins

def interglom_stimuli():
    """
    This reproduces SA small-world excitatory network to PGs via ETs.
    I HAVE NOT YET INCLUDED THE BACKGROUND DUE TO MEAN RESPIRATORY ACTIVITY
    """
    seed([600.0]) ##### Seed numpy's random number generator. If no parameter is given, it uses current system time
    # before taking mean activity over all glomeruli, we should clip negative values
    # because SA cells only receive _excitatory_ input from ET cells.
    frateOdorListClipped = clip(frateOdorList,0.0,1e6)
    # mean activity over gloms (axis=0) as a function of time
    fratemean = mean(frateOdorListClipped,axis=0)
    # Short axon network delays and neuronal integration is modelled as 25ms integrated excitation to PGs
    # the delay should be put in as part of the synapse
    integrate_window_length = int(SA_integrate_time/FIRINGFILLDT)
    # Make a normalized integration window
    # ones() makes array of floats by default, so no integer division problems
    rect = ones((integrate_window_length,))/integrate_window_length
    figure()
    title('SAs')
    ylabel('Hz')
    xlabel('time (s)')
    for odoridx,(odorA,odorB) in enumerate(inputList):
        # convolution takes moving average,
        # I must use mode='full' and take [0:numt], else the SA->(ET->)PG excitation is acausal.
        # of course, the initial part from t = 0 till 'SA_integrate_time' is invalid, but that is within SETTLETIME.
        interglom_rate = convolve(fratemean[odoridx], rect, mode='full')[0:numt]
        SAfirefilename = stim_firefiles_dirname+\
            +'/firetimes_SA_odor_'+str(odorA)+'_'+str(odorB)
        ornstimvector_merged = write_odor_files(NUM_ORN_FILES_PER_GLOM,
            interglom_rate, SAfirefilename, RUNTIME, firingtsteps)
        ratebins = [rate/NUM_ORN_FILES_PER_GLOM \
        for rate in plotBins(ornstimvector_merged, respbins, RUNTIME, SETTLETIME)]
        plot(tlist, ratebins, color=(odorA,odorB,0), marker=',')

def pulseinresp_stimuli(_showfigs):
    seed([stim_rate_seednum]) ##### Seed numpy's random number generator. If no parameter is given, it uses current system time
    stim_frate_filename = '../generators/firerates/firerates_pulseinresp_'+str(stim_rate_seednum)
    stim_frate_filename += '.pickle'
    f = open(stim_frate_filename,'r')
    pulseInRespList,kernels = pickle.load(f)
    f.close()

    pulseinresp_RUNTIME = (NUM_RESPS+1)*RESPIRATION
    pulseinresp_tsteps = arange(0,pulseinresp_RUNTIME+1e-10,FIRINGFILLDT)
    pulseinresp_bintsteps = arange(0,pulseinresp_RUNTIME,pulsebindt)

    for glomnum in range(NUM_GLOMS):
        (_,frates) = pulseInRespList[glomnum]
        if glomnum==0 and _showfigs: figure()
        for ratenum,frate in enumerate(frates):
            frate = NLcondition(frate,glomnum) # output non-linearity if switch is set
            for avgnum in range(MAXNUMAVG):
                mitralfirefilename = stim_firefiles_dirname+\
                    '/firetimes_pulseinresp_glom_'+str(glomnum)+\
                    '_odor'+['R','A','B'][ratenum]+'_avgnum'+str(avgnum)
                ornstimvector_merged = write_odor_files(NUM_ORN_FILES_PER_GLOM, frate,\
                    mitralfirefilename, pulseinresp_RUNTIME, pulseinresp_tsteps)
            ## plotBins here returns firing rate of NUM_ORN_FILES_PER_GLOM combined.
            ## So divide by NUM_ORN_FILES_PER_GLOM.
            ## We just plot ornstimvector_merged for the last avgnum 
            pulseinrespbins = int((pulseinresp_RUNTIME/pulsebindt))
            ratebins = [rate/NUM_ORN_FILES_PER_GLOM \
                for rate in plotBins(ornstimvector_merged, pulseinrespbins, pulseinresp_RUNTIME, 0.0)]
            if glomnum==0 and _showfigs:
                plot(pulseinresp_bintsteps,ratebins,['k','r','b'][ratenum])

def scaledpulses_stimuli(_showfigs):
    seed([stim_rate_seednum]) ##### Seed numpy's random number generator. If no parameter is given, it uses current system time
    stim_frate_filename = '../generators/firerates/firerates_scaledpulses_width'+\
        str(scaledWidth)+'_'+str(stim_rate_seednum)
    stim_frate_filename += '.pickle'
    f = open(stim_frate_filename,'r')
    scaledPulses,kernels = pickle.load(f)
    f.close()

    scaledpulses_tsteps = arange(0,SCALED_RUNTIME+1e-10,FIRINGFILLDT)
    scaledpulses_bintsteps = arange(0,SCALED_RUNTIME,pulsebindt)

    for glomnum in range(NUM_GLOMS):
        frates = scaledPulses[glomnum]
        if glomnum==0 and _showfigs: figure()
        for scalenum,scale in enumerate(scaledList):
            for fratenum,frate in enumerate(frates[scalenum]):
                frate = NLcondition(frate,glomnum) # output non-linearity if switch is set
                for avgnum in range(MAXNUMAVG):
                    mitralfirefilename = stim_firefiles_dirname+\
                        '/firetimes_scaledpulses_width'+str(scaledWidth)+'_glom_'+str(glomnum)+\
                        '_odor'+['A','B'][fratenum]+'_scale'+str(scalenum)+'_avgnum'+str(avgnum)
                    ornstimvector_merged = write_odor_files(NUM_ORN_FILES_PER_GLOM, frate,\
                        mitralfirefilename, SCALED_RUNTIME, scaledpulses_tsteps)
                if glomnum==0 and _showfigs:
                    ## plotBins here returns firing rate of NUM_ORN_FILES_PER_GLOM combined.
                    ## So divide by NUM_ORN_FILES_PER_GLOM.
                    ## We just plot ornstimvector_merged for the last avgnum 
                    scaledpulsebins = int((SCALED_RUNTIME/pulsebindt))
                    ratebins = [rate/NUM_ORN_FILES_PER_GLOM \
                        for rate in plotBins(ornstimvector_merged, scaledpulsebins, SCALED_RUNTIME, 0.0)]
                    plot(scaledpulses_bintsteps[:len(ratebins)],ratebins)

def interglom_pulsestimuli():
    """
    This reproduces SA small-world excitatory network to PGs via ETs.
    I HAVE NOT YET INCLUDED THE BACKGROUND DUE TO MEAN RESPIRATORY ACTIVITY
    """
    seed([800.0]) ##### Seed numpy's random number generator. If no parameter is given, it uses current system time
    # before taking mean activity over all glomeruli, we should clip negative values
    # because SA cells only receive _excitatory_ input from ET cells.
    # randomPulseList[glomnum][[pulses,frate],...]
    numpulses = len(randomPulseList[0])
    # the first pulse is air, then A and B alternate, the last has to be combined with the second last.
    frateList = empty([NUM_GLOMS,numpulses-2,len(randomPulseList[0][0][1])])
    for glomnum in range(NUM_GLOMS):
        frate_air = randomPulseList[glomnum][0][1]
        for pulsenum in range(1,numpulses-1):
            if pulsenum == 1:
                frateList[glomnum][pulsenum-1] =\
                    array( randomPulseList[glomnum][pulsenum][1] + ORN_BGND_MEAN )
            elif pulsenum < numpulses-2:
                frateList[glomnum][pulsenum-1] =\
                    array( randomPulseList[glomnum][pulsenum][1] + frate_air + ORN_BGND_MEAN)
            else:
                frateList[glomnum][pulsenum-1] = array( randomPulseList[glomnum][pulsenum+1][1] + \
                    randomPulseList[glomnum][pulsenum+1][1] + frate_air + ORN_BGND_MEAN)
    fratePulseListClipped = clip(frateList,0.0,1e6)
    # mean activity over gloms (axis=0) as a function of time
    fratemean = mean(fratePulseListClipped,axis=0)
    # Short axon network delays and neuronal integration is modelled as 25ms integrated excitation to PGs
    # the delay should be put in as part of the synapse
    integrate_window_length = int(SA_integrate_time/FIRINGFILLDT)
    # Make a normalized integration window
    # ones() makes array of floats by default, so no integer division problems
    rect = ones((integrate_window_length,))/integrate_window_length
    figure()
    title('SAs Pulses')
    ylabel('Hz')
    xlabel('time (s)')
    for pulse_i in range(numpulses-2):
        ## convolution takes moving average,
        ## I must use mode='full' and take [0:numt], else the SA->(ET->)PG excitation is acausal.
        ## of course, the initial part from t = 0 till 'SA_integrate_time' is invalid,
        ## but that is within SETTLETIME.
        interglom_rate = convolve(fratemean[pulse_i], rect, mode='full')[0:numtpulse]
        SAfirefilename = stim_firefiles_dirname+\
            +'/firetimes_SA_rndpulse'+'_pulse_'+str(pulse_i)
        ornstimvector_merged = write_odor_files(NUM_ORN_FILES_PER_GLOM,\
            interglom_rate, SAfirefilename, PULSE_RUNTIME, pulsetsteps)
        ratebins = [rate/NUM_ORN_FILES_PER_GLOM\
            for rate in plotBins(ornstimvector_merged, pulsebins, PULSE_RUNTIME, 0.0)]
        plot(pulsetlist, ratebins, color=(pulse_i/float(numpulses),1-pulse_i/float(numpulses),0), marker=',')

def write_odor_files(numfiles, frate, firefilename, runtime, tsteps, vary=None):
    frate = clip(frate,0.0,1e6)
    PEAK_FIRING_RATE = max(frate)*2
    if PEAK_FIRING_RATE < 1.0: PEAK_FIRING_RATE = 1.0
    ornstimvector_merged = []
    firefile = open(firefilename+'.txt','w')
    for i in range(numfiles):
        ## if variation is called for, then vary the firing rate by a factor for each file
        if not vary is None: frate_rand = frate * normal(vary[0],vary[1])/vary[0]
        else: frate_rand = frate
        ornstimvector = poissonTrainVaryingRate(runtime,PEAK_FIRING_RATE,REFRACTORY,tsteps,frate_rand)
        firefile.write(' '.join([str(t) for t in ornstimvector])+'\n')
        ornstimvector_merged.extend(ornstimvector)
    firefile.close()
    ornstimvector_merged.sort()
    print "wrote ", firefilename
    return ornstimvector_merged

if __name__ == "__main__":
    ### Seed only if called directly, else do not seed.
    ### Also seeding this way ensures seeding after importing other files that may set seeds.
    ### Thus this seed overrides other seeds.
    #seed([600.0]) ##### Seed numpy's random number generator. If no parameter is given, it uses current system time

    if 'NOSHOW' not in sys.argv: SHOWFIGS=True
    else: SHOWFIGS=False

    stim_frate_filename = '../generators/firerates/firerates_2sgm_'+str(stim_rate_seednum)
    stim_frate_filename += '.pickle'
    f = open(stim_frate_filename,'r')
    frateOdorList,fratePulseList,randomPulseList,\
    randomPulseStepsList,randomResponseList,kernels\
        = pickle.load(f)
    f.close()

    if mpisize == 1:
        ## seeds are set in each of the below functions:
        ## so you may comment / uncomment below functions
        ## depending on what you want to generate.
        #interglom_stimuli()
        if 'PULSEINRESP' in sys.argv: pulseinresp_stimuli(SHOWFIGS)
        elif 'SCALEDPULSES' in sys.argv: scaledpulses_stimuli(SHOWFIGS)
        else: odor_and_air_stimuli()
        #interglom_pulsestimuli()
        if SHOWFIGS: show()
    else:
        ## Random pulses, not doing two pulse stimulus
        ## since random pulses convey all the info.
        numpulses = RANDOM_PULSE_NUMS*2
        numpulses_fl = float(numpulses)
        if mpirank == boss: # boss collates
            for glomnum in range(NUM_GLOMS):
                if SHOWFIGS:
                    figure()
                    title("Pulses for glom "+str(glomnum))
                    ylabel('Hz')
                    xlabel('time (s)')
                for pulsenum in range(numpulses):
                    procnum = glomnum*numpulses + pulsenum + 1
                    print 'waiting for process '+str(procnum)+'.'
                    ratebins = mpicomm.recv(source=procnum, tag=0)
                    if SHOWFIGS:
                        plot(pulsetlist, ratebins, \
                            color=(pulsenum/numpulses_fl,1-pulsenum/numpulses_fl,0), marker=',')
                    print 'received from process '+str(procnum)+'.'
            if SHOWFIGS: show()
        else:
            seed([500.0*mpirank]) ##### Seed numpy's random number generator. If no parameter is given, it uses current system time
            glomnum = (mpirank-1)/numpulses
            pulsenum = (mpirank-1)%numpulses
            ratebins = random_pulse_stimuli(glomnum,pulsenum)
            mpicomm.send( ratebins, dest=boss, tag=0 )