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_concs.py ../results/odor_morphs/2011-01-13_odormorph_SINGLES_JOINTS_PGS.pickle
## This is only for concentration series, only one of odorAon or odorBon should be True

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 simset_odor import * # has NUMBINS
from networkConstants import * # has central_glom
from sim_utils import * # has rebin() to alter binsize

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

errcut = 1e-20
iterationnum = 1

NUMBINS = 17 # I override the NUMBINS in simset_odor above, and I rebin() below
bin_width_time = 2.0*RESPIRATION/NUMBINS # smooths / overlapping bins. non-overlapping would be RESPIRATION/NUMBINS

## index of the odor that is on, in inputList
if odorAon:
    odoridx = 0
    sortreverse = False
else:
    odoridx = 1
    sortreverse = True

NUMMIX = len(inputList)
NUMWTS = NUMMIX-3 # remove the two pure odors and one pure air weights
NUMPARAMS = 2*NUMBINS+1*NUMWTS+1 # one pure + one air response, weights of A or B for conc series, max of output sigmoid
firstrun = True

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

log81 = math.log(81)

def constrain0to1(x):
    return math.exp(x)/(1+math.exp(x))

# 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 chisqfunc(params, ydata, errdata):
    R = params[0:NUMBINS]
    Rair = params[NUMBINS:2*NUMBINS]
    #### 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
    inputs = [ constrain0to1(x) for x in params[2*NUMBINS:(2*NUMBINS+NUMWTS)] ]
    ## important to put these in else along with sort(),
    ## weights saturate at 0.9 or so rather than at 1.0
    inputs.extend([0.0,1.0]) # for pure odors
    inputs.sort(reverse=sortreverse) # in place sort
    inputs.append(0.0) # for air - keep this after sort!
    #### Mukund and Adil constrained sigmoidmax > ydatamax (note exp(x)>0.)
    sigmoidmax = ydata.max() + exp(params[2*NUMBINS+NUMWTS])

    global iterationnum
    chisqarray = [0.0]
    for i,input in enumerate(inputList):
        C = inputs[i]
        for bin in range(NUMBINS):
            Rmix = sigmoidmax*outputsigmoid( Rair[bin] + C*R[bin] )
            chisqarray.append( (ydata[i][bin] - Rmix)/errdata[i][bin] ) # divide by error to do chi-square fit
            
    chisqarray = array(chisqarray) / (ydata.size - NUMPARAMS) # normalized chi-sq to the number of dof
    if iterationnum%100==0: print 'iteration number =',iterationnum
    iterationnum += 1
    return chisqarray
    
def fit_morphs(filename, fitted_mitral):
    f = open(filename,'r')
    #### mitral_responses_list[avgnum][odornum][mitralnum][binnum]
    #### mitral_responses_binned_list[avgnum][odornum][mitralnum][binnum]
    mitral_responses_list, mitral_responses_binned_list = pickle.load(f)
    f.close()

    ###################### Input conditioning
    mitral_responses_binned_list = \
        rebin(mitral_responses_list, NUMBINS, bin_width_time)
    #### very important to convert to numpy array,
    #### else where() below returns empty list.
    mitral_responses_binned_list = array(mitral_responses_binned_list)
    mitral_responses_mean = mean(mitral_responses_binned_list, axis=0)
    mitral_responses_std = std(mitral_responses_binned_list, axis=0)
    # since I fit the mean response, I must use standard error/deviation of the _mean_
    # = standard deviation of a repeat / sqrt(num of repeats).
    NUMAVGs = len(mitral_responses_binned_list)
    mitral_responses_se= mitral_responses_std/sqrt(NUMAVGs)
    # take the odor responses of the mitral to be fitted
    firingbinsmeanList = mitral_responses_mean[:,fitted_mitral]
    firingbinserrList = mitral_responses_se[:,fitted_mitral]
    # put in a minimum error, else divide by zero problems.
    errmin = firingbinserrList[where(firingbinserrList>errcut)].min() # find the minimum error > errcut
    # numpy where(), replace by errmin, all those elements in firingbinsList which are less than errmin 
    firingbinserrList = where(firingbinserrList>errcut, firingbinserrList, errmin)

    ########################## Initial values for the parameters
    if firstrun:
        params0 = []
        spikesmax = firingbinsmeanList.max()
        if odorAon: R = firingbinsmeanList[-2] # odor A is last but one
        else: R = firingbinsmeanList[0] # odor B is first
        Rair = firingbinsmeanList[-1] # air response is last
        # The initial parameters are for odor 
        # the small value 0.001 should be put, else divide by zero errors in chi-sq!
        # extend(): Don't add the list as an element but add the elements of the list
        params0.extend([ ( inversesigmoid(0.998*R[i]/spikesmax+0.001) - \
            inversesigmoid(0.998*Rair[i]/spikesmax+0.001) )/log81 for i in range(NUMBINS) ])
        params0.extend([ inversesigmoid(0.998*Rair[i]/spikesmax+0.001)/log81 for i in range(NUMBINS) ]) # air is last
        # initial params for the air vector
        params0.extend([0.0]*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)
        ##### pure odor concentrations are not parameters.
        ##### They are set to (CA=1,CB=0) and (CA=0,CB=1) and act as normalization.
        # 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,input in enumerate(inputList[1:-2]):
            # to constrain weights between 0 and 1, sigmoid is used,
            # so use inversesigmoid to set initial value
            input = input[odoridx]
            params0[2*NUMBINS+i] = inversesigmoid(input)
    else:
        f = open(sys.argv[1]+'_params','r')
        params0 = pickle.load(f)
        f.close()

    ###################################### Fitting
    #params = array(params0) ## uncomment this line and comment next fitting line if you just want to plot parameters.
    # args is a tuple! if only one element write (elem, )
    params = optimize.leastsq( chisqfunc, params0,
        args=(firingbinsmeanList, firingbinserrList), full_output=1, maxfev=50000 )
    print params
    params = params[0] # leastsq returns a whole tuple of stuff - errmsg etc.

    paramsfile = open(sys.argv[1]+'_params','w')
    pickle.dump(params, paramsfile)
    paramsfile.close()

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

    #### Mukund and Adil constrained sigmoidmax > ydatamax (note exp(x)>0.)
    sigmoidmax = firingbinsmeanList.max() + math.exp(params[2*NUMBINS+NUMWTS])
    
    #### 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
    inputs = [ constrain0to1(x) for x in params[2*NUMBINS:(2*NUMBINS+NUMWTS)] ]
    inputs.extend([0.0,1.0])
    inputs.sort(reverse=sortreverse) # in place sort

    fitted_responses = []
    for inpnum,input in enumerate(inputList[:-1]):
        fitted_responses.append(\
            [ sigmoidmax*outputsigmoid( \
            inputs[inpnum]*params[i] + params[NUMBINS+i]\
            ) for i in range(NUMBINS) ] )
    fitted_responses.append( [ sigmoidmax*outputsigmoid( params[NUMBINS+i] )\
            for i in range(NUMBINS) ] )

    return (params,chisq,inputs,fitted_responses,firingbinsmeanList,firingbinserrList)

if __name__ == "__main__":
    if len(sys.argv) > 1:
        filename = sys.argv[1]
    else:
        print "Specify data file containing pickled list."
        sys.exit(1)

    for fitted_mitral in fitted_mitral_list:
        params,chisq,inputs,fitted_responses,firingbinsmeanList,firingbinserrList\
            = fit_morphs(filename, fitted_mitral)
        print "Mit",fitted_mitral,"chisq =",chisq

        ################# Plot simulated responses
        fig = figure(facecolor='w') # 'none' is transparent
        ax = fig.add_subplot(111)
        pure_i = 0
        brightness = 0.1
        for i,(inputA,inputB) in enumerate(inputList):
            ## The inputA/inputB acts to morph odor response from red/green to black (air)
            ## if not a pure odor/air, bring down its brightness.
            if odorAon: pureslist = [len(inputList)-2,len(inputList)-1]
            else: pureslist = [0,len(inputList)-1]
            if i not in pureslist:
                ax.errorbar(x=range(NUMBINS),y=firingbinsmeanList[i],yerr=firingbinserrList[i],\
                    color=( (1-brightness)+brightness*inputA, (1-brightness)+brightness*inputB, (1-brightness) ),\
                    marker='+',linestyle='solid')
            else:
                ax.errorbar(x=range(NUMBINS),y=firingbinsmeanList[i],yerr=firingbinserrList[i],\
                    color=(inputA,inputB,0),marker='+',linestyle='solid',\
                    label=['odor A','odor B','air'][pure_i])
                pure_i += 1

        ##################### Plot fitted responses.
        if odorAon:
            # RA + Rair
            ax.plot(fitted_responses[-2],\
                color=(1,0,1),marker='o',linestyle='dashed', linewidth=2.0, label='fit A')
        else:
            # RB + Rair
            ax.plot(fitted_responses[0],\
                color=(0,1,1),marker='o',linestyle='dashed', linewidth=2.0, label='fit B')
        # Rair
        ax.plot(fitted_responses[-1],\
            color=(0.5,0.5,0.5),marker='o',linestyle='dashed', linewidth=2.0, label='fit air')
        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 'weights =',inputs

        fig2 = figure(facecolor='w') # none is transparent
        ax2 = fig2.add_subplot(111)
        maxerror = sqrt(sum(array([0.8,0.6,0.4,0.2])**2)/4.0) # max rms error
        actualweights = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]
        if odorAon:
            ax2.plot(actualweights,inputs,color=(1,0,0),marker='+',linestyle='solid',label='weight odor A')
            ax2.plot(actualweights,arange(0.0,1.01,0.2),color=(1,0,1),marker='+',linestyle='dashed',label='linear A')
            ## normalized score = 1 - norm-ed rms error
            score = 1 - sqrt( sum( (inputs[1:-1]-arange(0.2,0.81,0.2))**2 )/4.0 )/maxerror
        else:
            ax2.plot(actualweights,inputs,color=(0,1,0),marker='+',linestyle='solid',label='weight odor B')
            ax2.plot(actualweights,arange(1.0,-0.01,-0.2),color=(0,1,1),marker='+',linestyle='dashed',label='linear B')
            ## normalized score = 1 - norm-ed rms error
            score = 1 - sqrt( sum( (inputs[1:-1]-arange(0.8,0.19,-0.2))**2 )/4.0 )/maxerror
        ax2.plot(actualweights,arange(0.0,1.01,0.2),color=(1,0,1),marker='+',linestyle='dashed',label='linear')
        #title( 'chisquare normalized = '+str(chisq) )
        title( 'Linearity mitral %d: \nscore=%.2f'%(fitted_mitral,score), fontsize=24 )
        axes_labels(ax2,'weight','fitted weight',adjustpos=True)
        legend(loc='center right')
    
    show()

Loading data, please wait...