Gamma-beta alternation in the olfactory bulb (David, Fourcaud-Trocmé et al., 2015)

 Download zip file 
Help downloading and running models
Accession:185014
This model, a simplified olfactory bulb network with mitral and granule cells, proposes a framework for two regimes of oscillation in the olfactory bulb: 1 - a weak inhibition regime (with no granule spike) where the network oscillates in the gamma (40-90Hz) band 2 - a strong inhibition regime (with granule spikes) where the network oscillates in the beta (15-30Hz) band. Slow modulations of sensory and centrifugal inputs, phase shifted by a quarter of cycle, possibly combined with short term depression of the mitral to granule AMPA synapse, allows the network to alternate between the two regimes as observed in anesthetized animals.
Reference:
1 . David F, Courtiol E, Buonviso N, Fourcaud-Trocmé N (2015) Competing Mechanisms of Gamma and Beta Oscillations in the Olfactory Bulb Based on Multimodal Inhibition of Mitral Cells Over a Respiratory Cycle. eNeuro [PubMed]
Citations  Citation Browser
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 cell; Olfactory bulb main interneuron granule MC cell;
Channel(s): I_Ks;
Gap Junctions:
Receptor(s): GabaA; AMPA;
Gene(s):
Transmitter(s):
Simulation Environment: Brian; Python;
Model Concept(s): Short-term Synaptic Plasticity; Gamma oscillations; Beta oscillations; Olfaction;
Implementer(s):
Search NeuronDB for information about:  Olfactory bulb main mitral cell; Olfactory bulb main interneuron granule MC cell; GabaA; AMPA; I_Ks;
# -*- coding: utf-8 -*-
"""
This script contains the main fucntion to simulate the network and the associated dictionary of default parameters.

The "main" section provide a way to simulate and record 
either a single run with lot of details 
or multiple runs (with multiprocessing) with only LFPs
"""

#~ import brian_no_units
from brian import *

from numpy.random import seed

from model_mitral_clean import *
from model_granule_clean import *
from populationstatemonitor import *

from oscillation_analysis import beta_gamma_detection, spectrum_analysis
from plot_single_run_from_file import single_LFP_analysis
from plot_multi_run_from_file import multi_LFP_analysis,LFP_frequency_content


# Dictionary of default parameters for simulations, used as a basis to generate dictionary of specific simulations
param_dict_default=dict(
                        # General simulation params
                        dt=0.05*ms,
                        sim_time=4.*second,
                        t_init_M=0.*ms, # time start of mitral sensory excitation (before, no input)
                        t_init_G=0.*ms, # time start of granule centrifugal modulation (before, base input is used)
                        
                        # Mitral intrinsic and input params
                        N_mitral = 100,
                        num_activated_mitral=100, # introduction of this parameter simplify the random selection of activated mitral cells at each run
                        M_gEinj= (linspace,(6.1,7.6,100),siemens*meter**-2), # (linspace(6.1,7.6,N_mitral))*siemens*meter**-2, # MC max sensory input conductances
                        M_gEinj_shift=0., # Used only to test gamma for distinct average M_gEinj
                        M_gEinj_base=None, # use None for no slow respiratory modulation of MC input, otherwise generally 4.0*siemens*meter**-2 
                        M_taumKs = 7.*ms, 
                        M_tauI=7.*ms, # weak inh decay time
                        M_tauI_G=7.*ms, # strong inh decay time
                        M_gI_cst= 20.*siemens*meter**-2, # mitral constant inhibitory input
                        recorded_mitrals=[10,20,30,40,50,60,70,80,90,99], # if return_details=True only !
                        
                        # Granule input params
                        use_granule=False,
                        N_granule=100,
                        G_input='constant', # choose between "constant", "ramp", "sinusoid"
                        G_I_base=-4.*nA, # before t_init or for "constant" input
                        G_I_min=-4.*nA, # min of sinusoid
                        G_I_max=-.1*nA, # max of sinusoid
                        recorded_granules=[10,20,30,40,50,60,70,80,90,99], # if return_details=True only !
                        
                        # Respiratory rhythm parameters
                        freq_modul=2.*Hz, 
                        MC_phase_dispersion=1.5, # SD of Gaussian distribution of phase shifts
                        GC_phase_dispersion=0.2, # SD of Gaussian distribution of phase shifts
                        GC_phase_shift=-pi/2, # Phase shift of granule centrifugal input vs mitral sensory input
                        
                        # Connectivity parameters
                        weakinh_connectproba=1.,
                        weakinh_weight=1.,
                        weakinh_mean_delay=9.*ms, # final delays: uniform in the range "mean +/- dispersion"
                        weakinh_delay_dispersion=4.*ms, # final delays: uniform in the range "mean +/- dispersion"
                        stronginh_connectproba=0.5,
                        stronginh_weight=1.,
                        exc_weight=1., # Automatically multiplied by 2.5 if use_AMPA_STP is set to True
                        use_AMPA_STP=False, # MC to GC short term plasticity (calibrated for a pure "depression")
                        )

def reseau_mitral_granule(param_dict,return_details=True,report='text'):
    """
    Main function to launch a single run of the network with a set of parameters given by param_dict
    Check param_dict_default for the list of parameters and their roles
    
    If return_details is set to True, a result_dict is produced with many details recorded during the simulations,
    otherwise only the LFP is returned
    
    "report" corresponds to the BRIAN "brian.run" function parameter
    """

    # Clear Brian objects in memory
    reinit_default_clock()
    clear(True)
    seed()

    # Parameters
    defaultclock.dt = param_dict['dt']
    sim_time = param_dict['sim_time']
    t_init_M = param_dict['t_init_M']
    t_init_G = param_dict['t_init_G']
    use_granule = param_dict['use_granule']
    
    # Mitral model and initialization
    N_mitral= param_dict['N_mitral']
    M = NeuronGroup(N_mitral, Mitral_eqs, threshold=-30*mV, reset=Mitral_reset,freeze=False,compile=True)
    M.V=(-60.+rand(N_mitral)*15.)*mV
    M.W=0.1*(1.+0.5*rand(N_mitral))
    M.X=0.3*(1.+0.5*rand(N_mitral))
    M.Y=0.
    func_name,func_param,func_unit=param_dict['M_gEinj']
    M.gEinj= func_name(*func_param)*func_unit+param_dict['M_gEinj_shift']*func_unit
    
    if param_dict['num_activated_mitral']<N_mitral:
        # Randomly put to 0 the activation of N_mitral-param_dict['num_activated_mitral']
        M.gEinj[np.random.permutation(N_mitral)[:100-int(param_dict['num_activated_mitral'])]]=0
    
    if param_dict['M_gEinj_base'] is None:
        gEinj_base=M.gEinj # no oscillation of MC sensory inputs
    else:
        gEinj_base=param_dict['M_gEinj_base']
        
    M.taumKS=param_dict['M_taumKs']
    M.tauI=param_dict['M_tauI']
    M.tauI_G=param_dict['M_tauI_G']
    M.gI_cst=param_dict['M_gI_cst']
    freq_modul = param_dict['freq_modul']
    phase_shifts = param_dict['MC_phase_dispersion']*randn(N_mitral) #3.0 # 0.3 // 2 // 5  (fig5C) #(fig3 and fig5 use 1.5) # fig4 uses 2.5 
    
    # M-M Weak inhibition
    dm,dd=param_dict['weakinh_mean_delay'],param_dict['weakinh_delay_dispersion']
    Cmm_weak_inh = Connection(M,M,'rI',
                                                  weight=param_dict['weakinh_weight'],
                                                  delay=(dm-dd,dm+dd), # mean delay +/- delay dispersion
                                                  sparseness=param_dict['weakinh_connectproba'])
    
    # Only for LFP recording
    Cmm_weak_inh2 = Connection(M,M,'rI2',weight=1.0)
    
    gEinj_max=M.gEinj.copy()
    @network_operation()
    def mitral_activation(clock):
        if clock.t>t_init_M:
            M.gEinj=(gEinj_max+gEinj_base)/2.+(gEinj_max-gEinj_base)/2.*cos(2*pi*freq_modul*Hz*clock.t+phase_shifts)
        else :
            M.gEinj=0.*siemens*meter**-2

    #  Granule params
    if use_granule:
        # Granule model and initialization
        N_granule=param_dict['N_granule']
        G_input=param_dict['G_input']
        phase_shifts_g=param_dict['GC_phase_dispersion']*randn(N_granule) # shift of phase to avoid artifact synchrony 0.2 base / 5 for vigile conditions
        G_I_base=param_dict['G_I_base']
        G_I_min=param_dict['G_I_min']
        G_I_max=param_dict['G_I_max']

        # Granule model and initialization
        G=NeuronGroup(N_granule,QIF_eqs,threshold=0.*mV,reset=-70.*mV,freeze=True,compile=True)
        G.V=-70.*mV
        G.Iinj=G_I_base

        # Dendro-dendritic connections
        if not(param_dict['use_AMPA_STP']):
            Cmg_AMPA = Connection(M,G,'sE',
                                                    weight=param_dict['exc_weight'],
                                                    sparseness=param_dict['stronginh_connectproba'],
                                                    delay=1.*ms)
        else: #withSTP
            Cmg_AMPA = Connection(M,G,'sE',
                                                    weight=2.5*param_dict['exc_weight'], # 2.5 factor to compensate for the depressed strength of AMPA in the beta regime
                                                    sparseness=param_dict['stronginh_connectproba'],
                                                    delay=1.*ms) #
            mystp=STP(Cmg_AMPA,taud=150*ms,tauf=1*ms,U=1.)
        
        Cgm_GABA = Connection(G,M,'sI_G',delay=1.*ms) # synaptic weight is given later
        connected=Cmg_AMPA.W.transpose().copy() # copy and transpose M->G connection array
        mask_connect=(connected>0)
        Cgm_GABA.connect(G,M,mask_connect*param_dict['stronginh_weight']) # finally create symmetrical connections
        
        # Time dependent activation of granule cells
        @network_operation()
        def granule_activation(clock):
            if clock.t<t_init_G:
                G.Iinj=G_I_base
            else:               
                if G_input=='sinusoid':
                    G.Iinj=(G_I_max+G_I_min)/2.+(G_I_max-G_I_min)/2.*cos(2*pi*freq_modul*Hz*clock.t+phase_shifts_g+param_dict['GC_phase_shift'])
                elif G_input=='ramp': # standard ramp : from -3.5nA to 0.5nA
                        G.Iinj=G_I_min* (1.-(clock.t-t_init_G)/(sim_time-t_init_G)) + G_I_max*(clock.t-t_init_G)/(sim_time-t_init_G) # current ramp, linear from first to second current value 
                else:
                    G.Iinj=G_I_base
                

    # Monitors
    LFP=PopulationStateMonitor(M,'LFP') # actually it is "fake" Isyn which omits random delays (only spikes convolved with an IPSC waveform, averaged across all cells)
    
    if return_details: # Additional detailed recordings
        recorded_mitrals=param_dict['recorded_mitrals']
        SpikesM = SpikeMonitor(M)
        Mpot=StateMonitor(M,'V',record=recorded_mitrals)
        avgIsyn=PopulationStateMonitor(M,'Isyn')
        Mvars=MultiStateMonitor(M,['Gsyn','Isyn','W','X','Y'],record=recorded_mitrals)
        if use_granule:
            recorded_granules=param_dict['recorded_granules']
            SpikesG = SpikeMonitor(G)
            Gpot=StateMonitor(G,'V',record=recorded_granules)
            Gvars=MultiStateMonitor(G,['sE'],record=recorded_granules)
            
    # Simulation
    run(sim_time,report=report)

    if return_details:
                
        result_dict={}
        result_dict['SpikesM_spikes']=SpikesM.spikes
        if use_granule:
            result_dict['SpikesG_spikes']=SpikesG.spikes
            result_dict['N_granule']=N_granule
            result_dict['Gpot_times']=Gpot.times
            result_dict['Gpot_values']=Gpot.values
        result_dict['N_mitral']=N_mitral
        result_dict['LFP_times']=LFP.times
        result_dict['LFP_values']=LFP.values
        result_dict['Mpot_times']=Mpot.times
        result_dict['Mpot_values']=Mpot.values
        result_dict['avgIsyn_times']=avgIsyn.times
        result_dict['avgIsyn_values']=avgIsyn.values
        result_dict['Mvars_times']=Mvars.times
        for name,mon in Mvars.items():
            result_dict['Mvars_'+name+'_values']=Mvars[name].values        
        result_dict['LFP_sr']=1./LFP.clock.dt
        result_dict['M_gEinj']=M.gEinj
        
        return result_dict
    else:
        print "One simulation done..."
        return LFP[0],1./LFP.clock.dt

if __name__=="__main__":
    
    
    from params_for_article_fig import *
    
    # Parameters to save simulation output
    save_output=False
    filename="test_dict_simplegamma" # _multi.out ou _single.out sont automatiquement ajouté après le nom suivant le type de simul     
    
    # Detection osc params
    osc_th=[0.1,0.1] # 0.2 for constant gamma or beta, 0.1 in presence of slow external modulations
    freq_cut=40.
    burn_time=1.*second
    
    if 1:  # To launch one network (with figures)

        # ## Choose params
        param_dict=param_dict_default.copy()
        # ## Or choose one of the predefind parameter set (defined for article figures)
        # param_dict=gamma_single_step_dict
        # param_dict=beta_single_ramp_dict
        # param_dict=competition_single_base_dict
        # param_dict=competitionSTP_single_lowintensity_dict
        # param_dict=competitionSTP_single_mediumintensity_dict
        # param_dict=competitionSTP_single_highintensity_dict
        # param_dict=competitionSTP_single_lowcentrifugal_dict
        # param_dict=competitionSTP_single_highcentrifugal_dict
        # param_dict=competitionSTP_single_awake
        # param_dict=competitionSTP_single_awake_lowcentrifugal
        
        # Run single network 
        result_dict = reseau_mitral_granule(param_dict)
        if save_output:
            print "Saving"
            import pickle
            fid=open(filename+"_single.out","wb")
            pickle.dump([param_dict,result_dict],fid)
            fid.close()
            print "Saved"
            
        # Plot detailed results
        single_LFP_analysis(result_dict,osc_th=osc_th,freq_cut=freq_cut,burn_time=burn_time)
    
    else: # To launch multiple network on different CPUs (only a synthesis figure is done)
        
        # Replace the parameter to vary with a numpy array of values to simulate
        param_dict=param_dict_default.copy()
        param_dict['weakinh_connectproba']=linspace(0.05,1.0,10)
        # ## Or choose one of the predefind parameter set (defined for article figures)
        # param_dict=gamma_multi_weakinhconnect_dict
        # param_dict=gamma_multi_numact_dict
        # param_dict=gamma_multi_weakinhweight_dict
        # param_dict=gamma_multi_tauGABA_dict
        # param_dict=gamma_multi_taumKs_dict
        # param_dict=gamma_multi_MgEinj_dict
        # param_dict=beta_multi_tauGABA_dict
        # param_dict=beta_multi_stronginhweight_dict
        # param_dict=beta_multi_GIinj_dict
        # param_dict=beta_multi_excweight_dict
        # param_dict=competition_multi_intensity
        # param_dict=competitionSTP_multi_intensity
        # param_dict=competitionSTP_multi_centrifugal_dict

        number_of_runs = 10 # Number of runs for each paramater (with different inputs and random connectivity)

        # Construct the list of param dictionaries for all simulations
        list_dicts=[]
        list_params=[]
        for key,val in param_dict.items():
            if (type(val)==ndarray)or((type(val)==list)and key[:8]!='recorded'):
                for param in val:
                    tmp_dict=param_dict.copy()
                    tmp_dict[key]=param
                    for i in range(number_of_runs):
                        list_dicts.append(tmp_dict)
                        list_params.append(param)
        
        # Third run multiple networks
        import multiprocessing as mp
        from functools import partial
        pool=mp.Pool(10)
        print "Start simulations (no idea of how long it will be)"
        results=pool.map(partial(reseau_mitral_granule,return_details=False,report=None),list_dicts)
        results=[(par,)+rr for par,rr in zip(list_params,results)] # Complete each result with its parameter
        print "finished : ",results

        if save_output:
            print "Saving"
            import pickle
            fid=open(filename+"_multi.out","wb")
            pickle.dump([param_dict,results],fid)
            fid.close()
            print "Saved"

        # Example of analysis from data in results
        distinct_figures=False
        tmp_results=[res[:3] for res in results]  # To keep only LFP and param
        multi_LFP_analysis(tmp_results,osc_th=osc_th,freq_cut=freq_cut,burn_time=burn_time,distinct_figures=distinct_figures)
        # Spectrum plot
        freq_min, freq_max =15.,90.
        LFP_frequency_content(tmp_results,burn_time=burn_time,freq_min=freq_min,freq_max=freq_max)
        
    show()