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

########## THIS FITTING PROGRAM IS MEANT TO ROUGHLY FOLLOW PRIYANKA'S ANALYSIS
########## This variant does not use an air kernel, rather a constant air rate.
## USAGE1: python2.6 fit_scaledpulses.py <../results/odor_pulses/scaledpulses_ ... .pickle> <stimseed> [SAVEFIG]

from scipy import optimize
from scipy import special
import scipy.interpolate
from scipy import signal
import scipy.io
from scipy.stats import linregress
from pylab import *
import pickle
import sys
import math
import copy as cp

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

from stimuliConstants import * # has SETTLETIME, SCALED_RUNTIME
from networkConstants import * # has central_glom
from data_utils import * # has axes_labels()
from analysis_utils import * # has load_morph...(), predict_respresp() and set_min_err()

## index of the reference response
ref_response_scalenum = 3
## if peak detect is True, peak scaling is plotted,
## rather than the scaling wrt reference response above
peak_detect = False#True

## rebin the spikes as below, irrespective of previous binning
pulserebindt = 50e-3#fitting_dt # 50ms as Priyanka uses
pulserebins = int(SCALED_RUNTIME/pulserebindt)
#bin_width_time = 2*pulserebindt ## overlapping bins with this bin width
pulsetlist = arange(pulserebindt/2.0,SCALED_RUNTIME,pulserebindt)

NOISE_ANALYSIS = False

## number of the mitral to be fitted.
fitted_mitrals = [2*central_glom,2*central_glom+1]

class fit_plot_scaledpulses():
    def __init__(self,filename,stimseed,scaledWidth=scaledWidth):
        self.filename = filename
        self.stimseed = stimseed
        self.scaledWidth = scaledWidth

    def fit_pulses(self,dirextn='',_noshow=True,_savefig=False, _test=False):
        if _test:
            ########## load in the stimuli
            ## scaledPulseList[glomnum][scalenum][odornum][binnum]
            scaledPulsesList,genORNkernels \
                = read_scaledpulses_stimuli(self.stimseed,self.scaledWidth)
            ## decimate pulseList by taking every convolutiondt/FIRINGFILLDT bins
            ## and decimate ORNkernels by taking every kerneldt/FIRINGFILLDT=50/1=50 bins
            ## pulserebinList[scalenum,binnum] numpy array
            pulserebinList,linORNkernelA,linORNkernelB = \
                decimate(scaledPulsesList[0,:,0],pulserebindt,genORNkernels,kerneldt,kernel_size)       
            numtrials = 1
        else:
            ######### load in the mitral pulse responses
            #### mitral_responses_list[avgnum][scalenum][mitralnum][spikenum]
            #### mitral_responses_binned_list[avgnum][scalenum][mitralnum][binnum]
            mitral_responses_list, mitral_responses_binned_list = read_pulsefile(self.filename)

            if NOISE_ANALYSIS:
                ##---------- Print/plot for best bin size, but this bin size is not used presently
                ## From Neural Computation 19, 1503–1527 (2007) Shimazaki & Shinomoto, pg 6
                DeltaList = arange(25.0e-3,500.0e-3,25.0e-3)
                SSvals = []
                SSvals_extended = []
                for Delta in DeltaList:
                    numtrials,mitral_responses_mean1,mitral_responses_std1 = \
                            rebin_mean(mitral_responses_list,Delta,SCALED_RUNTIME)
                    fitted_mitral = 0
                    firingbinsmeanList1 = mitral_responses_mean1[:,fitted_mitral]
                    firingbinserrList1 = mitral_responses_std1[:,fitted_mitral]
                    ## each val in firingbinsmeanList[scalenum][binnum] is k_i/(n*Delta)
                    kbar = mean(firingbinsmeanList1)*numtrials*Delta
                    kvariance = (std(firingbinsmeanList1)*numtrials*Delta)**2
                    shimazaki_shinomoto_val = (2*kbar-kvariance)/((numtrials*Delta)**2)
                    SSvals.append(shimazaki_shinomoto_val)
                    print 'Delta t =',Delta,"shimazaki_shinomoto_val =",shimazaki_shinomoto_val
                    ## Cost function i.e. shimazaki_shinomoto_val for different # of trials
                    Priyanka_numtrials = 12
                    ssval_extended = (1.0/Priyanka_numtrials - 1.0/numtrials)*kbar/numtrials/Delta**2 + \
                                        shimazaki_shinomoto_val
                    SSvals_extended.append(ssval_extended)
                if not _noshow:
                    fig = figure()
                    ax = fig.add_subplot(111)
                    ax.plot(DeltaList,SSvals,label='#='+str(numtrials))
                    ax.plot(DeltaList,SSvals_extended,label='#='+str(Priyanka_numtrials))
                    ax.legend()
                    axes_labels(ax,'Delta t (s)','SS Cost function (Hz^2)',fontsize=20)

            ##-------------------------- rebin the responses and pulses ------------------------------
            ## rebin sim responses to pulserebindt=50ms, then take mean
            numtrials,mitral_responses_mean,mitral_responses_std = \
                    rebin_mean(mitral_responses_list,pulserebindt,SCALED_RUNTIME)

        ## full fitting data for both mitrals
        fits_2mits = []
        peak_scales_2mits = []
        for mit_i,fitted_mitral in enumerate(fitted_mitrals):
            if _test:
                firingbinsmeanList = pulserebinList
                firingbinsmeanList += uniform(-0,0,shape(firingbinsmeanList))
                firingbinserrList = zeros(shape(firingbinsmeanList))
            else:
                ## take the odor responses of the mitral to be fitted
                firingbinsmeanList = mitral_responses_mean[:,fitted_mitral]
                ## The model predicts the individual response not the mean.
                ## Hence below fitting uses standard deviation, not standard error of the mean.
                firingbinserrList = mitral_responses_std[:,fitted_mitral]
            
            starti = int(PULSE_START/pulserebindt)
            endi = int((PULSE_START+scaledWidth+kernel_time)/pulserebindt)
            air_bgnd = firingbinsmeanList[0]
            air_bgnd_relevant = firingbinsmeanList[0][starti:endi]

            ##---------------------------- fit scaled pulse responses ------------------------------------------
                            
            ## define the reference response / scaling
            ref_scale = scaledList[ref_response_scalenum]
            ref_response = firingbinsmeanList[ref_response_scalenum][starti:endi]-air_bgnd_relevant
            fits = []
            peak_scales = []
            for scalenum in [1,2,3,4,5]: ## conc scaled pulses
                scale = scaledList[scalenum]
                response = firingbinsmeanList[scalenum][starti:endi]-air_bgnd_relevant
                ## http://www.jerrydallal.com/LHSP/slrout.htm for defn of std error of the estimate
                ## SEE is std error of data about the regression line
                slope, intercept, r_value, p_value, see = linregress(ref_response,response)
                ## SEE is called \hat{\sigma}_\epsilon i.e. sqrt(MSE) here:
                ## http://en.wikipedia.org/wiki/Regression_analysis
                ## formula for std error of slope is also from above Wikipedia article
                se_slope = see/std(ref_response)
                avgfrate = sum(firingbinsmeanList[scalenum][starti:endi])/float(endi-starti)
                fits.append((scale,slope,intercept,r_value,p_value,see,se_slope,avgfrate))
                ## peak scaling
                peak_scales.append(max(response))
            peak_scales =  array(peak_scales)#/peak_scales[ref_response_scalenum-1]*scaledList[ref_response_scalenum]
            peak_scales_2mits.append(peak_scales)
            fits_2mits.append(fits)

            ##---------------------------- plot scaled pulse responses -----------------------------------------

            if not _noshow:
                if peak_detect: print "BEWARE. Using peak scaling"
                else: print "BEWARE. Using fitted scaling"
                ############################### plot the responses and the fits
                fig = figure(figsize=(columnwidth,linfig_height/2),dpi=300,facecolor='w') # 'none' is transparent
                ## conc scaled pulses, leave the 0th pulse which is air_bgnd
                ax = plt.subplot2grid((1,3),(0,0),rowspan=1)
                for scaleiter,scale in enumerate(scaledList[1:]): ## conc scaled pulses
                    sister_ratio = (fitted_mitral%MIT_SISTERS)/float(MIT_SISTERS)
                    scaledpulsetime = array(pulsetlist[starti:endi]) - pulsetlist[starti] # start from t=0
                    ################### Plot the simulated responses
                    ## smooth the simulated response
                    ## windowsize=5 and SD=0.65 are defaults from matlab's smoothts() for gaussian smoothing
                    Gwindow = signal.gaussian(5,0.65)
                    ## help from http://www.scipy.org/Cookbook/SignalSmooth
                    simresponse = convolve(Gwindow/Gwindow.sum(),\
                        firingbinsmeanList[scaleiter+1]-air_bgnd,mode='same')
                    ## ditch the smoothing above for scaled pulses
                    simresponse = firingbinsmeanList[scaleiter+1][starti:endi]-air_bgnd_relevant
                    ## numpy array, hence adds element by element
                    scale_color = ['r','b','g','m','c'][scaleiter]
                    scale_label = ['1/3x','2/3x','1x','2x','5x'][scaleiter]+\
                            ' %1.2f'%(fits_2mits[-1][scaleiter][3]) # R-value = corr
                    print scale_label
                    fill_between(scaledpulsetime,
                        simresponse+firingbinserrList[scaleiter+1][starti:endi]/sqrt(numtrials),
                        simresponse-firingbinserrList[scaleiter+1][starti:endi]/sqrt(numtrials),
                        color=scale_color,alpha=0.3)
                    plot(scaledpulsetime,simresponse,linewidth=plot_linewidth,color=scale_color,label=scale_label)
                xmin,xmax,ymin,ymax = beautify_plot(ax,x0min=True,y0min=False,xticksposn='bottom',yticksposn='left')
                #biglegend(fontsize=label_fontsize-2, labelspacing=0., handletextpad=0.)
                ax.set_xticks([0,xmax])
                if ymin<10: ax.set_yticks([ymin,0,ymax])
                else: ax.set_yticks([ymin,ymax])
                plot([0.,scaledWidth],[ymin+2,ymin+2],linewidth=plot_linewidth*3,color='r')
                axes_labels(ax,'time (s)','rate (Hz)',adjustpos=False,xpad=0,ypad=-4)
                #ax.yaxis.set_label_coords(-0.4,1.2)
                    
                ## plot corr vs conc
                ax = plt.subplot2grid((1,3),(0,1))
                concratios,slopevsconc,_,_,_,_,se_slope,avg_frates = zip(*fits)
                ax.plot(scaledList[1:],array(fits_2mits[-1])[:,3],color='b',linewidth=linewidth,\
                    marker='o',ms=marker_size)
                beautify_plot(ax,x0min=False,y0min=False,xticksposn='bottom',yticksposn='left')
                ax.set_xlim(0,5)
                ax.set_xticks([0,1,2,5])
                ax.set_ylim(0,1)
                ax.set_yticks([0,1])
                axes_labels(ax,'conc (% SV)','corr to 1%',adjustpos=False,xpad=1,ypad=2)

                ## plot response scaling vs conc scaling
                ax = plt.subplot2grid((1,3),(0,2))
                concratios,slopevsconc,_,_,_,_,se_slope,avg_frates = zip(*fits)
                print "Average firing rates for mitral",fitted_mitral,\
                    "for different scales is",avg_frates
                print "Peak firing rates for mitral",fitted_mitral,\
                    "for different scales is",peak_scales
                ax.plot(scaledList,append([0],peak_scales),color='b',linewidth=linewidth,\
                    marker='o',ms=marker_size)
                beautify_plot(ax,x0min=False,y0min=True,xticksposn='bottom',yticksposn='left')
                ax.set_xlim(0,5)
                ax.set_xticks([0,1,2,5])
                axes_labels(ax,'conc (% SV)','peak (Hz)',adjustpos=False,xpad=1,ypad=-3)
                #ax.yaxis.set_label_coords(-0.3,1.2)

                fig.tight_layout()
                fig_clip_off(fig)
                #fig.subplots_adjust(top=0.94,left=0.1,right=0.99,hspace=0.4,wspace=0.5)

                if _savefig:
                    fig.savefig('../figures/scalelinearity_example_'+str(self.stimseed)+\
                        '_mit'+str(mit_i)+'.svg',dpi=fig.dpi)
                    fig.savefig('../figures/scalelinearity_example_'+str(self.stimseed)+\
                        '_mit'+str(mit_i)+'.png',dpi=fig.dpi)

                if NOISE_ANALYSIS:
                    ## plot the variance vs firing rate mean for each mitral
                    ## variance = mean/bintime of firng rate for Poisson process
                    fig2 = figure()
                    ax2 = fig2.add_subplot(111)
                    for scaleiter,scale in enumerate(scaledList): ## conc scaled pulses including air
                        ax2.scatter( firingbinsmeanList[scaleiter], \
                            firingbinserrList[scaleiter]**2, color='r' )
                    beautify_plot(ax2)
                    axes_labels(ax2,'mean rate (Hz)','variance (Hz^2)',fontsize=14)
                    
                    ## plot individual trials for a given response
                    fig3 = figure()
                    ax3 = fig3.add_subplot(111)
                    for trialspikelist in mitral_responses_list: 
                        plot( plotBins( trialspikelist[0][fitted_mitral],\
                            pulserebins, SCALED_RUNTIME, 0.0) )

        return fits_2mits, peak_scales_2mits
        

if __name__ == "__main__":
    NOSHOW = False
    if 'SAVEFIG' in sys.argv: SAVEFIG = True
    else: SAVEFIG = False
    if len(sys.argv) > 2:
        filename = sys.argv[1]
        stimseed = sys.argv[2]
        worker = fit_plot_scaledpulses(filename,stimseed)
        post_pulses = filename.split('odor_pulses')[1]
        dirextn = post_pulses.split('/')[0]
        print 'directory extension is',dirextn
        if 'TEST' in sys.argv: TEST=True
        else: TEST=False
        worker.fit_pulses(dirextn,NOSHOW,SAVEFIG,TEST)
        show()
    else:
        print "At least specify data file containing pickled mitral responses, and ORN frate seed."
        sys.exit(1)

Loading data, please wait...