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;
#!/usr/bin/env python
# -*- coding: utf-8 -*-

import os
import sys
import math
import datetime
import pickle

## from node000:
## mpiexec -machinefile ~/hostfile -n <numavgs*numodors+1> ~/Python-2.6.4/bin/python2.6 odor_morphs.py
## nohup mpiexec -machinefile ~/hostfile -n 57 ~/Python-2.6.4/bin/python2.6 odor_morphs.py < /dev/null &
## typical value for numavgs = 8 
## (depends on number of available processing nodes and number of odorfiles generated)
## typical value for numodors = len(inputList) = 7.
## OR for a single odor run; from any node (or laptop -- can use system python (2.7)):
## python2.6 odor_morphs.py [SAVEFORMOVIE]
## Set various option like NO_PGs or ONLY_TWO_MITS in simset_odor

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

from moose_utils import * # imports moose
from data_utils import * # has mpi import and variables also
from OBNetwork import *
from sim_utils import *

from stimuliConstants import * # has SETTLETIME, inputList and pulseList, GLOMS_ODOR, GLOMS_NIL
from simset_odor import * # has REALRUNTIME, NUMBINS

RUNTIME = REALRUNTIME + SETTLETIME

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

#-----------------------------------------------------------

class odorResponse:
    
    def __init__(self,mpirank=mpirank): # mpirank is defined in data_utils.py
        self.mpirank = mpirank
        self.context = moose.PyMooseBase.getContext()

    def setupStim(self,network,args,avgnum):
        odorA, odorB = args[0]
        self.setupOdor(network, odorA, odorB, avgnum)        
        print "Setup odorA =",odorA,"odorB =",odorB,"at",self.mpirank

    def setupOdor(self, network, odorA, odorB, avgnum):
        ### first figure out which PG belongs to which glom
        ### PG_glom_map[pgname] returns the glom num of the PG: needed for ORN to PG connections.
        PG_glom_map = {}
        for projname in network.projectionDict.keys():
            if 'PG_mitral' in projname:
                for i,proj in enumerate(network.projectionDict[projname][2]):
                    # get the glomnum from the post path proj[2]
                    # name of the mitral cell from '/mitrals_2/...'
                    mitname = string.split(proj[2],'/')[1]
                    # glomerulus number from 'mitrals_2' by integer division i.e. 2/2 = 1
                    glomnum = int(string.split(mitname,'_')[1]) / 2
                    # name of the PG cell from '/PGs_2/...'
                    pgname = string.split(proj[1],'/')[1]
                    PG_glom_map[pgname] = glomnum
        ### Now connect the ORNs
        for projname in network.projectionDict.keys():
            #### Calling attach_spikes() for each projection,
            #### would reconnect files to the same segment multiple times.
            #### But attach_files_uniquely() checks whether timetable.tableSize is zero or not
            #### i.e. files already attached or not.
            ### connect ORNs to mitrals
            if 'ORN_mitral' in projname:
                print "connecting ORN files to mitrals"
                for i,proj in enumerate(network.projectionDict[projname][2]):
                    # get the glomnum from the post path proj[2]
                    mitname = string.split(proj[2],'/')[1] # name of the mitral cell from '/mitrals_2/...'
                    glomnum = int(string.split(mitname,'_')[1]) / 2 # glomerulus number from 'mitrals_2' by integer division i.e. 2/2 = 1
                    filebase = ORNpathseedstr+'firetimes_2sgm_glom_'+str(glomnum)
                    self.attach_files_uniquely(filebase,proj[0],proj[2],odorA,odorB,avgnum)
            ### connect ORNs to PG
            if 'ORN_PG' in projname:
                print "connecting ORN files to PGs"
                for i,proj in enumerate(network.projectionDict[projname][2]):
                    pgname = string.split(proj[2],'/')[1] # name of the PG cell from '/PGs_2/...'
                    glomnum = PG_glom_map[pgname]
                    filebase = ORNpathseedstr+'firetimes_2sgm_glom_'+str(glomnum)
                    self.attach_files_uniquely(filebase,proj[0],proj[2],odorA,odorB,avgnum)
            ### connect SAs to PG
            if 'SA_PG' in projname:
                print "connecting SA files to PGs"
                for i,proj in enumerate(network.projectionDict[projname][2]):
                    pgname = string.split(proj[2],'/')[1] # name of the PG cell from '/PGs_2/...'
                    glomnum = PG_glom_map[pgname]
                    filebase = ORNpathseedstr+'firetimes_SA'
                    self.attach_files_uniquely(filebase,proj[0],proj[2],odorA,odorB)
            ###### I am back to 'extra-connecting' modelled mitral as extra sister mitrals excitation to granules
            ###### Previously, as below, I was connecting ORNs of the glom to granules
            ###### which caused inhibition even when the sister mitrals were not even firing!
            ### connect unmodelled extra sister mitrals as files to granules
            #if 'mitral_granule_extra' in projname:
            #    print "Connecting unmodelled sister excitation files to granules"
            #    for i,proj in enumerate(network.projectionDict[projname][2]):
            #        granulename = string.split(proj[2],'/')[1] # name of the granule cell from '/granules_singles_2/...'
            #        # glomnum from pre_path = proj[1] = 'file[+<glomnum>]_<filenumber1>[_<filenumber2>...]'
            #        glomstr = proj[1].split('+')[1].split('_',1)[0]
            #        filebase = ORNpathseedstr+'firetimes_2sgm_glom_'+glomstr
            #        self.attach_files_uniquely(filebase,proj[0]+'_'+glomstr,proj[2],odorA,odorB,avgnum)

    def attach_files_uniquely(self,filebase,synname,postsegpath,odorA,odorB,avgnum=None):
        ttpath = postsegpath+'/'+synname+'_tt'
        if self.context.exists(ttpath):
            # timetable already created by networkML reader - just wrap it below.
            tt = moose.TimeTable(ttpath) # post_segment_path+'/'+syn_name+'_tt'
        else:
            ## if timetable was not already created by networkML reader,
            ## it means that the synaptic weights must be zero!
            ## (no extra inhibition - only main inhibition)
            ## hence do not attach spikefiles
            return
        if tt.tableSize != 0: return # if files are already attached, do nothing!
        filebase += '_odor_'+str(odorA)+'_'+str(odorB)
        if avgnum is not None: filebase += '_avgnum'+str(avgnum)
        ## attach_spikes() accesses filenumbers to this segment
        ## from 'fileNumbers' field (of the timetable object in MOOSE)
        ## which is created while reading in networkML.
        attach_spikes(filebase, tt, uniquestr+str(self.mpirank))

    def run(self,network, binned):
        print "Resetting MOOSE."
        # from moose_utils.py sets clocks and resets
        resetSim(network.context, SIMDT, PLOTDT)
        print "Running at",self.mpirank
        network.context.step(RUNTIME)
        mitral_responses = []
        mitral_responses_binned = []
        if ONLY_TWO_MITS or NO_LATERAL: num_mits = MIT_SISTERS
        else: num_mits = NUM_GLOMS*MIT_SISTERS
        for mitnum in range(num_mits):
            mitral = network.mitralTable[mitnum]
            ## only the last respiration cycle is taken
            if binned: mitral_responses_binned.append(
                plotBins(mitral._vmTableSoma, NUMBINS, RUNTIME,\
                (NUM_RESPS-1)*RESPIRATION+SETTLETIME) )
            ## need to convert to numpy's array(),
            ## else MOOSE table cannot be pickled for mpi4py send()
            mitral_responses.append(array(mitral._vmTableSoma))
        return (mitral_responses,mitral_responses_binned)

#----------------------------------------------------------------

if __name__ == "__main__":

    ## uniquestr to put in every temp filename to avoid clashing with other processes
    if len(sys.argv)>2: uniquestr = sys.argv[2]+'_' # _ necessary, else say 'morphs2'+mpirank is screwed
    else: uniquestr = 'morphs_'
    
    numodors = len(inputList)
    #### if only one process is called, plot one odor directly
    if mpisize == 1:
        ## for ref: inputList = [ (0.0,1.0), (0.2,0.8), (0.4,0.6), (0.6,0.4), (0.8,0.2), (1.0,0.0), (0.0,0.0) ]
        odorA = 0.0
        odorB = 1.0
        avgnum = 0
        ## OBNetwork.py uses the proxympirank passed here
        proxympirank = avgnum*numodors + 0 + 1 # avgnum*numodors + odornum + 1
        sim =  odorResponse(proxympirank)
        ## 'PG' includes 'ORN_PG', 'PG_mitral', 'mitral_PG' and 'SA_PG'
        if ONLY_TWO_MITS and not NO_PGS: includeProjections = ['PG']
        else: includeProjections = []
        tweaks = build_tweaks(CLUB_MITRALS, NO_SPINE_INH, NO_SINGLES,\
            NO_JOINTS, NO_MULTIS, NO_PGS, ONLY_TWO_MITS,\
            includeProjections=includeProjections, nolateral=NO_LATERAL)
        BINNED = True#False # for mitrals only
        ## if not BINNED, save the full mitral Vm-s
        ## and not just their spiketimes by setting spiketable = False below.
        network = OBNetwork(OBNet_file, synchan_activation_correction, tweaks,\
            proxympirank, 'morphs', granfilebase, spiketable=BINNED)
        #printNetTree() # from moose_utils.py
        ## monitor those interneurons that are connected to mitral indices 0 and 1
        ## save only spiketimes by setting extras_spikes_only=True
        extras_spikes_only = True # for interneurons
        ## choose one of the below for only interneurons connected to mits0/1 vs all.
        interneurons_args = {}
        #extras_args = {'mitrals':[0,1]}
        tables = setupTables(network, NO_PGS, NO_SINGLES, NO_JOINTS, NO_MULTIS,\
            interneurons_args, spikes=extras_spikes_only)
        ## To watch the pre compartment of mit2 that inhibits soma of mit 1
        #mit2 = moose.Cell('/mitrals_2')
        #mit2.precomp = moose.Compartment(get_matching_children(mit2, ['Seg0_sec_dendd4_4_278'])[0])
        #mit2._vmTablePrecomp = setupTable("vmTablePrecomp",mit2.precomp,'Vm')
        sim.setupStim(network, ((odorA,odorB),), avgnum=avgnum)
        ## widely different resting potentials of mit0 and mit1
        if VARY_MITS_RMP:
            tweak_field('/mitrals_0/##[TYPE=Compartment]', 'Em', '-58e-3')
            tweak_field('/mitrals_1/##[TYPE=Compartment]', 'Em', '-70e-3')
        if "SAVEFORMOVIE" in sys.argv:
            SAVEFORMOVIE = True
            mitTables = setupMitralTables(network,BINNED)
        else: SAVEFORMOVIE = False
        mitral_responses,mitral_responses_binned = sim.run(network,BINNED)

        ## Save data for movie
        if SAVEFORMOVIE:
            moviedatafile = 'movie_data_netseed'+netseedstr+'_stimseed'+rateseedstr+'_directed'+str(frac_directed)
            timevec = arange(0.0,RUNTIME+1e-12,PLOTDT)
            ## number of colours in mitcolours must match number of mitral cells
            ## Each colour-entry below is a tuple of (baseline/initial colour, spiking/peak colour, colourmap)
            ## Each colour is a tuple of (r,g,b,a)
            ## Set in Moogli's config file, whether to change color, and/or alpha, or use colourmap.
            mitcolours = [ 
                ((0.3,0,0,0.3),(1,0,0,1),'jet'), ((0.3,0,0.3,0.3),(1,0,1,1),'jet'),
                ((0,0,0.3,0.3),(0,0,1,1),'jet'), ((0,0,0.3,0.3),(0,0,0.3,0.3),'jet'),
                ((0,0.3,0,0.3),(0,1,0,1),'jet'), ((0,0.3,0,0.3),(0,0.3,0,0.3),'jet') ]                
            dataTables = exportTables(network, \
                NO_PGS, NO_SINGLES, NO_JOINTS, NO_MULTIS, BINNED, \
                interneurons_args,mitcolours) ## args tells: take only interneurons connected to mits 0,1
            mitDataTables = exportMitralTables(mitTables,mitcolours,BINNED)
            ### Save colour info -- not implemented yet
            #f = open(moviedatafile+'_colours.pickle','w')
            #pickle.dump( colourTable, f )
            #f.close()
            ## Save movie data
            if BINNED: moviedatafile += '_mitspiketimes.pickle'
            else: moviedatafile += '_mitVms.pickle'
            f = open(moviedatafile,'w')
            pickle.dump( {'projections':network.projectionDict,\
                'sim_data':(timevec,dataTables,mitDataTables)}, f )
            f.close()
            print "Saved output file",moviedatafile
            sys.exit(0)
            
        ## plot and display
        if not extras_spikes_only:
            timevec = arange(0.0,RUNTIME+1e-12,PLOTDT)
            plot_extras(timevec, tables, NO_PGS, NO_SINGLES, NO_JOINTS, NO_MULTIS)
        else:
            deltabin = RESPIRATION/NUMBINS
            ## Only the last respiration cycle
            timevec = arange(SETTLETIME+(NUM_RESPS-1)*RESPIRATION+deltabin/2,RUNTIME,deltabin)
            plot_extras_spikes(timevec, tables, NO_PGS, NO_SINGLES, NO_JOINTS,\
                NO_MULTIS, NUMBINS, RUNTIME, SETTLETIME)
        figure()
        title('Glomerulus 0')
        if BINNED:
            deltabin = RESPIRATION/NUMBINS
            # Take only the last respiration cycle
            timevec = arange(SETTLETIME+(NUM_RESPS-1)*RESPIRATION+deltabin/2,RUNTIME,deltabin)
            mitral_responses = mitral_responses_binned
        else:
            timevec = arange(0.0,RUNTIME+1e-12,PLOTDT)
        plot(timevec,mitral_responses[0],color=(0.0,1.0,0.0))
        plot(timevec,mitral_responses[1],color=(0.0,1.0,0.5))
        ## plot soma; and precompartment of mit2 that inhibits mit0.
        #figure()
        #title('mitral 2')
        #plot(timevec,mitral_responses[2],color=(1,0,0))
        #plot(timevec,mit2._vmTablePrecomp,color=(0,0,0))
        show()

    #### if multiple processes are called, average over odor morphs
    else:
        ## construct the results filename
        today = datetime.date.today()
        if NO_SINGLES: singles_str = '_NOSINGLES'
        else: singles_str = '_SINGLES'
        if NO_JOINTS: joints_str = '_NOJOINTS'
        else: joints_str = '_JOINTS'
        if NO_PGS: pgs_str = '_NOPGS'
        else: pgs_str = '_PGS'
        if NO_LATERAL: lat_str = '_NOLAT'
        else: lat_str = '_LAT'
        if VARY_MITS_RMP: varmitstr = '_VARMIT'
        else: varmitstr = '_NOVARMIT'
        ## stable enough that time tags are not needed
        now =  ''#datetime.datetime.now().strftime("%Y_%m_%d_%H_%M")+'_'
        outfilename = '../results/odor_morphs/'+now+'odormorph'+\
            '_netseed'+netseedstr+'_stimseed'+rateseedstr
        if NONLINEAR_ORNS: outfilename += '_NL'+NONLINEAR_TYPE
        outfilename += singles_str+joints_str+pgs_str+lat_str+varmitstr+\
            '_numgloms'+str(NUM_GLOMS)
        if DIRECTED: outfilename += '_directed'+str(FRAC_DIRECTED)
        outfilename += '.pickle'
        
        ## if NOSHOW, then check if resultfile exists, proceed only if non-existent.
        if 'NOSHOW' in sys.argv:
            NOSHOW = True
            ## If NOSHOW, then automatic mode, hence don't overwrite resultfile, if exists beforehand.
            if os.path.exists(outfilename):
                ## activdep_inhibition_repeats.py searches for Wrote in first word,
                ## and filename in second word. so output that even if not simulating.
                if mpirank==boss:
                    for procnum in range(1,mpisize):
                        mpicomm.recv(source=procnum,tag=10)
                    print "ExistsSoNotWrote",outfilename
                else:
                    mpicomm.send('done',dest=boss,tag=10)
                sys.exit()
        else: NOSHOW = False

        if mpirank == boss:
            #### collate at boss process
            mitral_responses_list = []
            mitral_responses_binned_list = []
            numavgs = (mpisize-1)/numodors
            for avgnum in range(numavgs):
                response_odorset = []
                response_odorset_binned = []
                for odornum in range(numodors):
                    procnum = avgnum*numodors + odornum + 1
                    print 'waiting for process '+str(procnum)+'.'
                    #### you get a numpy array of rows=NUM_GLOMS*MIT_SISTERS and cols=NUMBINS
                    #### mitral responses has spike times, mitral_responses_binned has binned firing rates
                    mitral_responses,mitral_responses_binned = mpicomm.recv(source=procnum, tag=0)
                    response_odorset.append( mitral_responses )
                    response_odorset_binned.append( mitral_responses_binned )
                mitral_responses_list.append(response_odorset)
                mitral_responses_binned_list.append(response_odorset_binned)
            
            ## write results to a file
            f = open(outfilename,'w')
            pickle.dump((mitral_responses_list,mitral_responses_binned_list), f)
            f.close()
            print "Wrote", outfilename
            
            if not NOSHOW:
                plot_morphs(outfilename)
                show()

        else:
            #### run the slave processes
            sim =  odorResponse()
            avgnum = (mpirank-1)/numodors
            odorset = inputList[(mpirank-1)%numodors]
            odorA, odorB = odorset
            ## If CLUB_MITRAL=False, then extra exc from mitral sisters
            ## (to certain connected granules as proxy for unmodelled sisters) does NOT get used.
            ## Instead, here I connect extra baseline excitation to ALL granules if odor is non-zero.
            if not CLUB_MITRALS and not (odorA==0.0 and odorB==0.0):
                granfilebase += '_extra'
            ## includeProjections gets used only if ONLY_TWO_MITS is True:
            ## Keep below projections to 'second order cells'
            ## i.e. to cells (granules) connected to mits0&1.
            ## The connections between second order cell
            ## and mits0&1 are automatically retained of course.
            ## 'PG' includes 'ORN_PG', 'PG_mitral', 'mitral_PG' and 'SA_PG'
            includeProjections = ['PG','granule_baseline']
            tweaks = build_tweaks(CLUB_MITRALS, NO_SPINE_INH, NO_SINGLES,\
                NO_JOINTS, NO_MULTIS, NO_PGS, ONLY_TWO_MITS,\
                includeProjections=includeProjections, nolateral=NO_LATERAL)
            ## unique str = 'morphs_', etc so that temp files of morphs and pulses etc do not overlap
            network = OBNetwork(OBNet_file, synchan_activation_correction, tweaks,\
                mpirank, uniquestr, granfilebase, spiketable=True)
            ## widely different resting potentials of mit0 and mit1
            if VARY_MITS_RMP:
                tweak_field('/mitrals_0/##[TYPE=Compartment]', 'Em', '-58e-3')
                tweak_field('/mitrals_1/##[TYPE=Compartment]', 'Em', '-70e-3')
            #printNetTree() # from moose_utils.py

            sim.setupStim(network, (odorset,), avgnum)
            mitral_responses_both = sim.run(network, binned=True)
            mpicomm.send( mitral_responses_both, dest=boss, tag=0 )
            print 'sent from process',mpirank

Loading data, please wait...