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;
#!/usr/bin/env python
# -*- coding: utf-8 -*-

import os
import sys
import math
import pickle

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

from moose.utils import * # imports moose

from stimuliConstants import * # has SETTLETIME
from simset_activinhibition import * # has REALRUNTIME
from data_utils import *
from plot_mitral_connectivity_NetworkML import * # get connectivity from netfile

RUNTIME = REALRUNTIME + SETTLETIME

from pylab import * # part of matplotlib that depends on numpy but not scipy

ASYM_TEST = False
DIRECTED = True
FRAC_DIRECTED = 0.01#0.003#0.03
IN_VIVO = False#True
if IN_VIVO: invivo_str = '_invivo'
else: invivo_str = ''
if DIRECTED: dirt_str = '_directed'+str(FRAC_DIRECTED)
else: dirt_str = ''

USE_BEST_HALF = True

## inh =  (no_singles,no_joints,no_PGs)
inh = (False,False,False)
netseeds = arange(100.0,5000.0,100.0)

if IN_VIVO:
    distance_list = [50.0, 200.0, 400.0, 600.0]#, 800.0] # microns
    reverse_list = [True,False]
else:
    distance_list = [50.0] # microns
    reverse_list = [False]
## whether to calculate and print the number of joint granules connected
SHOW_CONN = False

########## You need to run: python2.6 average_activdep_inhibition.py

## Read in the .csv file having the average ADI of Arevian et al.
## numpy's genfromtxt() to read csv files
## Arevian_mean_ADI['x'], Arevian_mean_ADI['Curve1'] contain x and y lists
Arevian_mean_ADI = genfromtxt('arevianetal_average_ADI.csv', delimiter=',', dtype=None, names=True)


def get_filename(netseed,inh,mitdistance,reverse,latmitnum=2):
    ### 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[2]: pgs_str = '_NOPGS'
    else: pgs_str = '_PGS'
    if reverse: rev_str = '_reversed'
    else: rev_str = ''
    if ASYM_TEST: asym_str = '_asym'
    else: asym_str = ''
    ## stable enough that time tags are not needed
    outfilename = '../results/ADI/ADI_'+str(latmitnum)+\
        '_seed'+str(netseed)+'_mitdist'+str(mitdistance)+\
        singles_str+joints_str+pgs_str+invivo_str+dirt_str+\
        rev_str+asym_str+'.pickle'
    return outfilename, (singles_str, joints_str, pgs_str)

def get_rates_inh(filename):
    f = open(filename,'r')
    Ainjectarray, both_firingratearrays = pickle.load(f)
    f.close()
    mit_alone = array(both_firingratearrays[0])
    mit_inhibited = array(both_firingratearrays[1])
    only_reasonable_inh = where(Ainjectarray<1.5e-9)
    diff_frates = mit_alone[only_reasonable_inh] - mit_inhibited[only_reasonable_inh]
    max_redux = max(diff_frates)
    avg_redux = mean(diff_frates)
    ## take max f rate redux from 6 to 55Hz
    only_reasonable_55Hz = where(Ainjectarray<1.2e-9)
    mit_alone_55Hz = mit_alone[only_reasonable_55Hz]
    only_to_55Hz = where(mit_alone_55Hz<=55.0)
    only_6Hz_to_55Hz = where(mit_alone_55Hz[only_to_55Hz]>=6.0)
    diff_frates = mit_alone_55Hz[only_6Hz_to_55Hz] - mit_inhibited[only_6Hz_to_55Hz]
    max_asym_redux = max(diff_frates)
    return Ainjectarray, both_firingratearrays, max_redux, avg_redux, max_asym_redux

def plot_avg_inh(reverse, mit_distance, title_str):
    ## load the data
    full_firingratearray = []
    redux_frate_array = []
    ## 0 to 140 Hz into 10 Hz bins
    maxfrate = 140
    fratebinsize = 10
    redux_vs_frate_array = [ [] for i in range(maxfrate/fratebinsize) ]
    max_redux_frate_array = []
    avg_redux_frate_array = []
    max_asym_redux_frate_array = []
    mit_somadend_joints_avg = array([0,0])
    num_available = 0
    for netseed in netseeds:
        filename,_ = get_filename(netseed,inh,mit_distance,reverse)
        ## if the result file for these seeds & tweaks doesn't exist,
        ## then carry on to the next.
        print 'checking for',filename
        if not os.path.exists(filename): continue
        num_available += 1
        print 'found',filename

        Ainjectarray, both_firingratearrays, max_redux, avg_redux, max_asym_redux = \
            get_rates_inh(filename)
        full_firingratearray.append(both_firingratearrays)
        redux_frate_array_this = array(both_firingratearrays[1])-array(both_firingratearrays[0])
        redux_frate_array.append(redux_frate_array_this)
        max_redux_frate_array.append(max_redux)
        avg_redux_frate_array.append(avg_redux)
        max_asym_redux_frate_array.append(max_asym_redux)
        
        if SHOW_CONN:
            ## get connectivity details about this netfile
            bits = filename.split('_')
            netseedstr = bits[1]
            mitdistancestr = bits[2]
            OBNet_file = '../netfiles/syn_conn_array_10000_singlesclubbed20_jointsclubbed1'\
                '_numgloms2_'+netseedstr+'_'+mitdistancestr
            if directed: OBNet_file += '_directed0.0_proximal'
            OBNet_file += '_2GLOMS'
            if not IN_VIVO: OBNet_file += '_INVITRO.xml'
            else: OBNet_file += '.xml'
            netml = NetworkML(twogloms=True)
            tweaks = {}
            ## mit_posn_grannos = 
            ## (mit_prim_joints,mit_somadend_joints,mit_sec_joints,mit_prim_singles,mit_sec_singles)
            (mitralsDict,connectivityDict,mit_posn_grannos) = \
                netml.readNetworkMLFromFile(OBNet_file,params=tweaks)
            mit_somadend_joints = mit_posn_grannos[1] # dict {0:0,mitB:0}, mitB=2 if TWOGLOMS
            mit_somadend_joints_avg[0] += mit_somadend_joints[0]
            mit_somadend_joints_avg[1] += mit_somadend_joints[2]

    ## select 50% of the strongest ADI : may or may not use for plots: see below
    ## Currently not using only strongest.
    ## get indices of the best ADI, sorted according to max redux in frate, to keep them
    low_idx = [i[0] for i in sorted(enumerate(max_redux_frate_array), key=lambda x:x[1], reverse=True)]
    low_idx = low_idx[0:num_available/2]
    keep_arrays = []
    keep_redux_frate_array = []
    for idx in low_idx:
        keep_arrays.append(full_firingratearray[idx])
        keep_redux_frate_array.append(redux_frate_array[idx])
    keep_arrays = array(keep_arrays)
    full_firingratearray = array(full_firingratearray)

    ## calculate means and std errs
    ## comment/uncomment one of the two blocks below to take mean of top 50% OR all
    ## mean inh trough of 50% is approx 25% lower than that of all
    ## mean of 50%
    if USE_BEST_HALF:
        print "Using means of top 50% of the mitral pairs having most inh redux."
        mean_firingratearray = mean(keep_arrays, axis=0)
        std_firingratearray = std(keep_arrays, axis=0)
        redux_frate_array = keep_redux_frate_array
        nummitpairs = len(keep_arrays)
    else:
    ## mean of all
        print "Using means of all mitral pairs, rather than top half."
        mean_firingratearray = mean(full_firingratearray, axis=0)
        std_firingratearray = std(full_firingratearray, axis=0)
        nummitpairs = len(full_firingratearray)

    print "Total number of mitral pairs used for means =",nummitpairs
    ## bin the reduction so that you get reduction vs firing rate of A (rather than I_A)
    ## some bins will have more elements than others
    for redux_frate_array_this in redux_frate_array:
        for i,frateA in enumerate(both_firingratearrays[0]):
            redux_vs_frate_array[int(frateA/fratebinsize)].append(redux_frate_array_this[i])

    if SHOW_CONN:
        mit_somadend_joints_avg /= float(len(filelist))
        print "Avg # of granules to soma and nearest segments to mits A and B are",\
            mit_somadend_joints_avg

    ## plot the frate redux and percentage redux in firing rates
    frateA_array = []
    frateA_running = fratebinsize/2.0
    frate_redux_mean_array = []
    frate_redux_std_array = []
    for arr in redux_vs_frate_array:
        if arr:
            frate_redux_mean_array.append( mean(arr) )
            frate_redux_std_array.append( std(arr) )
            frateA_array.append(frateA_running)
        frateA_running += fratebinsize
    redux_perc_mean_array = array(frate_redux_mean_array)/array(frateA_array)*100
    redux_perc_std_array = array(frate_redux_std_array)/array(frateA_array)*100
    
    ## simulation of fig 2e of Arevian et al
    fig = figure(figsize=(columnwidth*3/4, linfig_height/2.0), dpi=fig_dpi, facecolor='w')
    ### plot reduction in percentage
    #ax2 = fig.add_subplot(111,frameon=False)
    #ax2.fill_between(frateA_array, \
    #    redux_perc_mean_array + redux_perc_std_array/sqrt(nummitpairs), \
    #    redux_perc_mean_array - redux_perc_std_array/sqrt(nummitpairs), \
    #    color='b',alpha=0.4,linewidth=0)
    #ax2.plot(frateA_array, redux_perc_mean_array, \
    #    color='b',marker='s',ms=marker_size,linewidth=linewidth)
    #ax2.set_ylim(-30,0)
    #ax2.set_yticks([-30,-20,-10,0])
    #ax2.add_artist(Line2D((0, 0), (-30, 0), \
    #    color='black', linewidth=axes_linewidth))
    #ax2.get_xaxis().set_ticks_position('bottom')
    #axes_labels(ax2,"mitral A firing rate (Hz)","change in rate of A (%)")
    
    ## plot delta rate in Hz
    #ax = ax2.twinx()
    ax = fig.add_subplot(111)
    ax.plot(Arevian_mean_ADI['x'],Arevian_mean_ADI['Curve1'],\
        color='b',dashes=(1.5,0.5),linewidth=linewidth)
    ax.fill_between(frateA_array, \
        array(frate_redux_mean_array) + array(frate_redux_std_array)/sqrt(nummitpairs), \
        array(frate_redux_mean_array) - array(frate_redux_std_array)/sqrt(nummitpairs), \
        color='r',alpha=0.4,linewidth=0)
    ax.plot(frateA_array, frate_redux_mean_array, \
        color='r',marker='o',ms=marker_size,linewidth=linewidth)
    ax.set_xlim(0,maxfrate)
    ax.set_xticks([0,maxfrate/4,maxfrate/2,maxfrate*3/4,maxfrate])
    ax.set_ylim(-15,0)
    ax.set_yticks([-15,-10,-5,0])
    axes_labels(ax,"rate of A (Hz)","$\Delta$ rate of A (Hz)")
    ax.xaxis.set_label_position('top')
    ax.get_yaxis().set_ticks_position('left')
    ax.get_xaxis().set_ticks_position('top')
    ## set xticks text sizes
    for label in ax.get_xticklabels():
        label.set_fontsize(label_fontsize)
    ## hide the top and left axes
    for loc, spine in ax.spines.items(): # items() returns [(key,value),...]
        spine.set_linewidth(axes_linewidth)
        if loc in ['right','bottom']:
            spine.set_color('none') # don't draw spine
    fig.tight_layout()
    fig.subplots_adjust(right=0.8)
    fig.savefig('../figures/ADI/ADI'+invivo_str+dirt_str+'.svg',dpi=fig.dpi)
    
    ## Plot the best ADI-s separately also:
    for ploti,idx in enumerate(low_idx):
        mainfig = figure(figsize=(columnwidth/2.0,linfig_height),dpi=300,facecolor='w')
        mainaxes = mainfig.add_subplot(111,frameon=False)
        best_firingratearray = full_firingratearray[idx]
        for mitralBinject in [0,1]:
            if ASYM_TEST:
                if mitralBinject == 0: mitB_str = '0.0 A'
                else: mitB_str = 'A inject'
            else: mitB_str = str(mitralBinject*onInject)+' A'
            mainaxes.plot(Ainjectarray*1e9, best_firingratearray[mitralBinject],\
            color=(mitralBinject,0,0), marker=['s','o'][mitralBinject],\
            label="mitral B inject = "+mitB_str, linewidth=linewidth, markersize=marker_size)
        mainaxes.get_yaxis().set_ticks_position('left')
        mainaxes.get_xaxis().set_ticks_position('bottom')
        mainaxes.set_xticks([0,1,2.5])
        mainaxes.set_ylim(0,120)
        mainaxes.set_yticks([0,20,40,60,80,100,120])
        axes_labels(mainaxes,"mitral A inject (nA)","mitral A firing rate (Hz)") # enlarges x and y ticks and labels
        xmin, xmax = mainaxes.get_xaxis().get_view_interval()
        ymin, ymax = mainaxes.get_yaxis().get_view_interval()
        mainaxes.add_artist(Line2D((0, 0), (0, ymax), color='black', linewidth=axes_linewidth))
        mainaxes.add_artist(Line2D((0, xmax), (0, 0), color='black', linewidth=axes_linewidth))
        fig_clip_off(mainfig)
        mainfig.tight_layout()
        mainfig.savefig('../figures/ADI/ADI' \
            +invivo_str+dirt_str+rev_str+'_'+str(idx)+'.png', \
            dpi=mainfig.dpi)
        mainfig.savefig('../figures/ADI/ADI' \
            +invivo_str+dirt_str+rev_str+'_'+str(idx)+'.svg', \
            dpi=mainfig.dpi)

    ## plot the mean ADI firing rate vs I curves with B on and off
    mainfig = figure(figsize=(columnwidth/2.0,linfig_height),dpi=300,facecolor='w')
    mainaxes = mainfig.add_subplot(111,frameon=False)
    for mitralBinject in [0,1]:
        if ASYM_TEST:
            if mitralBinject == 0: mitB_str = '0.0 A'
            else: mitB_str = 'A inject'
        else: mitB_str = str(mitralBinject*onInject)+' A'
        mainaxes.errorbar(x=Ainjectarray*1e9,\
        y=mean_firingratearray[mitralBinject],yerr=std_firingratearray[mitralBinject],\
        color=(mitralBinject,0,0), marker=['s','o'][mitralBinject], capsize=cap_size,\
        label="mitralB I = "+mitB_str, linewidth=linewidth, markersize=marker_size)
    #if IN_VIVO: biglegend("lower right")#"upper left")
    #else: biglegend("lower right")
    axes_labels(mainaxes,"mitral A inject (nA)","mitral A firing rate (Hz)")
    mainaxes.set_title('mean '+title_str,fontsize=label_fontsize)
    mainaxes.get_yaxis().set_ticks_position('left')
    mainaxes.get_xaxis().set_ticks_position('bottom')
    mainaxes.set_xticks([0,1,2.5])
    mainaxes.set_ylim(0,120)
    mainaxes.set_yticks([0,20,40,60,80,100,120])
    xmin, xmax = mainaxes.get_xaxis().get_view_interval()
    ymin, ymax = mainaxes.get_yaxis().get_view_interval()
    mainaxes.set_xlim(0,xmax)
    mainaxes.set_ylim(0,ymax)
    mainaxes.add_artist(Line2D((0, 0), (0, ymax), color='black', linewidth=axes_linewidth))
    mainaxes.add_artist(Line2D((0, xmax), (0, 0), color='black', linewidth=axes_linewidth))
    fig_clip_off(mainfig)
    mainfig.tight_layout()
    mainfig.savefig('../figures/ADI/ADI'+invivo_str+'_'+title_str+'.png',\
        bbox_inches='tight',dpi=mainfig.dpi)

    ## print the max reduction in firing rates
    print 'Max redux '+title_str+' till IA=1.5nA =',max_redux_frate_array
    print 'Avg redux '+title_str+' till IA=1.5nA =',avg_redux_frate_array
    print 'Max redux '+title_str+' 6Hz to 55Hz =',max_asym_redux_frate_array
    return max_redux_frate_array, avg_redux_frate_array, max_asym_redux_frate_array

if __name__ == "__main__":
    redux = []
    redux_std = []
    for i,reverse in enumerate(reverse_list):
        print "Reversed =", reverse
        if reverse: rev_str = "Reversed"
        else: rev_str = "Forward"
        redux_directed = []
        redux_directed_std = []
        for j,mit_distance in enumerate(distance_list):
            print 'Distance =',mit_distance,'microns:'
            max_redux_frate_array, avg_redux_frate_array, max_asym_redux_frate_array = \
                plot_avg_inh(reverse, mit_distance, rev_str+' @ '+str(mit_distance)+' microns')
            redux_directed.append(mean(avg_redux_frate_array))
            redux_directed_std.append(std(avg_redux_frate_array))
        redux.append((rev_str,redux_directed))
        redux_std.append(redux_directed_std)
    if IN_VIVO:
        fig = figure(figsize=(columnwidth,linfig_height),dpi=300,facecolor='w')
        ax = fig.add_subplot(111,frameon=False)
        errorbar(distance_list,redux[0][1],yerr=redux_std[0],\
            color=(1,0,0),linewidth=linewidth,marker='o',label=redux[0][0])
        errorbar(distance_list,redux[1][1],yerr=redux_std[1],\
            color=(0,0,1),linewidth=linewidth,marker='x',label=redux[1][0])
        biglegend('upper right')#'lower left')
        axes_labels(ax,'separation (microns)','mean frate redux (Hz)') # (0 to 1.5nA) in 400ms
        ax.set_xticks(distance_list)
        ax.get_yaxis().set_ticks_position('left')
        ax.get_xaxis().set_ticks_position('bottom')
        xmin, xmax = ax.get_xaxis().get_view_interval()
        ymin, ymax = ax.get_yaxis().get_view_interval()
        ax.set_ylim(0,ymax)
        ax.set_yticks([0,ymax])
        ax.add_artist(Line2D((0, 0), (0, ymax), color='black', linewidth=axes_linewidth))
        ax.add_artist(Line2D((0, xmax), (0, 0), color='black', linewidth=axes_linewidth))
        fig.tight_layout()
        fig.savefig('../figures/ADI/ADI_asym_invivo.png',bbox_inches='tight',dpi=fig.dpi)
    if ASYM_TEST:
        max_redux_frate_array_rev, max_asym_redux_frate_array_rev =  \
            plot_avg_inh(filelist_invitro_reversed, 'reversed ')
    
        mainfig = figure(facecolor='w')
        mainaxes = mainfig.add_subplot(111)
        for i,redux in enumerate(max_asym_redux_frate_array):
            mainaxes.plot([1,2],[-max_asym_redux_frate_array_rev[i],-redux],
                marker='+',linewidth='2.0')
        mainaxes.set_xticks([0.8,1,2,2.2])
        mainaxes.set_xticklabels(['','A','B',''])
        axes_labels(mainaxes,"mitral cell","reduction in mitral firing (Hz)")

    ## plot the fig 2c of Arevian et al.
    expdata_Boff = [(0,0),(79.429,0),(159.286,5.067),(239.143,9.966),(317.464,17.748),\
    (400.393,19.765),(481.786,30.141),(560.107,44.839),(639.964,59.826),(719.821,82.306),\
    (801.214,97.292),(879.536,109.973),(959.393,116.890),(1039.250,129.571),(1122.179,136.776)]
    expdata_Bon = [(0,0),(79.429,4.779),(160.821,7.661),(239.143,12.272),(322.071,17.460),\
    (398.857,25.241),(480.250,29.853),(560.107,39.940),(639.964,52.332),(719.821,65.013),\
    (801.214,84.899),(879.536,104.497),(959.393,119.196),(1040.786,127.265),(1120.643,139.658)]
    ## remove the two highest points
    expdata_Boff = expdata_Boff[:-2]
    expdata_Bon = expdata_Bon[:-2]
    fig = figure(figsize=(columnwidth/2.0,linfig_height),dpi=300,facecolor='none') # transparent
    ax = fig.add_subplot(111,frameon=False)
    expdata_BoffX,expdata_BoffY = zip(*expdata_Boff)
    ax.plot(array(expdata_BoffX)*1e-3,expdata_BoffY,\
        color='k', linewidth=linewidth, marker='s',markersize=marker_size)
    expdata_BonX,expdata_BonY = zip(*expdata_Bon)
    ax.plot(array(expdata_BonX)*1e-3,expdata_BonY,\
        color='r', linewidth=linewidth, marker='o',markersize=marker_size)
    axes_labels(ax,'mitral A inject (nA)','mitral A firing rate (Hz)')
    ax.get_yaxis().set_ticks_position('left')
    ax.get_xaxis().set_ticks_position('bottom')
    ax.set_ylim(0,120)
    ax.set_yticks([0,20,40,60,80,100,120])
    ax.set_xticks([0,0.5,1.2])
    xmin, xmax = ax.get_xaxis().get_view_interval()
    ymin, ymax = ax.get_yaxis().get_view_interval()
    ax.add_artist(Line2D((0, 0), (0, ymax), color='black', linewidth=axes_linewidth))
    ax.add_artist(Line2D((0, xmax), (0, 0), color='black', linewidth=axes_linewidth))
    fig_clip_off(fig)
    fig.tight_layout()
    fig.savefig('../figures/ADI/ADI_exp.png',dpi=fig.dpi,transparent=True)
    fig.savefig('../figures/ADI/ADI_exp.svg',dpi=fig.dpi,transparent=True)

    show()

Loading data, please wait...