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 -*-

########## Need to run:
## ~/Python2.6.4/bin/python2.6 CellTest_synapses.py <cellname> <synapsename> \
##      <Iclamp/Vclamp> <numsyns_tobeplotted> <timeforeachsynapsesim> \
##      <simultaneous|staggered> <compartment_name>
## the option <simultaneous|stagerred> only needs to be given if numsyns_tobeplotted is > 1.
## synapsefactor would be say GABA_factor or AMPA_factor
## (used in the network -- see networkConstants.py), typically 1.0, etc.
## Presently only event based synapses are supported.
## Typical usage: python2.6 CellTest_synapses.py mitral ORN_mitral <Vclamp|Iclamp> 1 500e-3 simultaneous
## Eg1: python2.6 CellTest_synapses.py mitral granule_mitral <Vclamp|Iclamp> 1 500e-3 simultaneous Seg0_prim_dend_1_16
## Eg2: python2.6 CellTest_synapses.py mitral ORN_mitral Iclamp 1 500e-3 simultaneous Seg0_glom_81_102
## Eg3: python2.6 CellTest_synapses.py granule mitral_granule Iclamp 100 6e-3 staggered
## Eg4: python2.6 CellTest_synapses.py mitral granule_mitral Iclamp 50 100e-3 staggered

import os
import sys
import math

sys.path.extend(["..","../channels","../synapses","../neuroml","../simulations","../networks"])

from moose.utils import *
from synapseConstants import *
from simset_inhibition import * # has SIMDT and PLOTDT but SETTLETIME is overridden below
SETTLETIME = 500e-3 # give a large settletime for measuring post synaptic potentials

from load_channels import *
from load_synapses import *
from networkConstants import * # for GABA_factor, AMPA_factor and NMDA_factor

from pylab import * # part of matplotlib that depends on numpy but not scipy
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes # for inset
from mpl_toolkits.axes_grid1.inset_locator import mark_inset

from data_utils import *

from moose.neuroml.MorphML import *

SATURATING_MITGRAN = False
if SATURATING_MITGRAN:
    mitgran_ampa = 'mitral_granule_saturatingAMPA'
    mitgran_nmda = 'mitral_granule_saturatingNMDA'
else:
    mitgran_ampa = 'mitral_granule_AMPA'
    mitgran_nmda = 'mitral_granule_NMDA'

#### inject a current into the tuft of mitral.
#### This to check effect of IPSP in prim dend vs sec dend.
#injtuft = 200e-12#1e-10 # A
injtuft = 0.0 # A
#injsoma = 400e-12 # A
injsoma = 0.0 # A

## nerve_shock is for injecting a sharp current into the mitral tuft after SETTLETIME
nerve_shock = False
num_tuft_comps_shock = 53 # number of tuft compartents to inject current into
tuft_shock_inject = 1e-12 # this current is injected into each of the num_tuft_comps_shock above

## set 0 mV for IPSCs, -65mV for EPSCs
VCLAMP_V = -65e-3#0e-3 # mV

class CellTest:

    def __init__(self, cellname, pseudosynname, libsynname,
        synfactor, Vclamp, numsyns, syntime, stagger, comp_name=None):
        load_channels()
        ## synchan_activation_correction is for graded synapses -- not useful for this script
        load_synapses(synchan_activation_correction)
        MML = MorphML({'temperature':CELSIUS})
        if cellname == 'PG':
            #filename = 'PG_aditya2010_neuroML_L1_L2_L3.xml'
            #filename = 'PG_aditya2012_neuroML_L1_L2_L3.xml'
            #filename = 'PG_aditya2013_neuroML_L1_L2_L3.xml'
            #cellname = 'PG_LTS'
            ## Choose one the two below
            ## set the ORN->PG and mitral->PG to 0.45 vs 1.25 nS for plateauing vs LTS
            #filename = 'PG_aditya2010unified_neuroML_L1_L2_L3.xml' # plateauing i.e. non-LTS
            filename = 'PG_aditya2013unified_neuroML_L1_L2_L3.xml' # LTS
            self.somaName = 'soma_0'
        elif cellname == 'mitral':
            #filename = 'mitral_bbmit1993davison_neuroML_L1_L2_L3_mod.xml'
            #filename = 'mitral_bbmit1993davison_neuroML_L1_L2_L3.xml'
            filename = 'mitral_bbmit1993davison_neuroML_L1_L2_L3_mod_withspikeinit.xml'
            #filename = 'mitral_bbmit1993davison_neuroML_TEST_L1_L2_L3.xml'
            #filename = 'mitral_bbmit1993davisonMS_neuroML_L1_L2_L3.xml'
            #filename = 'mitral_migliore_etal_2007.xml'
            self.somaName = 'Seg0_soma_0'
        elif cellname == 'granule':
            filename = 'granule_granadityaMS2007_neuroML_L1_L2_L3.xml'
            self.somaName = 'soma_0'
        self.cellName = cellname
        self.synName = libsynname
        self.pseudoSynName = pseudosynname
        self.synFactor = synfactor
        self.Vclamp = Vclamp
        self.numSyns = numsyns
        self.synTime = syntime
        self.stagger = stagger
        self.compName = comp_name
        if self.stagger:
            self.synRUNTIME = self.synTime*(self.numSyns)+2*SETTLETIME
        else:
            self.synRUNTIME = self.synTime+2*SETTLETIME

        self.cellsDict = MML.readMorphMLFromFile(filename,{})
        self.context = moose.PyMooseBase.getContext()
        self.cell = moose.Cell( self.context.deepCopy(\
            self.context.pathToId('/library/'+cellname),self.context.pathToId('/'),cellname) )
        self.soma = moose.Compartment('/'+self.cellName+'/'+self.somaName)
        self.somaVm = setupTable('soma_vm', self.soma,'Vm')
        if cellname=='mitral':
            ## inject a current into tuft of mitral
            tuftcomp = moose.Compartment('/'+self.cellName+'/Seg0_glom_81_102')
            tuftcomp.inject = injtuft
            ## inject a current into soma of mitral
            self.soma.inject = injsoma
            ## if nerve shock, inject a sharp current into tuft compartments
            if nerve_shock:
                ## We want a single sharp current injection at SETTLETIME, so cannot use soma.inject
                ## segment info - see MorphML_reader.py
                ## Only connect to 25 tuft compartments
                for seginfo in (self.cellsDict[self.cellName].values())[0:num_tuft_comps_shock]:
                    ## seginfo[5] is a list of potential synapses at this segment
                    if 'ORN_mitral' in seginfo[5]:
                        ## wrap this segment, seginfo[0] is segment name
                        tuftcomp = moose.Compartment('/'+self.cellName+'/'+seginfo[0])
                    ## compartment, name_extn, start_time, duration, current (all SI)
                    setup_iclamp(tuftcomp, '_nerveshock',\
                        SETTLETIME, 1e-3, tuft_shock_inject)

        self.spikeTableList = []
        self.attachSpikeTables()
        #printCellTree(self.cell)

        if cellname=='mitral':
            ## monitor Vm at the base of the tuft
            tuftbase_seg = moose.Compartment('/'+self.cellName+\
                #'/Seg0_tuftden_19_23') # for migliore and shepherd 2007 cell
                '/Seg0_prim_dend_5_20') # tuft base
            self.tuftBaseTable = setupTable('mitdendTable',tuftbase_seg,'Vm')

            ## monitor Vm at secondary dendrite
            dend_seg = moose.Compartment('/'+self.cellName+\
                #'/Seg0_sec_dendp4_0_254') # at 43.86 microns from soma
                #'/Seg0_sec_dendd3_0_204') # at 190.56 microns from soma
                '/Seg0_sec_dendd4_2_269') # at 1004.47 microns from soma
            self.secDendTable = setupTable('mitdendTable',dend_seg,'Vm')

            ## monitor Vm in the tuft compartment
            ## same as current injection above
            self.tuftTable = setupTable('mitTuftTable',tuftcomp,'Vm')

        ### testing - channels
        ## setup only one of the Tables below
        #self.var = setupTable('PG_mitral', moose.SynChan('/mitral/Seg0_glom_0_21/PG_mitral'), 'Gk')
        #self.var = setupTable('mitral_granule_NMDA', \
        #    moose.SynChan('/'+self.cellName+'/'+self.compName+'/mitral_granule_NMDA'), 'Gk')
        #self.var = setupTable('PG_mitral', moose.SpikeGen('/PG/soma_0/PG_mitral_spikegen'), 'state')
        #self.var = setupTable('soma_vm', moose.Compartment('/PG/soma_0'), 'Vm')
        #self.var = setupTable('mitral_glom', moose.Compartment('/mitral/Seg0_glom_0_21'), 'Vm')
        #self.var = setupTable('mitral_secdend', moose.Compartment('/mitral/Seg0_sec_dendd2_7_200'), 'Vm')
        #self.var = setupTable('mitral_secdend', moose.Compartment('/mitral/Seg0_sec_dendd2_0_165'), 'Vm')
        #self.var = setupTable('mitral_glom', moose.Compartment('/mitral/Seg0_glom_81_102'), 'Vm')
        #self.var = setupTable('soma_vm', moose.Compartment('/mitral/Seg0_soma_0'), 'Vm')
        #self.var = setupTable('soma_ca', moose.CaConc(self.soma.path+'/Ca_mit_conc'), 'Ca')
        
    def attachSpikeTables(self):
        ## If user did not supply a compartment_name, make a list of synapses
        ## connected to allowed compartments in the cell's xml file
        if self.compName is None:
            synapse_list = []
            synapse_list2 = []
            ## segment info - see MorphML_reader.py
            for seginfo in self.cellsDict[self.cellName].values():
                ## seginfo[5] is a list of potential synapses at this segment
                if self.pseudoSynName in seginfo[5]:
                    ## wrap this segment, seginfo[0] is segment name
                    compartment = moose.Compartment('/'+self.cellName+'/'+seginfo[0])
                    synapse = connectSynapse(self.context, compartment, self.synName, self.synFactor)
                    ## NMDA is also created if AMPA mit->gran is made.
                    if self.synName == mitgran_ampa:
                        synapse2 = connectSynapse(self.context, compartment, mitgran_nmda, NMDA_factor)
                        synapse_list2.append(synapse2)
                    synapse_list.append(synapse)
                    self.compName = seginfo[0] # the last synapse's comp becomes self.compName
        ## If user supplied a compartment_name,
        ## just have a synapse connected to the specified compartment
        else:
            #printCellTree(self.cell)
            comp_path = '/'+self.cellName+'/'+self.compName
            if self.context.exists(comp_path): compartment = moose.Compartment(comp_path)
            else:
                print "Error: non-existent compartment", comp_path
                sys.exit(1)
            synapse_list = [ connectSynapse(self.context, compartment, self.synName, self.synFactor) ]
            if self.synName == mitgran_ampa:
                synapse_list2 = [ connectSynapse(self.context, compartment, mitgran_nmda, NMDA_factor) ]
        for i in range(self.numSyns):
            ## Bad practice - This should be part of NetworkML
            ## - no random numbers in simulations, only in generators.
            randomsynindex = int(uniform(0,len(synapse_list)))
            synapse = synapse_list[randomsynindex]
            print "Connecting ", synapse.path
            ## put in i also, because randomly a synapse might get connected twice.
            spiketable = moose.TimeTable(synapse.path+'/tt_'+str(i))
            #### SynChan's synapse MsgDest takes time as its argument.
            ## Thus spiketable should contain a list of spike times.
            spiketable.connect("event", synapse,"synapse")
            if self.synName == mitgran_ampa:
                synapse2 = synapse_list2[randomsynindex]
                spiketable.connect("event", synapse2,"synapse")
            self.spikeTableList.append(spiketable)
            ## Presently only method 4 i.e. loading from file is supported for TimeTable.
            ## Hence the need to write a single value in a file.
            fn = 'temp_spikefile.txt'
            f = open(fn,'w')
            if self.stagger:
                f.write(str(i*self.synTime+SETTLETIME))
            else:
                f.write(str(SETTLETIME))
            f.close()
            spiketable.filename = fn
            os.remove(fn)
            if self.Vclamp:
                    self.synITable = setupTable('synITable',synapse,'Ik')
        if self.Vclamp:
            #### I block Ca and K channels in the cell
            #### We don't want these channels to be active when cell is held at 0mV; Na inactivates.
            blockChannels(self.cell, ['K','Ca']) # in moose_utils.py
            ## PID gain: I think for Davison 4-comp-mitral/granule: 0.5e-5 # optimal gain
            ## too high 0.5e-4 drives it to oscillate at high frequency,
            ## too low 0.5e-6 makes it have an initial overshoot (due to Na channels?)
            ## But for BBmit1993, gain of 1e-6 is optimal
            self.PIDITable = setup_vclamp(self.soma, '_somavclamp', 0,
                self.synRUNTIME, VCLAMP_V, gain=1e-6)
            print "Connected Vclamp to soma at",VCLAMP_V,"mV."

    def run(self):
        #printCellTree(self.cell)
        resetSim(self.context, SIMDT, PLOTDT)
        self.context.step(self.synRUNTIME) # each synapse is given 500ms!

        ## Paper figure 2 for granule cell electrophysiology
        fig = figure(figsize=(columnwidth/2.0,linfig_height/2.0),dpi=300,facecolor='w') # 'none' is transparent
        ax = fig.add_subplot(111)
        starttime = SETTLETIME*1000-10
        timevec = linspace(0.0,self.synRUNTIME,len(self.somaVm))*1000 - starttime
        self.somaVm = array(self.somaVm)*1000
        plot(timevec, self.somaVm,'-k', label='soma', linewidth=plot_linewidth, linestyle='-') # ms and mV
        xmin,xmax,ymin,ymax = beautify_plot(ax,x0min=True,y0min=False,drawxaxis=True,drawyaxis=True,\
            xticks=arange(0,self.synRUNTIME*1000-starttime,200),yticks=[])
        ax.set_xlim(0,800)
        ax.set_yticks([ymin,0,ymax])
        axes_labels(ax,'time (ms)','Vm(mV)',fontsize=label_fontsize,xpad=1,ypad=-5.5) # sets fontsize for ticks and labels
        #fig.text(0.04,0.7,'Vm (mV)', fontsize=label_fontsize, rotation='vertical', transform = fig.transFigure)
        
        if self.cellName == 'granule':
            ## Paper figure 2 for granule cell electrophysiology
            ## zoomed inset showing synaptic integration
            ## bbox_to_anchor can be instance of BboxBase, (x, y, width, height of the bbox), or (x, y) with width=height=0
            #axins = zoomed_inset_axes(ax, zoom=10, bbox_to_anchor=(0.43,1.0), bbox_transform=ax.transAxes)
            axins = inset_axes(ax, width=0.4, height=0.6, bbox_to_anchor=(0.5,1.2), bbox_transform=ax.transAxes)
            ## set the linewidth of the axes (spines)
            for s in axins.spines.values():
                 s.set_linewidth(axes_linewidth)
            ## Important: you need to plot the original plot again!!!
            plot(timevec, self.somaVm,'-k', label='soma', linewidth=plot_linewidth, linestyle='-') # ms and mV        
            zoom_xmin,zoom_xmax = 210, 250
            axins.set_xlim(zoom_xmin,zoom_xmax)
            ## take out the Vm data that is in the zoom region
            ## timevec may not have zoom_xmin exactly, so use index for value just below it, similar for zoom_xmax
            zoom_data = self.somaVm[where(timevec<zoom_xmin)[0][-1]:where(timevec>zoom_xmax)[0][0]]
            axins.set_ylim(min(zoom_data),max(zoom_data))
            ### autoscale will not work -- as the full data, not a subset is plotted !!
            #axins.autoscale(enable=True, axis='y', tight=True)
            #axins.autoscale_view(tight=True, scaley=True)
            #ymin, ymax = axins.get_yaxis().get_view_interval()
            #axins.set_ylim(ymin,ymax)
            axins.set_xticks([])
            axins.set_yticks([])
            mark_inset(ax, axins, loc1=3, loc2=4, fc="none", ec="0.5", linewidth=axes_linewidth)
            #add_scalebar(axins,matchx=False,matchy=False,hidex=False,hidey=False,\
            #    sizex=20,labelx='20 ms',sizey=1,labely='  1 mV',\
            #    bbox_to_anchor=[1.3,0.0],bbox_transform=axins.transAxes)
            axes_labels(axins,' {:3.0f} ms'.format(zoom_xmax-zoom_xmin),\
                '     {:3.2f} mV'.format(max(zoom_data)-min(zoom_data)),\
                xpad=-6,ypad=-6,fontsize=label_fontsize-2)
        elif self.cellName=='mitral':
            #plot(timevec, array(self.var)*1000,',-b', label='Sec dend Vm')
            plot(timevec, array(self.tuftBaseTable)*1000,'-b',\
                label='tuft base', linewidth=plot_linewidth, linestyle='--')
            plot(timevec, array(self.tuftTable)*1000,'-g',\
                label='tuft', linewidth=plot_linewidth, linestyle='-.')
            biglegend()
            
            ################ Paper figure 2
            if nerve_shock:
                fig2 = figure(figsize=(columnwidth/2.0,linfig_height/2.0),dpi=300,facecolor='w') # 'none' is transparent
                ax = fig2.add_subplot(111)
                plot(timevec, self.somaVm,'-k', label='soma', linewidth=linewidth, linestyle='-') # ms and mV
                plot(timevec, array(self.tuftBaseTable)*1000,'-b',\
                    label='tuft base', linewidth=linewidth, linestyle='--')
                beautify_plot(ax,drawxaxis=False,drawyaxis=False,xticks=[],yticks=[])
                add_scalebar(ax,matchx=False,matchy=False,hidex=True,hidey=True,\
                    sizex=3,labelx='3 ms',sizey=20,labely='20 mV',\
                    bbox_to_anchor=[0.3,0.6],bbox_transform=ax.transAxes)
                ax.set_xlim(SETTLETIME*1000-10,(SETTLETIME+15e-3)*1000) #ms
                ax.set_ylim(-70,50)
                fig2.tight_layout()
                fig2.savefig('../figures/connectivity/cells/mitral_spikeinit.png',dpi=fig2.dpi)
                fig2.savefig('../figures/connectivity/cells/mitral_spikeinit.svg',dpi=fig2.dpi)

        fig.tight_layout()
        savefig('../figures/connectivity/cells/'+self.pseudoSynName+'.png',dpi=fig.dpi)
        savefig('../figures/connectivity/cells/'+self.pseudoSynName+'.svg',dpi=fig.dpi)

        #figure(facecolor='w')
        #plot(timevec, self.var,',-b', label='Gk')
        if self.Vclamp:
            figure(facecolor='w')
            plot(timevec, self.synITable,'g-,')
            title("synapse current")
            xlabel('time (ms)',fontsize='large')
            ylabel('I (A)',fontsize='large')
            figure(facecolor='w')
            plot(timevec, self.PIDITable,'g-,')
            title("V clamp current")
            xlabel('time (ms)',fontsize='large')
            ylabel('I (A)',fontsize='large')
        show()
        
if __name__ == "__main__":
    seed([122.0])
    numsyns = int(sys.argv[4])
    stagger = False
    if numsyns > 1 and sys.argv[6] == 'staggered':
        stagger = True
    if sys.argv[3]=='Vclamp': Vclamp = True
    elif sys.argv[3]=='Iclamp': Vclamp = False
    else:
        print "You must specify Vclamp or Iclamp. Using Iclamp by default."
        Vclamp = False
    if sys.argv[2]=='granule_mitral':
        libsynname = 'granule_mitral_GABA'
        synfactor = GABA_factor*strongsynfactorinh
    elif sys.argv[2]=='mitral_granule':
        libsynname = mitgran_ampa
        synfactor = AMPA_factor
    else:
        libsynname = sys.argv[2]
        synfactor = 1.0
    if len(sys.argv)>=8: comp_name = sys.argv[7]
    else: comp_name = None
    cell = CellTest(cellname=sys.argv[1],pseudosynname=sys.argv[2],libsynname=libsynname,
        synfactor=synfactor,Vclamp=Vclamp, numsyns=numsyns,
        syntime=float(sys.argv[5]),stagger=stagger, comp_name=comp_name)
    cell.run()

Loading data, please wait...