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 BE A CLONE OF MUKUND'S AND ADIL'S MATLAB ONE
## USAGE: python2.6 fit_odor_morphs.py ../results/odor_morphs/2011-01-13_odormorph_SINGLES_JOINTS_PGS.pickle [CHISQ_HIST] [SAVEFIG]

from scipy import optimize
from scipy.special import * # has error function erf() and inverse erfinv()
from pylab import *
import pickle
import sys
import math

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

from stimuliConstants import * # has SETTLETIME, inputList and pulseList, GLOMS_ODOR, GLOMS_NIL
from networkConstants import * # has central_glom
from sim_utils import * # has rebin() to alter binsize
from analysis_utils import * # has read_morphfile() and NUM_REBINS, etc.

## use error function(x) for x>=0 (zero for x<0),
## OR use sigmoid(x) (non-zero for -ve x)
USE_ERF = False#True

iterationnum = 1

## I don't use the NUMBINS in simset_odor.py, rather I rebin() with below NUM_REBINS
## Adil used 17 bins for a 1s rat respiration cycle.
## I'm using 9 bins to get the same binwidth, else there are oscillations ~ 35 Hz gamma?
NUM_REBINS = 9#17

NUMMIX = len(inputList)
## remove the two pure odors and one pure air weights
NUMWTS = NUMMIX-3
firstrun = False#True

### numbers of mitral to be fitted.
fitted_mitral_list = [2*central_glom+0, 2*central_glom+1]

## Fit type: 'lin' : linear or 'arb' : monotonic arbitrary
## if arbitrary fit_type, weights are also free params,
## if linear fit_type, weights are not free params.
## This param is passed to fit_morphs()
fit_type = 'arb'

log81 = math.log(81)
 
def constrain0to1(x):
    try:
        return exp(x)/(1+exp(x)) # use numpy's exp
    except OverflowError as overflowerr:
        print overflowerr
        print x
        return 1.0

# define sigmoid which runs from (-0.5,0.1) to (+0.5,0.9)    
# Ideally the fitted sigmoid should be shifted by 0.5 i.e.
# exp((x-0.5)*log81)/(1+exp((x-0.5)*log81))
# This will overlap most of the linear part.
# But for fitting it doesn't matter,
# the fit routine will shift the parameters as required.
# But while plotting the internal response parameters,
# shift by 0.5 and plot -- see below
def outputsigmoid(x):
    if USE_ERF:
        if x<0: return 0
        else: return erf(x)
    else:
        try:
            return exp(x*log81)/(1+exp(x*log81)) # use numpy's exp
        except OverflowError as overflowerr:
            print overflowerr
            print x
            return 1.0
    
def inversesigmoid(x):
    if USE_ERF:
        if x<0: return x
        else: return erfinv(x)
    else:
        ## just to set initial values, value doesn't matter too much when x tends to 0
        if x>1e-200: return math.log(x/(1-x))
        else: return -5e2
        
def rectifier(x):
    x[where(x<0)[0]]=0
    return x

def chisqfunc(params, ydata, errdata, fit_type):
    RA = params[0:NUM_REBINS]
    RB = params[NUM_REBINS:2*NUM_REBINS]
    Rair = params[2*NUM_REBINS:3*NUM_REBINS]
    if fit_type == 'arb':
        #### for the weights also, we use exactly what is done by Mukund and Adil in matlab:
        #### constrain weights to be between 0 and 1
        #### sort the weights to ensure monotonicity
        inputsA = [ constrain0to1(x) for x in params[3*NUM_REBINS:(3*NUM_REBINS+NUMWTS)] ]
        ## important to put these in else along with sort(),
        ## weights saturate at 0.9 or so rather than at 1.0
        inputsA.extend([0.0,1.0]) # for pure odors
        inputsA.sort() # in place sort
        inputsA.append(0.0) # for air - keep this after sort!
        inputsB = [ constrain0to1(x) for x in params[(3*NUM_REBINS+NUMWTS):(3*NUM_REBINS+2*NUMWTS)] ]
        ## important to put these in else along with sort(),
        ## weights saturate at 0.9 or so rather than at 1.0
        inputsB.extend([0.0,1.0]) # for pure odors
        inputsB.sort(reverse=True) # weights of odor B need to be used in reverse
        inputsB.append(0.0) # for air - keep this after sort!
        #### Mukund and Adil constrained sigmoidmax > ydatamax (note exp(x)>0.)
        sigmoidmax = ydata.max() + exp(params[3*NUM_REBINS+2*NUMWTS])
    else:
        ## *-operator unpacks the list which become args of zip()
        ## zip collects the i-th elements together of all the args.
        inputsA,inputsB = zip(*inputList) # keep the last (0,0) air input        

    global iterationnum
    if iterationnum%1000==0: print 'iteration number =',iterationnum
    #if iterationnum%100==0: print inputsA, inputsB
    chisqarray = [0.0]
    for i,(inputA,inputB) in enumerate(inputList):
        CA = inputsA[i]
        CB = inputsB[i]
        if fit_type == 'arb':
            for bin in range(NUM_REBINS):
                Rmix = sigmoidmax*outputsigmoid( Rair[bin] + CA*RA[bin] + CB*RB[bin] )
                chisqarray.append( (ydata[i][bin] - Rmix)/errdata[i][bin] ) # divide by error to do chi-square fit
        else:
            Rmix = rectifier(Rair + CA*RA + CB*RB)
            chisqarray.extend( (ydata[i] - Rmix)/errdata[i] )
            
    ## not yet squared, so normalize 'chi' to sqrt of number of dof
    chisqarray = array(chisqarray) / sqrt(ydata.size-params.size)
    iterationnum += 1
    return chisqarray # misnomer -- actually individual chi array

def fit_morphs(filename, fitted_mitral, fit_type='arb', refit=True):
    ## The model predicts the individual response not the mean.
    ## Hence below fitting uses standard deviation, not standard error of the mean.
    numavgs,firingbinsmeanList,firingbinserrList = read_morphfile(filename,fitted_mitral,NUM_REBINS)

    ########################## Initial values for the parameters
    if fit_type=='arb':
        params_filename = filename+'_params'+str(fitted_mitral)
    else:
        params_filename = filename+'_paramsFULLlin'+str(fitted_mitral)
    
    if firstrun or refit:
        params0 = []
        spikesmax = firingbinsmeanList.max()
        RA = firingbinsmeanList[-2] # odor A is last but one
        RB = firingbinsmeanList[0] # odor B is first
        Rair = firingbinsmeanList[-1] # air response is last
        # The initial parameters are for odor A followed by odor B
        # extend(): Don't add the list as an element but add the elements of the list
        
        if fit_type == 'arb':
            # the small value 0.001 should be put, else divide by zero errors in chi-sq!
            params0.extend([ ( inversesigmoid(0.998*RA[i]/spikesmax+0.001) - \
                inversesigmoid(0.998*Rair[i]/spikesmax+0.001) )/log81 for i in range(NUM_REBINS) ])
            params0.extend([ ( inversesigmoid(0.998*RB[i]/spikesmax+0.001) - \
                inversesigmoid(0.998*Rair[i]/spikesmax+0.001) )/log81 for i in range(NUM_REBINS) ])
            # initial params for the air vector # air is last
            params0.extend([ inversesigmoid(0.998*Rair[i]/spikesmax+0.001)/log81 for i in range(NUM_REBINS) ])
            params0.extend([0.0]*2*NUMWTS) # weights of mixtures
            # argument for the exp in sigmoidmax as per Mukund and Adil.
            # -1 gives match for generated data, -3 went into local minimum.
            params0.append(-1)
        else:
            params0.extend(RA - Rair)
            params0.extend(RB - Rair)
            params0.extend(Rair)

        ##### pure odor concentrations are not parameters.
        ##### They are set to (CA=1,CB=0) and (CA=0,CB=1) and act as normalization.
        ## if arbitrary fit_type, weights are also free params,
        ## if linear fit_type, weights are not free params.
        if fit_type == 'arb':
            ## take only the mixture values, not the start and end-1 points which are pure odors,
            ## nor end point which is pure air
            for i,(inputA,inputB) in enumerate(inputList[1:-2]):
                # to constrain weights between 0 and 1, sigmoid is used,
                # so use inversesigmoid to set initial value for weights
                params0[3*NUM_REBINS+i] = inversesigmoid(inputA)
                params0[3*NUM_REBINS+NUMWTS+i] = inversesigmoid(inputB)
    else:
        f = open(params_filename,'r')
        params0,chisq = pickle.load(f)
        f.close()

    ###################################### Fitting
    if not refit:
        params = array(params0) ## only use params, do not fit again
    else:
        ## args is a tuple! if only one element write (elem, )
        params = optimize.leastsq( chisqfunc, params0,
            args=(firingbinsmeanList, firingbinserrList, fit_type), full_output=1, maxfev=50000)
        params = params[0] # leastsq returns a whole tuple of stuff - errmsg etc.

    ## Calculate sum of squares of the chisqarray
    chisqarraysq = [i**2 for i in chisqfunc(params, firingbinsmeanList, firingbinserrList, fit_type)]
    chisq = reduce(lambda x, y: x+y, chisqarraysq)

    if refit:
        paramsfile = open(params_filename,'w')
        pickle.dump((params,chisq), paramsfile)
        paramsfile.close()

    ############################## Calculate fitted responses and return them
    
    if fit_type == 'arb':
        #### for the weights also, we use exactly what is done by Mukund and Adil in matlab:
        #### constrain weights to be between 0 and 1
        #### sort the weights to ensure monotonicity
        inputsA = [ constrain0to1(x) for x in params[3*NUM_REBINS:(3*NUM_REBINS+NUMWTS)] ]
        inputsA.extend([0.0,1.0])
        inputsA.sort() # in place sort
        inputsB = [ constrain0to1(x) for x in params[(3*NUM_REBINS+NUMWTS):(3*NUM_REBINS+2*NUMWTS)] ]
        inputsB.extend([0.0,1.0])
        inputsB.sort(reverse=True) # weights of odor B need to be used in reverse
        #### Mukund and Adil constrained sigmoidmax > ydatamax (note exp(x)>0.)
        sigmoidmax = firingbinsmeanList.max() + math.exp(params[3*NUM_REBINS+2*NUMWTS])
    else:
        ## *-operator unpacks the list which become args of zip()
        ## zip collects the i-th elements together of all the args.
        inputsA,inputsB = zip(*(inputList[:-1])) # leave out the last (0,0) air input

    fitted_responses = []
    Rair = params[2*NUM_REBINS:3*NUM_REBINS]
    for inpnum,(inputA,inputB) in enumerate(inputList[:-1]):
        if fit_type == 'arb':
            fitted_responses.append(\
                [ sigmoidmax*outputsigmoid( \
                inputsA[inpnum]*params[i] + inputsB[inpnum]*params[NUM_REBINS+i] + Rair[i]\
                ) for i in range(NUM_REBINS) ] )
        else:
            fitted_responses.append( rectifier( \
                inputsA[inpnum]*params[:NUM_REBINS] + \
                inputsB[inpnum]*params[NUM_REBINS:2*NUM_REBINS] + Rair ) )
    if fit_type == 'arb':
        fitted_responses.append([ sigmoidmax*outputsigmoid( Rair[i] ) \
                for i in range(NUM_REBINS) ] )
    else:
        fitted_responses.append( rectifier(Rair) )

    return (params,chisq,inputsA,inputsB,fitted_responses,numavgs,firingbinsmeanList,firingbinserrList)

def plot_example_onemit(ax1,ax2,fitted_mitral,mit_fit_params):
    bindt = RESPIRATION/float(NUM_REBINS)
    respiration2time = arange(RESPIRATION,2*RESPIRATION,bindt) + bindt/2.0

    params,chisq,inputsA,inputsB,fitted_responses,numavgs,firingbinsmeanList,firingbinserrList =\
        mit_fit_params
    print "Mit",fitted_mitral,"normalized chisq =",chisq
    brightness = 0.2
    num_morphs = len(inputList)-1
    for i,(inputA,inputB) in enumerate(inputList):
        ## The inputA acts to morph odor response from red to blue color
        ## air response in black
        ## if not a pure odor/air, bring down its brightness.
        if i==0: color,alpha = 'b',1.0
        elif i==num_morphs-1: color,alpha = 'r',1.0
        elif i==num_morphs: color,alpha = 'k',1.0
        else: color,alpha = (i/float(num_morphs),0,1.0-i/float(num_morphs)),brightness
        if i in [0,num_morphs-1,num_morphs]:
            simresponse = firingbinsmeanList[i]
            ## For the plots, show std error of the mean
            simerr = firingbinserrList[i]/sqrt(numavgs)
            ax1.fill_between(respiration2time,simresponse+simerr,simresponse-simerr,
                color=color,alpha=alpha*0.4,linewidth=0)
            ax1.plot(respiration2time,simresponse,\
                color=color,alpha=alpha,marker='.',markersize=marker_size,\
                linewidth=linewidth,clip_on=False)

    ##################### Plot fitted responses.
    ## RA + Rair
    line, = ax1.plot(respiration2time,fitted_responses[-2],\
        color='m',marker='+',markersize=marker_size,\
        linestyle='dashed', linewidth=linewidth, label='fit A',clip_on=False)
    line.set_dashes((3,1))
    ## RB + Rair
    line, = ax1.plot(respiration2time,fitted_responses[0],\
        color='c',marker='+',markersize=marker_size,\
        linestyle='dashed', linewidth=linewidth, label='fit B',clip_on=False)
    line.set_dashes((3,1))
    ## Rair
    line, = ax1.plot(respiration2time,fitted_responses[-1],\
        color=(0.5,0.5,0.5),marker='+',markersize=marker_size,\
        linestyle='dashed', linewidth=linewidth, label='fit air',clip_on=False)
    line.set_dashes((3,1))
    #title('Mitral %d responses & linear-sigmoid fit'%fitted_mitral,fontsize=24 )
    #axes_labels(ax,'respiratory phase bin','firing rate (Hz)',adjustpos=True)
    #ylim(ymin=-6, ymax=4)
    #legend()

    ################### Linearity of weights plot
    print 'weightsA =',inputsA
    print 'weightsB =',inputsB

    actualweights = [ wts[0] for wts in inputList[:-1]]
    ax2.plot(actualweights,arange(0.0,1.01,0.2),color='r',\
        marker='.',markersize=marker_size,clip_on=False,\
        linestyle='solid',linewidth=linewidth,label='linear A')
    ax2.plot(actualweights,arange(1.0,-0.01,-0.2),color='b',\
        marker='.',markersize=marker_size,clip_on=False,\
        linestyle='solid',linewidth=linewidth,label='linear B')
    line, = ax2.plot(actualweights,inputsA,color='m',clip_on=False,\
        marker='+',linestyle='dashed',markersize=marker_size,\
        linewidth=linewidth,label='weight odorA')
    line.set_dashes((3,1))
    line, = ax2.plot(actualweights,inputsB,color='c',clip_on=False,\
        marker='+',linestyle='dashed',markersize=marker_size,\
        linewidth=linewidth,label='weight odorB')
    line.set_dashes((3,1))
    #title( 'chisquare normalized = '+str(chisq) )
    maxerror = sqrt(sum(array([0.8,0.6,0.4,0.2])**2)/4.0) # max rms error
    ## normalized score = 1 - norm-ed rms error
    scoreA = 1 - sqrt( sum( (inputsA[1:-1]-arange(0.2,0.81,0.2))**2 )/4.0 )/maxerror
    scoreB = 1 - sqrt( sum( (inputsB[1:-1]-arange(0.8,0.19,-0.2))**2 )/4.0 )/maxerror
    #title( 'Linearity mitral %d: \nscoreA=%.2f, scoreB=%.2f'%(fitted_mitral,scoreA,scoreB), fontsize=24 )
    #axes_labels(ax2,'weight','fitted weight',adjustpos=True)
    #legend(loc='center right')

    ## beautify plots
    for ax in [ax1,ax2]:
        xmin,xmax,ymin,ymax = \
            beautify_plot(ax,x0min=False,y0min=True,xticksposn='bottom',yticksposn='left')

def plot_example_chisq():
    fig = figure(figsize=(columnwidth,columnwidth/2.0),dpi=300,facecolor='w') # 'none' is transparent
    if 'CHISQ_HIST' in sys.argv:
        axgrid = (2,3)
        ax5 = plt.subplot2grid(axgrid,(0,2),frameon=False)
        ax6 = plt.subplot2grid(axgrid,(1,2),frameon=False)
        import average_odor_morphs as chisq_hist
        ## chi-sq histograms for both non-lin and lin weights
        chisq_hist.plot_chisq_hist_paperfigure(ax5,ax6,'../results/odor_morphs'+dirextn)
    else: axgrid = (2,2)
    for fitted_mitral in fitted_mitral_list:
        mit_fit_params = fit_morphs(filename, fitted_mitral, fit_type=fit_type)

        ################# Plot simulated responses
        ax1 = plt.subplot2grid(axgrid,(fitted_mitral,0),frameon=False)
        #text(0.1,1.0,['A','C'][fitted_mitral],fontsize=label_fontsize,transform=ax1.transAxes)
        ax2 = plt.subplot2grid(axgrid,(fitted_mitral,1),frameon=False)
        #text(0.1,1.0,['B','D'][fitted_mitral],fontsize=label_fontsize,transform=ax2.transAxes)
        
        plot_example_onemit(ax1,ax2,fitted_mitral,mit_fit_params)

    fig.tight_layout()
    fig_clip_off(fig)
    subplots_adjust(left=0.1,top=0.92,bottom=0.15,wspace=0.5)
    fig.text(0.015,0.7,'firing rate (Hz)',fontsize=label_fontsize, rotation='vertical', transform=fig.transFigure)
    fig.text(0.15,0.025,'time (s)',fontsize=label_fontsize, transform=fig.transFigure)
    fig.text(0.35,0.7,'fitted weight',fontsize=label_fontsize, rotation='vertical', transform=fig.transFigure)
    fig.text(0.43,0.025,'ORN weight',fontsize=label_fontsize, transform=fig.transFigure)
    if 'SAVEFIG' in sys.argv:
        fig.savefig('../figures/sim_morphs.svg', bbox_inches='tight',dpi=fig.dpi)
        fig.savefig('../figures/sim_morphs.png', bbox_inches='tight',dpi=fig.dpi)

def plot_example(refit=True):
    fig = figure(figsize=(columnwidth*2./3.,columnwidth/2.0),dpi=300,facecolor='w') # 'none' is transparent
    axgrid = (2,2)
    for fitted_mitral in fitted_mitral_list:
        mit_fit_params = fit_morphs(filename, fitted_mitral, fit_type=fit_type, refit=refit)

        ################# Plot simulated responses
        ax1 = plt.subplot2grid(axgrid,(fitted_mitral,0),frameon=False)
        #text(0.1,1.0,['A','C'][fitted_mitral],fontsize=label_fontsize,transform=ax1.transAxes)
        ax2 = plt.subplot2grid(axgrid,(fitted_mitral,1),frameon=False)
        #text(0.1,1.0,['B','D'][fitted_mitral],fontsize=label_fontsize,transform=ax2.transAxes)
        
        plot_example_onemit(ax1,ax2,fitted_mitral,mit_fit_params)
        axes_labels(ax1,['','time (s)'][fitted_mitral],\
            ['firing rate (Hz)',''][fitted_mitral],xpad=0,ypad=0)
        axes_labels(ax2,['','ORN weight'][fitted_mitral],\
            ['fitted weight',''][fitted_mitral],xpad=0,ypad=0)
        ax1.yaxis.set_label_coords(-0.25,-0.3)
        ax2.yaxis.set_label_coords(-0.16,-0.3)

    fig.tight_layout()
    fig.subplots_adjust(hspace=0.2,wspace=0.4)
    fig_clip_off(fig)
    if 'SAVEFIG' in sys.argv:
        fig.savefig('../figures/sim_morphs.svg', bbox_inches='tight',dpi=fig.dpi)
        fig.savefig('../figures/sim_morphs.png', bbox_inches='tight',dpi=fig.dpi)

if __name__ == "__main__":
    if len(sys.argv) > 1:
        filename = sys.argv[1]
        post_pulses = filename.split('odor_morphs')[1]
        dirextn = post_pulses.split('/')[0]
        print 'directory extension is',dirextn
    else:
        print "Specify data file containing pickled list."
        sys.exit(1)

    ### old paper figure that shows example and chisq distribs
    #plot_example_chisq()
    ### OBSOLETE: PAPER FIGURE supplementary fig 1 that shows only example
    ### use average_odor_morphs.py to get the R2N distribution plot
    ##plot_example()
    ## For PAPER FIGURE supplementary Fig 5
    ## use average_odor_morphs.py to get the sqrt(R2N) distribution plot
    ## and it also calls plot_example_onemit() to plot one mitral example
    
    ## Plot the fitting result of the commandline data file, do not refit
    plot_example(refit=False)

    show()

Loading data, please wait...