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 -*-
# ALL SI UNITS
# milliMolar is same as mol/m^3

## USAGE: python2.6 average_odor_scaledpulses_noregress.py

import os,sys,math,string
import os.path
import pickle
import subprocess
cwd = os.getcwd() # current working directory

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

from stimuliConstants import * # has RESPIRATION
from pylab import * # part of matplotlib that depends on numpy but not scipy
from scipy import stats
from scipy import signal # for gaussian smoothing
from scipy import interpolate # for interp1d
from scipy import optimize # for fsolve
from sim_utils import * # has rebin() and imports data_utils.py for axes_off()
from analysis_utils import * # has rebin_mean() and read_pulsefile()
from calc_corrs import * # has calc_corrs()
## BE CAREFUL: use _constair_separateodors below, others are obsolete.
import fit_scaledpulses_noregress as fit_sp # has fit_plot_scaledpulses and ref pulse constants
import average_odor_morphs as corr_utils # has get_filename() for morphs

IN_VIVO = True
directed = True
FRAC_DIRECTED = 0.01 ## overrides the one in networkConstants.py (came in via sim_utils.py)
## Below two overide the variables that came in from stimuliConstantsMinimal.py via stimuliConstants.py
NONLINEAR_ORNS = False
NONLINEAR_TYPE = 'P' # P for primary glom non-linear, L for lateral gloms non-linear
_scaledWidth = 0.2 # s # overrides that in stimuliConstantsMinimal.py
num_scalings = len(pulseList)-1

#fullfilename = '../results/odor_morphs/morphs_random'
#if NONLINEAR_ORNS: fullfilename += 'NL'+NONLINEAR_TYPE
#fullfilename += '.pickle'
#fullfile = open(fullfilename,'r')
#morphs = pickle.load(fullfile)
#fullfile.close()

PLOTRESP_NUM = 1 # whether to plot 2 respiration cycles or 1
NUMBINS = 5
BIN_WIDTH_TIME = RESPIRATION/NUMBINS
bindt = RESPIRATION/float(NUMBINS)
## I take the last PLOTRESP_NUM of respiration cycles out of NUM_RESPS simulated
responsetlist = arange( SETTLETIME+(NUM_RESPS-PLOTRESP_NUM)*RESPIRATION+bindt/2.0, \
    SETTLETIME+NUM_RESPS*RESPIRATION, RESPIRATION/NUMBINS )

min_frate_cutoff = 0.25 # Hz

salient = False#True
if salient:
    stim_seeds = [-1,-19]#range(-1,-37,-1)#[-1,-2,-3,-4,-5,-6,-7,-8,-9]
    num_gloms_list = [3]
    ## inh_options = [ (no_singles,no_joints,no_lat,no_PGs,varyRMP), ... ]
    ## in order,below options are:
    ## all cells; no lat; no joints, varyRMP; no PGs; no singles + no joints, only mitrals
    #inh_options = [ (0,(False,False,False,False,False)), (1,(False,False,True,False,False)), \
    #    (2,(False,True,False,False,False)), (3,(False,False,False,False,True)), \
    #    (4,(False,False,False,True,False)), (5,(True,True,False,False,False)), \
    #    (6,(True,True,False,True,False))]
    inh_options = [ (0,(False,False,False,False,False)), (1,(False,False,True,False,False)) ]
else:
    #stim_seeds = [157.0,160.0,190.0,191.0,212.0,441.0]
    #num_gloms_list = [5,2]
    stim_seeds = arange(750.0,800.0,1.0)#[157.0,160.0,190.0,191.0]
    num_gloms_list = [3]
    ## inh_options = [ (no_singles,no_joints,no_lat,no_PGs,varyRMP), ... ]
    ## in order,below options are: all cells; no lat; no joints, varyRMP; no PGs; only mitrals
    inh_options = [ (0,(False,False,False,False,False)), (1,(False,False,True,False,False)), \
        (2,(True,False,False,False,False)), (3,(True,True,False,False,True)), \
        (6,(False,False,False,True,False)), (8,(True,True,False,True,False))]
net_seeds = [200.0]

def get_filename(netseed,stimseed,inh,ngi,\
        nl_orns=NONLINEAR_ORNS,nl_type=NONLINEAR_TYPE,\
        resultsdir='../results/odor_pulses'):
    ### read filename from the output file of automated run
    ## construct the filename
    if inh[0]: singles_str = '_NOSINGLES'
    else: singles_str = '_SINGLES'
    if inh[1]: joints_str = '_NOJOINTS'
    else: joints_str = '_JOINTS'
    if inh[3]: pgs_str = '_NOPGS'
    else: pgs_str = '_PGS'
    if inh[2]: lat_str = '_NOLAT'
    else: lat_str = '_LAT'
    if inh[4]: varmitstr = '_VARMIT'
    else: varmitstr = '_NOVARMIT'
    ## stable enough that time tags are not needed
    filename = resultsdir+'/scaledpulses_width'+str(_scaledWidth)+\
        '_netseed'+str(netseed)+'_stimseed'+str(stimseed)
    if nl_orns: filename += '_NL'+nl_type
    filename += singles_str+joints_str+pgs_str+lat_str+varmitstr+\
        '_numgloms'+str(num_gloms_list[ngi])
    if directed: filename += '_directed'+str(FRAC_DIRECTED)
    filename += '.pickle'
    return filename, (singles_str, joints_str, pgs_str, lat_str, varmitstr)

def plot_lin_contribs_oldpaperfigure():
    """ plot linearity scores for all options
    """
    fig1 = figure(figsize=(columnwidth,linfig_height),dpi=300,facecolor='w')
    ax1a = plt.subplot2grid((2,5),(0,0)) # full
    ax1b = plt.subplot2grid((2,5),(1,0)) # full
    ax2a = plt.subplot2grid((2,5),(0,1)) # no lat
    ax2b = plt.subplot2grid((2,5),(1,1)) # no lat
    ax3a = plt.subplot2grid((2,5),(0,2)) # no self
    ax3b = plt.subplot2grid((2,5),(1,2)) # no self
    ax4a = plt.subplot2grid((2,5),(0,3)) # no inh
    ax4b = plt.subplot2grid((2,5),(1,3)) # no inh
    ax5a = plt.subplot2grid((2,5),(0,4)) # non-lin
    ax5b = plt.subplot2grid((2,5),(1,4))
    ## inh = (no_singles,no_joints,no_lat,no_PGs,varyRMP)
    inh_options = [
        ('',0,(False,False,False,False,False),'all',False,None,(ax1a,ax1b)), \
        ('',1,(False,False,True,False,False),'no lat-inh',False,None,(ax2a,ax2b)), \
        ('',2,(True,False,False,True,False),'no self-inh',False,None,(ax3a,ax3b)), \
        ('',3,(True,True,False,True,False),'no inh',False,None,(ax4a,ax4b)), \
        ('',4,(False,False,False,False,False),'non-lin ORNs',True,'P',(ax5a,ax5b)) ]
    Rsqs_all = [ [] for i in range(len(inh_options)) ]
    R_all = [ [] for i in range(len(inh_options)) ]
    maxy_hist = 0.0
    for ploti,(dirextn,inhi,inh,inhstr,nl_orns,nl_type,(axa,axb)) in enumerate(inh_options):
        slopes_all = []
        peaks_all = []
        avg_frates_all = []
        for stimi,stimseed in enumerate(stim_seeds):
            filename, switch_strs \
                = get_filename(stimseed,stimseed,inh,0,nl_orns,nl_type,\
                    resultsdir='../results/odor_pulses'+dirextn)
            ## if the result file for these seeds & tweaks doesn't exist,
            ## then carry on to the next.
            if not os.path.exists(filename): continue
            worker = fit_sp.fit_plot_scaledpulses(filename,stimseed,_scaledWidth)
            fits_2mits,peaks_2mits = worker.fit_pulses(dirextn,True,False,False) # (dirextn,noshow,savefig,test)
            for mitrali in [0,1]:
                _,slopes,intercepts,r_values,_,_,se_slope,avg_frates = zip(*fits_2mits[mitrali])
                ## leave out the fit for reference with itself
                r_values = delete(r_values,fit_sp.ref_response_scalenum-1)
                R_all[ploti].extend(r_values)
                Rsqs_all[ploti].extend(r_values**2)
                nan_present = False
                for slope in slopes:
                    if isnan(slope): nan_present = True
                intercepts_bad = False
                #for intercept in intercepts:
                #    if abs(intercept)>2.0: intercepts_bad = True
                if not nan_present and not intercepts_bad:
                    slopes_all.append(slopes)
                    peaks_all.append(peaks_2mits[mitrali])
                    print filename
                else:
                    print 'Left out :',filename
                avg_frates_all.append(avg_frates)
        slopes_mean = append([0],mean(slopes_all,axis=0))
        slopes_std = append([0],std(slopes_all,axis=0))
        peaks_mean = append([0],mean(peaks_all,axis=0))
        peaks_std = append([0],std(peaks_all,axis=0))
        avg_frates_mean = mean(avg_frates_all,axis=0)
        print "The average firing rates for the different scaling is",avg_frates_mean
        print "Peaks mean is",peaks_mean
        
        ## plots
        axa.hist(R_all[ploti],10,range=(0,1.0),normed=True,histtype='step',\
            linewidth=linewidth,color='k',ls='solid')
        xmin,xmax,ymin,ymax = \
            beautify_plot(axa,x0min=True,y0min=True,xticksposn='bottom',yticksposn='left',yticks=[])
        axa.set_ylim(0,ymax) # later this is set to the max over all histograms maxy_hist
        maxy_hist = max(maxy_hist,ymax)
        ## x-axis label for full first row
        if inhi==2:
            axes_labels(axa,"correlation with 1x","",xpad=2)
            #axa.xaxis.set_label_coords(1.3,-0.25)
        ### inset plots within subplots
        ### passing transform=ax.transAxes to add_axes() doesn't work, hence jugglery from
        ### http://matplotlib.1069221.n5.nabble.com/Adding-custom-axes-within-a-subplot-td20316.html
        #Bbox = matplotlib.transforms.Bbox.from_bounds(1.6, 0.4, 0.6, 0.6) 
        #trans = ax.transAxes + fig1.transFigure.inverted() 
        #l, b, w, h = matplotlib.transforms.TransformedBbox(Bbox,trans).bounds
        #axinset = fig1.add_axes([l, b, w, h])
        xconcs = array(scaledList[1:])
        ref_scale = scaledList[fit_sp.ref_response_scalenum]
        print "Slopes are =",slopes_mean*ref_scale
        #axb.plot(range(6),range(6),color=(0,0,0.7,0.5),dashes=(2.0,1.0)) # linear reference
        axb.errorbar(x=append([0],xconcs),y=peaks_mean,yerr=peaks_std,\
            color='k',linewidth=linewidth,capsize=cap_size)
        #for slopes in slopes_all:
        #    axb.plot(xconcs,array(slopes)*ref_scale,linewidth=linewidth)
        _,_,_,ymax = beautify_plot(axb,x0min=True,y0min=True,\
            xticksposn='bottom',yticksposn='left',yticks=[])
        axb.set_xlim(0,5)
        axb.set_ylim(0,180)
        axb.set_xticks([0,1,5])
        axb.set_yticks([])
        ### Draw the twin y axis (turned off always by beautify_plot)
        #for loc, spine in axb.spines.items(): # items() returns [(key,value),...]
        #    spine.set_linewidth(axes_linewidth)
        #    if loc in ['right']:
        #        spine.set_color('k') # draw spine in black

        ## x-axis label for full second row
        if inhi==2:
            axes_labels(axb,"ORN scaling","",xpad=2)
            #axb.xaxis.set_label_coords(1.3,-0.25)
        if inhi==0: # laebl twin y axes of rightmost plot
            axb.set_yticks([0,180])
            axes_labels(axb,'','peak (Hz)',ypad=-6) # to set font size for twin y ticklabels

    for i,(_,_,_,_,_,_,(axa,axb)) in enumerate(inh_options):
        axa.set_ylim(0,maxy_hist)
        if i==0: # label y axis for left-most plots
            axa.set_yticks([0,maxy_hist])
            axes_labels(axa,'','density',ypad=2)
    fig1.tight_layout()
    fig_clip_off(fig1)
    fig1.subplots_adjust(top=0.95,left=0.1,bottom=0.15,right=0.98,wspace=0.4,hspace=0.4)
    #fig1.text(0.31,0.65,'density',fontsize=label_fontsize,\
    #        rotation='vertical', transform=fig1.transFigure)
    #fig1.text(0.6,0.025,'$R^2$',fontsize=label_fontsize,transform=fig1.transFigure)
    fig1.savefig('../figures/lin_contribs_scaledpulses.svg',dpi=fig1.dpi)
    fig1.savefig('../figures/lin_contribs_scaledpulses.png',dpi=fig1.dpi)

def plot_lin_contribs_paperfigure():
    """ plot linearity scores for all options
    """
    fig1 = figure(figsize=(columnwidth,linfig_height),dpi=300,facecolor='w')
    ax1a = plt.subplot2grid((2,5),(0,0)) # full
    ax1b = plt.subplot2grid((2,5),(1,0)) # full
    ax2a = plt.subplot2grid((2,5),(0,1)) # no lat
    ax2b = plt.subplot2grid((2,5),(1,1)) # no lat
    ax3a = plt.subplot2grid((2,5),(0,2)) # no self
    ax3b = plt.subplot2grid((2,5),(1,2)) # no self
    ax4a = plt.subplot2grid((2,5),(0,3)) # no inh
    ax4b = plt.subplot2grid((2,5),(1,3)) # no inh
    ax5a = plt.subplot2grid((2,5),(0,4)) # non-lin
    ax5b = plt.subplot2grid((2,5),(1,4))
    ## inh = (no_singles,no_joints,no_lat,no_PGs,varyRMP)
    inh_options = [
        ('',0,(False,False,False,False,False),'all',False,None,(ax1a,ax1b)), \
        ('',1,(False,False,True,False,False),'no lat-inh',False,None,(ax2a,ax2b)), \
        ('',2,(True,False,False,True,False),'no self-inh',False,None,(ax3a,ax3b)), \
        ('',3,(True,True,False,True,False),'no inh',False,None,(ax4a,ax4b)), \
        ('',4,(False,False,False,False,False),'non-lin ORNs',True,'P',(ax5a,ax5b)) ]
    Rsqs_all = [ [] for i in range(len(inh_options)) ]
    R_all = [ [] for i in range(len(inh_options)) ]
    maxy_R = 0.0
    for ploti,(dirextn,inhi,inh,inhstr,nl_orns,nl_type,(axa,axb)) in enumerate(inh_options):
        slopes_all = []
        peaks_all = []
        avg_frates_all = []
        R_thisinh = []
        for stimi,stimseed in enumerate(stim_seeds):
            filename, switch_strs \
                = get_filename(stimseed,stimseed,inh,0,nl_orns,nl_type,\
                    resultsdir='../results/odor_pulses'+dirextn)
            ## if the result file for these seeds & tweaks doesn't exist,
            ## then carry on to the next.
            if not os.path.exists(filename): continue
            worker = fit_sp.fit_plot_scaledpulses(filename,stimseed,_scaledWidth)
            fits_2mits,peaks_2mits = worker.fit_pulses(dirextn,True,False,False) # (dirextn,noshow,savefig,test)
            for mitrali in [0,1]:
                _,slopes,intercepts,r_values,_,_,se_slope,avg_frates = zip(*fits_2mits[mitrali])
                R_thisinh.append(r_values)
                ## leave out the fit for reference with itself
                r_values = delete(r_values,fit_sp.ref_response_scalenum-1)
                R_all[ploti].extend(r_values)
                Rsqs_all[ploti].extend(r_values**2)
                nan_present = False
                for slope in slopes:
                    if isnan(slope): nan_present = True
                intercepts_bad = False
                #for intercept in intercepts:
                #    if abs(intercept)>2.0: intercepts_bad = True
                if not nan_present and not intercepts_bad:
                    slopes_all.append(slopes)
                    peaks_all.append(peaks_2mits[mitrali])
                    print filename
                else:
                    print 'Left out :',filename
                avg_frates_all.append(avg_frates)
        slopes_mean = append([0],mean(slopes_all,axis=0))
        slopes_std = append([0],std(slopes_all,axis=0))
        peaks_mean = append([0],mean(peaks_all,axis=0))
        peaks_std = append([0],std(peaks_all,axis=0))
        avg_frates_mean = mean(avg_frates_all,axis=0)
        print "The average firing rates for the different scaling is",avg_frates_mean
        print "Peaks mean is",peaks_mean
        
        ## plots
        R_mean = mean(R_thisinh,axis=0)
        R_std = std(R_thisinh,axis=0)
        xconcs = scaledList[1:]
        axa.errorbar(xconcs,R_mean,R_std,linewidth=linewidth,color='b',\
            ls='solid',marker='o',ms=marker_size,capsize=cap_size)
        #axa.plot(xconcs,R_std,linewidth=linewidth,color='r',ls='solid',marker='x',ms=marker_size)
        xmin,xmax,ymin,ymax = \
            beautify_plot(axa,x0min=True,y0min=True,xticksposn='bottom',yticksposn='left',yticks=[])
        axa.set_ylim(0,ymax) # later this is set to the max over all histograms maxy_hist
        axa.set_xticks([0,1,2,5])
        maxy_R = max(maxy_R,ymax)
        ## x-axis label for full first row
        if inhi==2:
            axes_labels(axa,"","",xpad=2) # ORN scaling is on the next row x-axis
            #axa.xaxis.set_label_coords(1.3,-0.25)
        ### inset plots within subplots
        ### passing transform=ax.transAxes to add_axes() doesn't work, hence jugglery from
        ### http://matplotlib.1069221.n5.nabble.com/Adding-custom-axes-within-a-subplot-td20316.html
        #Bbox = matplotlib.transforms.Bbox.from_bounds(1.6, 0.4, 0.6, 0.6) 
        #trans = ax.transAxes + fig1.transFigure.inverted() 
        #l, b, w, h = matplotlib.transforms.TransformedBbox(Bbox,trans).bounds
        #axinset = fig1.add_axes([l, b, w, h])
        xconcs = array(scaledList[1:])
        ref_scale = scaledList[fit_sp.ref_response_scalenum]
        print "Slopes are =",slopes_mean*ref_scale
        #axb.plot(range(6),range(6),color=(0,0,0.7,0.5),dashes=(2.0,1.0)) # linear reference
        axb.errorbar(x=append([0],xconcs),y=peaks_mean,yerr=peaks_std,\
            color='b',linewidth=linewidth,capsize=cap_size)
        #for slopes in slopes_all:
        #    axb.plot(xconcs,array(slopes)*ref_scale,linewidth=linewidth)
        _,_,_,ymax = beautify_plot(axb,x0min=True,y0min=True,\
            xticksposn='bottom',yticksposn='left',yticks=[])
        axb.set_xlim(0,5)
        axb.set_ylim(0,180)
        axb.set_xticks([0,1,2,5])
        axb.set_yticks([])
        ### Draw the twin y axis (turned off always by beautify_plot)
        #for loc, spine in axb.spines.items(): # items() returns [(key,value),...]
        #    spine.set_linewidth(axes_linewidth)
        #    if loc in ['right']:
        #        spine.set_color('k') # draw spine in black

        ## x-axis label for full second row
        if inhi==2:
            axes_labels(axb,"conc (% SV)","",xpad=2)
            #axb.xaxis.set_label_coords(1.3,-0.25)
        if inhi==0: # laebl twin y axes of rightmost plot
            axb.set_yticks([0,180])
            axes_labels(axb,'','peak (Hz)',ypad=-6) # to set font size for twin y ticklabels

    for i,(_,_,_,_,_,_,(axa,axb)) in enumerate(inh_options):
        axa.set_ylim(0,1)
        if i==0: # label y axis for left-most plots
            axa.set_yticks([0,1])
            axes_labels(axa,'','corr to 1%',ypad=2)
    fig1.tight_layout()
    fig_clip_off(fig1)
    fig1.subplots_adjust(top=0.95,left=0.1,bottom=0.15,right=0.98,wspace=0.4,hspace=0.4)
    #fig1.text(0.31,0.65,'density',fontsize=label_fontsize,\
    #        rotation='vertical', transform=fig1.transFigure)
    #fig1.text(0.6,0.025,'$R^2$',fontsize=label_fontsize,transform=fig1.transFigure)
    fig1.savefig('../figures/lin_contribs_scaledpulses.svg',dpi=fig1.dpi)
    fig1.savefig('../figures/lin_contribs_scaledpulses.png',dpi=fig1.dpi)

def plot_responseshapes_vs_conc():
    fig = figure(figsize=(columnwidth/3.*2,linfig_height/2.),dpi=300,facecolor='w')
    ax1 = fig.add_subplot(1,2,1)
    ax2 = fig.add_subplot(1,2,2)
    xconcs = array(scaledList[1:])
    ## inh = (no_singles,no_joints,no_lat,no_PGs,varyRMP)
    inh_option = ('',0,(False,False,False,False,False),'all',False,None)
    tpeaks_all = []
    tpeaks_normed_all = []
    durations_normed_all = []
    risetimes_all = []
    durations_all = []
    deltapeaks_all = []
    latencies_all = []
    latencies_normed_all = []
    (dirextn,inhi,inh,inhstr,nl_orns,nl_type) = inh_option
    for stimi,stimseed in enumerate(stim_seeds):
        filename, switch_strs \
            = get_filename(stimseed,stimseed,inh,0,nl_orns,nl_type,\
                resultsdir='../results/odor_pulses'+dirextn)
        ## if the result file for these seeds & tweaks doesn't exist,
        ## then carry on to the next.
        if not os.path.exists(filename): continue
        print filename
        ######### 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(filename)
        ##-------------------------- 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,fit_sp.pulserebindt,SCALED_RUNTIME)
        for mitrali in [0,1]:
            ## take the odor responses of the mitral to be fitted
            firingbinsmeanList = mitral_responses_mean[:,mitrali]
            
            starti = int(PULSE_START/fit_sp.pulserebindt)
            endi = int((PULSE_START+scaledWidth+kernel_time)/fit_sp.pulserebindt)
            air_bgnd = firingbinsmeanList[0]
            air_bgnd_relevant = firingbinsmeanList[0][starti:endi]
            tpeaks = []
            risetimes = []
            durations = []
            deltapeaks = []
            latencies = []
            for scalenum in [1,2,3,4,5]: ## conc scaled pulses
                ## find the mean first spike latency after odor onset
                latency = 0
                for avgi in range(len(mitral_responses_list)):
                    for spiketime in mitral_responses_list[avgi][scalenum][mitrali]:
                        if spiketime > PULSE_START:
                            latency += spiketime - PULSE_START
                            break
                latency_mean = latency / float(len(mitral_responses_list))
                latencies.append(latency_mean)
                ## find peak times and durations of binned responses
                scale = scaledList[scalenum]
                response = firingbinsmeanList[scalenum][starti:endi]-air_bgnd_relevant
                tpeak = argmax(response)
                tpeaks.append(tpeak)
                ##--------- Calculate half-max duration of the response
                ## f is an interpolated 1D function that gives
                ## down-shifted response that is zero at halfmax for fsolve/bisect below
                f = interpolate.interp1d(range(len(response)),response-max(response)/2.,\
                        bounds_error=False,fill_value=0.0)
                ### fsolve to find half-max points on either side of tpeak
                ### but often gives negative values!! so ditch this
                #halfresp_low = optimize.fsolve(f,tpeak-1)[0]
                #halfresp_high = optimize.fsolve(f,tpeak+1)[0]
                ## use bisect to find halfmax pt within an interval
                ## but bisect need opposite signs at the ends of the interval
                ## else ValueError, so increase interval, till you get it
                tmin = tpeak-1
                while True:
                    try:
                        halfresp_low = optimize.bisect(f,tpeak,tmin-1)
                        break
                    except ValueError:
                        tmin -= 1
                tmax = tpeak+1
                while True:
                    try:
                        halfresp_high = optimize.bisect(f,tpeak,tmax)
                        break
                    except ValueError:
                        tmax += 1
                risetimes.append(halfresp_low)
                duration = (halfresp_high-halfresp_low)
                durations.append(duration)
                
            ## take only the reduction in the time to peak vs concentration
            tpeaks = array(tpeaks)*fit_sp.pulserebindt*1000
            risetimes = array(risetimes)*fit_sp.pulserebindt*1000
            durations = array(durations) *fit_sp.pulserebindt*1000
            tpeaks_all.append(tpeaks)
            tpeaks_normed_all.append(tpeaks/mean(tpeaks))
            risetimes_all.append(risetimes)
            durations_all.append(durations)
            durations_normed_all.append(durations/mean(durations))
            deltapeaks = [tpeaks[-1]-tpeaks[-2],tpeaks[-2]-tpeaks[-3]]
            deltapeaks_all.extend(deltapeaks)
            latencies_all.append(latencies)
            latencies_normed_all.append(latencies/mean(latencies))
            #ax1.plot(xconcs,tpeaks)
            #ax2.plot(xconcs,durations)
    
    tpeaks_mean = mean(tpeaks_all,axis=0)
    tpeaks_std = std(tpeaks_all,axis=0)
    tpeaks_normed_mean = mean(tpeaks_normed_all,axis=0)
    tpeaks_normed_se = std(tpeaks_normed_all,axis=0) / sqrt(len(tpeaks_normed_all))
    risetimes_mean = mean(risetimes_all,axis=0)
    risetimes_std = std(risetimes_all,axis=0)
    durations_mean = mean(durations_all,axis=0)
    durations_std = std(durations_all,axis=0)
    durations_normed_mean = mean(durations_normed_all,axis=0)
    durations_normed_se = std(durations_normed_all,axis=0) / sqrt(len(durations_normed_all))
    latencies_mean = mean(latencies_all,axis=0)
    latencies_std = std(latencies_all,axis=0)
    latencies_normed_mean = mean(latencies_normed_all,axis=0)
    latencies_normed_se = std(latencies_normed_all,axis=0) / sqrt(len(latencies_normed_all))
    ## plots
    ax1.errorbar(xconcs,latencies_normed_mean,latencies_normed_se,color='b',linewidth=linewidth)
    #ax1.errorbar(xconcs,risetimes_mean,risetimes_std,color='b',linewidth=linewidth)
    ax2.errorbar(xconcs,tpeaks_normed_mean,tpeaks_normed_se,color='b',linewidth=linewidth)
    #ax2.errorbar(xconcs,durations_normed_mean,durations_normed_se,color='b',linewidth=linewidth)
    for i,ax in enumerate([ax1,ax2]):
        _,_,ymin,ymax = beautify_plot(ax,x0min=True,y0min=False,\
            xticksposn='bottom',yticksposn='left')
        ax.set_xlim(0,5)
        ax.set_xticks([0,1,2,5])
        ax.set_yticks([ymin,1,ymax])
        axes_labels(ax,'conc (% SV)',['norm-t-spike','norm-t-peak'][i])
    ### Draw the twin y axis (turned off always by beautify_plot)
    #for loc, spine in axb.spines.items(): # items() returns [(key,value),...]
    #    spine.set_linewidth(axes_linewidth)
    #    if loc in ['right']:
    #        spine.set_color('k') # draw spine in black

    fig_clip_off(fig)
    fig.tight_layout()
    fig.savefig('../figures/scaledpulses_latency.svg',dpi=fig.dpi)
    fig.savefig('../figures/scaledpulses_latency.png',dpi=fig.dpi)
    
    fig = figure()
    ax = fig.add_subplot(111)
    ax.hist(deltapeaks_all)
    axes_labels(ax,'delta peak times','#')

if __name__ == "__main__":
    ## plot linearity with parts of the network removed / tweaked
    #plot_lin_contribs_oldpaperfigure()
    ## PAPER figure. Similar to above, but plotting mean and sd of corr with 1x
    ## versus concentration instead of a combined distrib.
    plot_lin_contribs_paperfigure()
    ## Compare pulse shapes versus concentration
    plot_responseshapes_vs_conc()
    
    show()

Loading data, please wait...