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
#### OBSOLETE -- Priyanka does not fit air kernel! Use fit_odor_pulses_constair.py instead
#### This program also has no rebinning with 50ms bin and no padding of init kernel length with zeros.
## USAGE1: python2.6 fit_odor_pulses.py <../results/odor_pulses/pulseresult.pickle> 191.0
## USAGE2: python2.6 fit_odor_pulses.py <../results/odor_pulses/pulseresult.pickle> 191.0 <../results/odor_morphs/morphresult.pickle>

########## OBSOLETE ######### OBSOLETE ########## OBSOLETE ######### OBSOLETE ######## OBSOLETE ########
####### Use fit_odor_pulses_constair.py instead. Priyanka does not fit air kernel, nor A+B as in _constair.

from scipy import optimize
from scipy import special
import scipy.interpolate
from pylab import *
import pickle
import sys
import math

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

from stimuliConstants import * # has SETTLETIME, PULSE_RUNTIME, pulsebins, pulsebindt
from networkConstants import * # has central_glom
from data_utils import * # has axes_labels()
from analysis_utils import * # has read_morphfile() and various constants
import generate_firerates_odors as gen_frate_odors # has respiration_function()

iterationnum = 1
sigmoid = False
FitScaledORNkernels = False#True
SMOOTH_KERNELS = False

numpulses = RANDOM_PULSE_NUMS*2
pulsetlist = arange(pulsebindt/2.0,PULSE_RUNTIME,pulsebindt)
## very IMPORTANT to have the same dt for pulse and kernel
## when doing discrete convolution!
kerneldt = pulsebindt
kernel_size = int(kernel_time/kerneldt)
kerneltlist = arange(0.0,kernel_time,kerneldt)[:kernel_size]
#NUMWTS = NUMMIX-3 # remove the two pure odors and one pure air weights
#NUMPARAMS = 3*NUMBINS+2*NUMWTS+1 # two pure + one air response,
# weights of A and B for mixtures, max of output sigmoid
firstrun = True

## If you want to plot in style of Priyanka
fill_below = True

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

global iternum
iternum = 0

## for matching morphs
match_morph = False
## for the morphs, I don't use the NUMBINS in simset_odor.py,
## rather I rebin() with an overlap
NUM_REBINS = 20
resp_pulses = gen_frate_odors.respirationPulses
st_idx = gen_frate_odors.startidx
numt = gen_frate_odors.numt

def responseconvolve(pulses, kernel):
    ## pulsebindt is typically unchanged, pulse_time can change
    ## changes number of pulsebins in stimuliConstants.py
    decimatefactor = int(round(pulsebindt/FIRINGFILLDT)) # ensure that these are multiples
    pulses = pulses[::decimatefactor]
    ## don't use mode='same', as it takes only the valid part of the convolution.
    ## always multiply convolution_result by dt
    ## to leave response invariant on decimating stim and kernel simultaneously.
    responsecalc = convolve( pulses, kernel, mode='full' ) * pulsebindt
    ## bin0 of the simulated response is at time pulsebindt
    ## so the fitted response should also be shifted by a pulsebindt
    return array( responsecalc[1:pulsebins+1] )
    
def sigmoid_fn(valarray, max, half_point):
    ## ensure that valarray is a numpy array
    ## else where() returns empty tuple, and below exp() statement won't work.
    valarray = array(valarray)
    if sigmoid:
        ## exp() - logistic function - not used
        #return max/(1.0+exp(-valarray+half_point)) # numpy's exp works for arrays.
        ## thresholded erf() - used
        ## erf(0.5)~=0.5; set half_point as center of sigmoid.
        valarray[where(valarray<0)[0]] = 0
        posindices = where(valarray<0)[0]
        valarray[posindices] = max * special.erf( 0.5*valarray[posindices]/half_point )
        return valarray
    else:
        ## where value in valarray < 0, set it to zero.
        valarray[where(valarray<0)[0]]=0
        return valarray

#def inv_sigmoid_fn(valarray, max, half_point):
#    ## ensure that valarray is a numpy array
#    ## else where() returns empty tuple, and below exp() statement won't work.
#    valarray = array(valarray)
#    if sigmoid:
#        ## numpy has overloaded operators, also numpy ln is overloaded for arrays.
#        return half_point - log(max/valarray-1)
#    else:
#        return valarray

def SGfilter3(kernelR,kernelA,kernelB):
    ## savitsky-golay filter maintains high frequency signal, unlike moving average filter
    ## window_size must be odd, savitsky_golay copied from scipy cookbook online
    if int(kernel_size/2)%2==0: window_size = int(kernel_size/2)+1
    else: window_size = int(kernel_size/2)
    kernelR = savitzky_golay(kernelR, window_size, order=5)
    kernelA = savitzky_golay(kernelA, window_size, order=5)
    kernelB = savitzky_golay(kernelB, window_size, order=5)
    return kernelR, kernelA, kernelB

def chisqfunc(params, firingbinsmeanList, firingbinserrList,\
        pulseList, ret_array=True, ORNkernels=()):
    if not FitScaledORNkernels:
        kernelR = params[0:kernel_size]
        kernelA = params[kernel_size:2*kernel_size]
        kernelB = params[2*kernel_size:3*kernel_size]
        bgnd = params[3*kernel_size:3*kernel_size+1]
    else:
        scaleR = params[0]
        if scaleR<0: scaleR = 0.0
        scaleA = params[1]
        if scaleA<0: scaleA = 0.0
        scaleB = params[2]
        if scaleB<0: scaleB = 0.0
        kernelR = scaleR*ORNkernels[0]
        kernelA = scaleA*ORNkernels[1]
        kernelB = scaleB*ORNkernels[2]
        bgnd = params[3]
    if sigmoid:
        sigmoid_max = params[-2]
        sigmoid_threshold = params[-1]
    else:
        sigmoid_max = None
        sigmoid_threshold = None
    if SMOOTH_KERNELS:
        kernelR, kernelA, kernelB = SGfilter3(kernelR, kernelA, kernelB)
    ## These below should be numpy arrays,
    ## else + acts like append in python (not element-wise addition)!!!
    pulseair = responseconvolve(pulseList[0],kernelR)
    responsecalcair = sigmoid_fn(
        responseconvolve(pulseList[1],kernelR) + bgnd,\
        sigmoid_max, sigmoid_threshold)
    responsecalcA1 = sigmoid_fn(
        responseconvolve(pulseList[2],kernelA) + pulseair + bgnd,\
        sigmoid_max, sigmoid_threshold)
    responsecalcB1 = sigmoid_fn(
        responseconvolve(pulseList[3],kernelB) + pulseair + bgnd,\
        sigmoid_max, sigmoid_threshold)
    responsecalcA2 = sigmoid_fn(
        responseconvolve(pulseList[4],kernelA) + pulseair + bgnd,\
        sigmoid_max, sigmoid_threshold)
    responsecalcB2 = sigmoid_fn(
        responseconvolve(pulseList[5],kernelB) + pulseair + bgnd,\
        sigmoid_max, sigmoid_threshold)
    responsecalcAB = sigmoid_fn(
        responseconvolve(pulseList[6],kernelA) + responseconvolve(pulseList[7],kernelB) + \
        + pulseair + bgnd, sigmoid_max, sigmoid_threshold)
    ## fitted response is also shifted by a pulsebindt, to compare index to index
    chisqarray = [ (val-firingbinsmeanList[0][i])/firingbinserrList[0][i]\
        for i,val in enumerate(responsecalcair) ]
    chisqarray.extend( [ (val-firingbinsmeanList[1][i])/firingbinserrList[1][i]\
        for i,val in enumerate(responsecalcA1) ] )
    chisqarray.extend( [ (val-firingbinsmeanList[2][i])/firingbinserrList[2][i]\
        for i,val in enumerate(responsecalcB1) ] )
    chisqarray.extend( [ (val-firingbinsmeanList[3][i])/firingbinserrList[3][i]\
        for i,val in enumerate(responsecalcA2) ] )
    chisqarray.extend( [ (val-firingbinsmeanList[4][i])/firingbinserrList[4][i]\
        for i,val in enumerate(responsecalcB2) ] )
    chisqarray.extend( [ (val-firingbinsmeanList[5][i])/firingbinserrList[5][i]\
        for i,val in enumerate(responsecalcAB) ] )
        
    global iternum
    iternum += 1
    ## normalize chi-sq to the number of dof
    num_dof = float(firingbinsmeanList.size - params.size)
    if ret_array:
        if iternum % 10000 == 0:
            chisqarraysq = [i**2 for i in chisqarray]
            chisq = reduce(lambda x, y: x+y, chisqarraysq) / num_dof
            print chisq
        ## len(chisqarray) must be >= len(params):
        ## since leastsq() assumes it is solving M eqns in N unknowns where M>=N.
        return array(chisqarray) / sqrt(num_dof) # this must be squared and added for chi-sq
    else:
        chisq = reduce(lambda x, y: x+y**2, chisqarray) / num_dof
        if iternum % 10000 == 0:
            print chisq
        return chisq

def match_morphs_from_kernels(kernels,bgnd):
    kernelR = kernels[0]
    kernelA = kernels[1]
    kernelB = kernels[2]
    ## resp_pulses is peak-normalized to 1; the on-off random pulses also go to 1.
    ## Thus, no factor in response for scaling of stimulus
    ## However the pulse, response, kernels were using dt=kerneldt=pulsebindt
    ## But the kernel is a binned kernel, now interpolate to have dt = FIRINGFILLDT
    ## FIRINGFILLDT is the dt for respirationPulses and odor ORN frates
    interpol_type = 'cubic' # can have 'cubic', 'linear', etc.
    kernelR_fn = scipy.interpolate.interp1d(kerneltlist, kernelR, kind=interpol_type)
    kernelA_fn = scipy.interpolate.interp1d(kerneltlist, kernelA, kind=interpol_type)
    kernelB_fn = scipy.interpolate.interp1d(kerneltlist, kernelB, kind=interpol_type)
    ## make a finer dt for kernel = FIRINGFILLDT,
    ## need to subtract kerneldt from endpoint to keep within interpolation function range
    kerneltlist_finer = arange(0.0,kernel_time-kerneldt,FIRINGFILLDT)
    kernelR_full = [kernelR_fn(t) for t in kerneltlist_finer]
    kernelA_full = [kernelA_fn(t) for t in kerneltlist_finer]
    kernelB_full = [kernelB_fn(t) for t in kerneltlist_finer]
    ## always multiply convolution_result by dt to leave response invariant
    ## on decimating/un-decimating stim and kernel simultaneously.
    respirationPulses = gen_frate_odors.respirationPulses
    responseR = convolve(respirationPulses,kernelR_full,mode='full') * FIRINGFILLDT
    responseR = responseR[st_idx:st_idx+numt]
    #responseR = responseR[:numt]
    responseA = convolve(respirationPulses,kernelA_full,mode='full') * FIRINGFILLDT
    responseA = responseA[st_idx:st_idx+numt]
    #responseA = responseA[:numt]
    responseB = convolve(respirationPulses,kernelB_full,mode='full') * FIRINGFILLDT
    responseB = responseB[st_idx:st_idx+numt]
    #responseB = responseB[:numt]

    matched_responses = []
    for inpnum,(inputA,inputB) in enumerate(inputList[:-1]):
        ## ideally I should use the fitted input weights for this morph
        matched_response = inputA*responseA + inputB*responseB + responseR + bgnd
        ## clip values going negative
        matched_response[ where(matched_response<0)[0] ] = 0
        matched_responses.append( matched_response )
    ## clip values going negative
    responseR[ where(responseR<0)[0] ] = 0
    matched_responses.append( responseR +bgnd )
    return matched_responses

if __name__ == "__main__":
    if len(sys.argv) > 2:
        ######### load in the mitral pulse responses
        f = open(sys.argv[1],'r')
        #### mitral_responses_list[avgnum][pulsenum][mitralnum][binnum]
        #### mitral_responses_binned_list[avgnum][pulsenum][mitralnum][binnum]
        mitral_responses_list, mitral_responses_binned_list = pickle.load(f)
        f.close()
        ########## load in the stimuli
        seednum = sys.argv[2]
        f = open('../generators/firerates/firerates_2sgm_'+seednum+'.pickle','r')
        frateOdorList,fratePulseList,randomPulseList,\
            randomPulseStepsList,randomResponseList,kernels\
            = pickle.load(f)
        #del frateOdorList, fratePulseList
        pulseList = array(randomPulseList)
        pulseStepsList = array(randomPulseStepsList)
        ## randomResponseList[glomnum][frate,...]
        ORNfrateList = array(randomResponseList)[central_glom]
        f.close()
        ######### should we match mitral morph responses using kernel
        if len(sys.argv)>3:
            match_morph = True
            morph_filename = sys.argv[3]
    else:
        print "Specify data file containing pickled mitral responses, and ORN frate seed."
        sys.exit(1)
    ## very important to convert to numpy array(),
    ## else where() below does not find any element satisfying condition.
    mitral_responses_binned_list = array(mitral_responses_binned_list)
    numtrials = len(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).
    ## NO! The model is for an individual response, not for the mean response!
    #mitral_responses_se = mitral_responses_std/sqrt(numtrials)
    
    for fitted_mitral in fitted_mitrals:
        ## take the odor responses of the mitral to be fitted
        firingbinsmeanList = mitral_responses_mean[:,fitted_mitral]
        firingbinserrList = mitral_responses_std[:,fitted_mitral]

        ### force the errors to be 1.2*sqrt(mean) which is what Adil observes.
        ### The model seems to have constant noise independent of the mean firing rate.
        ### standard error of the mean, hence /sqrt(numtrials)
        #firingbinserrList = 1.2*firingbinsmeanList/sqrt(numtrials)

        errcut = 1e-20
        #errcut = 1.0/pulsebindt/sqrt(numtrials) # errcut = 3.33Hz!
        #errcut = 10.0 #Hz
        ## put in a minimum error, else divide by zero problems.
        ## find the minimum error >= errcut
        largeerrors = firingbinserrList[where(firingbinserrList>errcut)]
        if largeerrors!=[]: errmin = largeerrors.min()
        else: errmin = errcut
        ## numpy where(), replace by errmin,
        ## all those elements in firingbinsList which are less than errmin 
        firingbinserrList = where(firingbinserrList>errcut, firingbinserrList, errmin)

#---------------------------- fitting ------------------------------------------------

        decimatefactor = int(kerneldt/FIRINGFILLDT)
        ## kernels are only kernel_size long.
        ## numpy append() is like python's extend() - makes flat not nested lists
        ## I take only kernelsize length after decimating the ORN kernel 
        ## by taking only every int(kerneldt/FIRINGFILLDT) element
        ## since convolve_result is always multiplied by dt,
        ## the response is invariant of decimating the stimulus and the kernel simultaneously.
        ## NOTE: We are using ORN kernel as mitral kernel, typical mitral
        ## firing rates are 5 times larger. So put in factor of 5!!!
        MitralORNfiringratio = 5.0
        linORNkernelR = array(kernels[0][::decimatefactor][:kernel_size])*MitralORNfiringratio
        linORNkernelA = array(kernels[1][::decimatefactor][:kernel_size])*MitralORNfiringratio
        linORNkernelB = array(kernels[2][::decimatefactor][:kernel_size])*MitralORNfiringratio
        ## Initial values for the parameters
        if firstrun:
            ORNkernels = (linORNkernelR, linORNkernelA, linORNkernelB)
            if not FitScaledORNkernels:
                params = linORNkernelR # ORN kernelR is initial kernelR
                params = append(params,linORNkernelA) # ORN kernelA is initial kernelA
                params = append(params,linORNkernelB) # ORN kernelB is initial kernelB
                ## constant background mitral firing w/o air and odor due to ORN bgnd.
                params = append(params,ORN_BGND_MEAN*MitralORNfiringratio)
            else:
                ## initial scaling of the ORN kernels to fit mitral kernels
                params = array([1.0,1.0,1.0])
                params = append(params,0.0) # mitral background firing rate
            #figure()
            #plot(linORNkernelR,color=(0,0,0))
            #plot(linORNkernelA,color=(1,0,0))
            #plot(linORNkernelB,color=(0,1,0))
            #figure()
            #plot(responseconvolve(pulseList[5],linORNkernelB))
            ##show()
            ##sys.exit()
            if sigmoid:
                ## I don't have a steepness param,\
                ## else there will be redundancy in kernel scaling vs steepness.
                params = append(params, [120.0,60.0]) # 120Hz max and 60Hz threshold for sigmoid
        else:
            f = open(sys.argv[1]+'_params'+str(fitted_mitral),'r')
            sigmoid,chisq,params = pickle.load(f)
            f.close()
            ORNkernels = ()

        ## comment next fitting line if you just want to plot parameters.
        ## args is a tuple! if only one element write (elem, )

        #########
        ### fmin_powell is faster than fmin (simplex) but slower than leastsq
        ### However it does not enter into the very bad local minima of leastsq
        ### WARNING! But often it is worse at fitting than leastsq()
        #iternum = 0
        #params = optimize.fmin_powell( chisqfunc, params, \
        #    args=(firingbinsmeanList, firingbinserrList, pulseList, False),\
        #    full_output=1, maxfun=50000, maxiter=50000, ftol = 1 )
        #chisq = params[1]
        #print "The normalized chisq for mitral",fitted_mitral,"after fmin_powell =",chisq
        #params = params[0]

        #########
        ## leastsq() uses a modified version of Levenberg-Marquardt algorithm
        ## it optimizes M equations in N unknowns, where M>=N.
        ## Hence chisqfunc must return M numbers in an array
        ## which must be >= the number of params N in params0, 
        ## else: 'TypeError: Improper input parameters.'
        iternum = 0
        params = optimize.leastsq( chisqfunc, params,\
            args=(firingbinsmeanList, firingbinserrList, pulseList, True, ORNkernels),\
            full_output=1, maxfev=100000, ftol=1e-100, xtol=1e-100 )
        print params[3] # print the status message
        params = params[0] # take only fitted params; leastsq returns a tuple with errmsg, etc.
        ### Calculate sum of squares of the chisqarray
        chisqarraysq = [i**2 for i in chisqfunc(params,\
            firingbinsmeanList, firingbinserrList, pulseList, True, ORNkernels)]
        #chisqnumpy = array(chisqarraysq)
        #pos = where(chisqnumpy>100)[0]
        #lenf = len(firingbinsmeanList[0])
        #print "The [(pulsenum,index),...] where chisq>100 are",[ (i/lenf,i%lenf) for i in pos]
        chisq = reduce(lambda x, y: x+y, chisqarraysq)
        print "The normalized chisq for mitral",fitted_mitral,"after leastsq =",chisq
        
        ##########
        ### minimize is a wrapper for all the various fitting algorithms!
        #iternum = 0
        #fit_method = "Nelder-Mead"#"CG"#"Powell"#"Nelder-Mead"
        #params,info = optimize.minimize( chisqfunc, params,\
        #    args=(firingbinsmeanList, firingbinserrList, pulseList, False, ORNkernels),\
        #    method = fit_method, full_output=True,\
        #    options={'maxfev':100000, 'ftol':1e-100, 'xtol':1e-100, 'maxfun':500000, 'maxiter':500000} )
        #print info['message']
        #chisq = info['fun']
        #print "The normalized chisq for mitral",fitted_mitral,"after",fit_method,"=",chisq

        if not FitScaledORNkernels:
            kernelR = params[0:kernel_size]
            kernelA = params[kernel_size:2*kernel_size]
            kernelB = params[2*kernel_size:3*kernel_size]
            bgnd = params[3*kernel_size:3*kernel_size+1]
        else:
            kernelR = ORNkernels[0]*params[0]
            kernelA = ORNkernels[1]*params[1]
            kernelB = ORNkernels[2]*params[2]
            print params[0],params[1],params[2]
            bgnd = params[3]
        if SMOOTH_KERNELS:
            kernelR, kernelA, kernelB = SGfilter3(kernelR, kernelA, kernelB)

        if sigmoid:
            sigmoid_max = params[-2]
            sigmoid_threshold = params[-1]
        else:
            sigmoid_max = None
            sigmoid_threshold = None
        pulseair = responseconvolve(pulseList[0],kernelR)
        
        #chisq = chisqfunc(params,firingbinsmeanList, firingbinserrList, pulseList)

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


#---------------------------- matching morphs using kernels -----------------------------------


        if match_morph:
            fitkernels = (kernelR,kernelA,kernelB)
            matched_responses = match_morphs_from_kernels(fitkernels,bgnd)

#---------------------------- plot matched responses -------------------------------------------

            mfiringbinsmeanList,mfiringbinserrList = \
                read_morphfile(morph_filename,fitted_mitral,NUM_REBINS)
            bindt = RESPIRATION/float(NUM_REBINS)
            responsetlist = arange( SETTLETIME+RESPIRATION+bindt, \
                SETTLETIME+2*RESPIRATION+bindt/2.0, bindt)
            figure()
            firingtsteps = gen_frate_odors.firingtsteps
            plot(firingtsteps, matched_responses[0], color=(0,1,0))
            plot(firingtsteps, matched_responses[-2], color=(1,0,0))
            plot(firingtsteps, matched_responses[-1], color=(0,0,0))
            plot(responsetlist, mfiringbinsmeanList[0], color=(0,1,0.5))
            plot(responsetlist, mfiringbinsmeanList[-2], color=(1,0,0.5))
            plot(responsetlist, mfiringbinsmeanList[-1], color=(0,0,0.5))
            #plot(firingtsteps, frateOdorList[central_glom][0], color=(0,1,0.5))
            #plot(firingtsteps, frateOdorList[central_glom][-2], color=(1,0,0.5))
            #plot(firingtsteps, frateOdorList[central_glom][-1], color=(0,0,0.5))
            #plot(firingtsteps, gen_frate_odors.respirationPulses[st_idx:st_idx+numt], color=(1,1,0.5))

#---------------------------- plot kernels and pulse responses -----------------------------------------


        ############################## plot the kernels
        fig = figure(facecolor='w') # 'none' is transparent
        ## A super axes to set a common ylabel
        bigAxes = axes(frameon=False) # hide frame
        xticks([]) # don't want to see any ticks on this axis
        yticks([])
        text(-0.125,0.3,'arb units', fontsize=24, rotation='vertical')

        #ax = fig.add_subplot(111)
        ax = fig.add_subplot(3,1,1)
        ax.set_yticks([0,75])
        ax.set_yticklabels(['0','75'])
        for label in ax.get_yticklabels():
            label.set_fontsize(24)
        axes_off(ax,x=True,y=False)
        ax.plot(kerneltlist, kernelA, color=(1,0,0), marker=',',\
            linestyle='solid',linewidth=2,label='OdorA kernel')
        ax.plot(kerneltlist, linORNkernelA, color=(1,0,1), marker=',',\
            linestyle='solid',linewidth=2,label='5 * ORN kernel')
        #ax.set_ylim(-10,25)
        #ax.set_xlim(0,2)
        biglegend()

        ax = fig.add_subplot(3,1,2)
        ax.set_yticks([-50,0,50])
        ax.set_yticklabels(['-50','0','50'])    
        for label in ax.get_yticklabels():
            label.set_fontsize(24)
        axes_off(ax,x=True,y=False)
        ax.plot(kerneltlist, kernelB, color=(0,1,0), marker=',',\
            linestyle='solid',linewidth=2,label='OdorB kernel')
        ax.plot(kerneltlist, linORNkernelB, color=(0,1,1), marker=',',\
            linestyle='solid',linewidth=2,label='5 * ORN kernel')
        #ax.set_xlim(0,2)
        #ax.set_ylim(-20,20)
        biglegend()

        ax = fig.add_subplot(3,1,3)
        ax.set_yticks([0,50])
        ax.set_yticklabels(['0','50'])    
        for label in ax.get_yticklabels():
            label.set_fontsize(24)
        ax.set_xticks([0,2])
        ax.set_xticklabels(['0','2'])    
        ax.plot(kerneltlist, kernelR, color=(0,0,0), marker=',',\
            linestyle='solid',linewidth=2,label='Air kernel')
        ax.plot(kerneltlist, linORNkernelR, color=(0,0,1), marker=',',\
            linestyle='solid',linewidth=2,label='5 * ORN kernel')
        axes_labels(ax,'time (s)','',adjustpos=False)
        #ax.set_xlim(0,2)
        #ax.set_ylim(-10,20)
        biglegend()
        fig.suptitle('Mitral '+str(fitted_mitral)+\
            ': odor and air kernels',fontsize=24)

        ############################### plot the responses and the fits
        fig = figure(facecolor='w')
        ## A super axes to set a common ylabel
        bigAxes = axes(frameon=False) # hide frame
        xticks([]) # don't want to see any ticks on this axis
        yticks([])
        #ylabel('firing rate (Hz)', fontsize=24)
        text(-0.125,0.25,'firing rate (Hz)', fontsize=28, rotation='vertical')
        for pulsenum in range(numpulses):
            sister_ratio = 0#(mitnum%MIT_SISTERS)/float(MIT_SISTERS)
            ax = fig.add_subplot(6,1,pulsenum+1)
            ################# random pulses and ORN responses
            xpulse = 50
            randompulsetime = arange(0,PULSE_RUNTIME+1e-10,FIRINGFILLDT)
            if pulsenum in [0]: # air response does not have air pedestal
                if fill_below:
                    fill_between(randompulsetime,xpulse*pulseStepsList[pulsenum+1]+
                        ORN_BGND_MEAN,linewidth=0,color=(0,0.9,0.9),alpha=0.4)
                else:
                    plot(randompulsetime,2*pulseList[pulsenum+1]+ORN_BGND_MEAN,\
                        linewidth=2, label='Air waveform')
                    plot(randompulsetime,ORNfrateList[pulsenum+1]+ORN_BGND_MEAN,\
                        linewidth=2, label='Receptors response')
            elif pulsenum in [1,2,3,4]: # odor A or B
                if fill_below:
                    if pulsenum in [1,3]: col = (1,0,0)
                    if pulsenum in [2,4]: col = (0,1,0)
                    fill_between(randompulsetime,xpulse*pulseStepsList[pulsenum+1]+
                        pulseList[0]+ORN_BGND_MEAN,linewidth=0,color=col,alpha=0.4)
                else:
                    plot(randompulsetime,2*pulseList[pulsenum+1]+pulseList[0]+ORN_BGND_MEAN,
                        linewidth=2,label='Air+Odor waveform')
                    plot(randompulsetime,ORNfrateList[pulsenum+1]+ORNfrateList[0]+ORN_BGND_MEAN,
                        linewidth=2,label='Receptors response')
            elif pulsenum in [5]: # odor A & odor B
                if fill_below:
                    fill_between(randompulsetime,xpulse*pulseStepsList[pulsenum+1],\
                        color=(1,0,0),linewidth=0,alpha=0.4)
                    fill_between(randompulsetime,xpulse*pulseStepsList[pulsenum+2],\
                        color=(0,1,0),linewidth=0,alpha=0.4)
                else:
                    plot(randompulsetime,pulseList[0]+ORN_BGND_MEAN,\
                        linewidth=2, label='Air waveform')
                    plot(randompulsetime,2*pulseList[pulsenum+1],\
                        linewidth=2, label='Odor A waveform')
                    plot(randompulsetime,2*pulseList[pulsenum+2],\
                        linewidth=2, label='Odor B waveform')
                    plot(randompulsetime,ORNfrateList[pulsenum+2]+ORNfrateList[pulsenum+1]\
                        +ORNfrateList[0]+ORN_BGND_MEAN, linewidth=2, label='Receptors response')
            ################### Plot the simulated responses
            if fill_below:
                ## numpy array, hence adds element by element
                fill_between(pulsetlist,
                    firingbinsmeanList[pulsenum]+firingbinserrList[pulsenum],
                    firingbinsmeanList[pulsenum]-firingbinserrList[pulsenum],
                    color=(0.7,0.7,0.7))
                plot(pulsetlist,firingbinsmeanList[pulsenum],linewidth=2,color=(0.3,0.3,0.3))
            else:
                errorbar(pulsetlist,y=firingbinsmeanList[pulsenum],\
                    yerr=firingbinserrList[pulsenum],\
                    color=(pulsenum/float(numpulses),1-pulsenum/float(numpulses),sister_ratio),\
                    marker='+',linestyle='solid',linewidth=2,label='Mitral response')
            ################## Plot the fitted responses
            if pulsenum == 0:
                response = sigmoid_fn(responseconvolve(pulseList[pulsenum+1],kernelR) + bgnd,\
                    sigmoid_max, sigmoid_threshold)
                titlestr = 'Air random pulse'
            elif pulsenum in [1,3]:
                response = sigmoid_fn(responseconvolve(pulseList[pulsenum+1],kernelA) + pulseair + bgnd,\
                    sigmoid_max, sigmoid_threshold)
                titlestr = 'Odor A random pulse'
            elif pulsenum in [2,4]:
                response = sigmoid_fn(responseconvolve(pulseList[pulsenum+1],kernelB) + pulseair + bgnd,\
                    sigmoid_max, sigmoid_threshold)
                titlestr = 'Odor B random pulse'
            else:
                response = responseconvolve(pulseList[pulsenum+1],kernelA) + pulseair
                response += responseconvolve(pulseList[pulsenum+2],kernelB)
                response += bgnd
                response = sigmoid_fn(response, sigmoid_max, sigmoid_threshold)
                titlestr = 'Odor A and odor B random pulse'
            if fill_below:
                plot(pulsetlist,response,linestyle='solid',linewidth=2.0,color=(0,0,0),
                    label=['air','A','B','A','B','A&B'][pulsenum])
                lgd = legend()
                for k in lgd.get_lines():
                    k.set_linewidth(0)
                lgd.draw_frame(False)
                ltext  = lgd.get_texts()
                for l in ltext:
                    l.set_fontsize(20)
                #fig.setp(ltext, fontsize=20)
            else:
                plot(pulsetlist,response,marker='o',linestyle='dashed',
                    linewidth=2.0, label='Linear fit')
                #title('Odor morphs w/ smooth fit',fontsize=24 )
                #title( 'pulsenum = '+str(pulsenum)+', chisquare normalized = '+str(chisq) )
                title(titlestr, fontsize=24)
                lgd = legend()
                ltext  = lgd.get_texts()
                for l in ltext:
                    l.set_fontsize(8)
                ax.set_ylim(0,xpulse)
            ax.set_yticks([xpulse])
            ax.set_yticklabels([str(xpulse)])
            for label in ax.get_yticklabels():
                label.set_fontsize(24)
            ax.set_xlim(0,9)
            ax.set_xticks([])
        ax.set_xticks([0,9])
        ax.set_xticklabels(['0','9'])    
        axes_labels(ax,'time (s)','',adjustpos=False)

    show()

Loading data, please wait...