ModelDB is moving. Check out our new site at The corresponding page is

MEC PV-positive fast-spiking interneuron network generates theta-nested fast oscillations

 Download zip file 
Help downloading and running models
We use a computational model of a network of Fast-Spiking Parvalbumin-positive Basket Cells to study its synchronizing properties. The intrinsic properties of neurons, properties of chemical synapses and of gap junctions are calibrated using electrophysiological recordings in mice Medial Entorhinal Cortex slices. The neurons synchronize, generating Fast Oscillations nested in an external theta drive. We show how gap junctions are necessary for the generation of the oscillations, how hyperpolarizing chemical synapses give rise to more robust fast oscillations, compared to shunting ones, and how short-term depression in the chemical synapses confine the fast oscillation on a narrow range of phases from the external theta drive.
Model Information (Click on a link to find other models with that property)
Model Type: Connectionist Network;
Brain Region(s)/Organism: Entorhinal cortex; Mouse;
Cell Type(s): Entorhinal cortex fast-spiking interneuron;
Channel(s): I Na,t; I K;
Gap Junctions: Gap junctions;
Receptor(s): GabaA;
Transmitter(s): Gaba;
Simulation Environment: Brian 2;
Model Concept(s): Brain Rhythms; Excitability; Gamma oscillations; Theta oscillations; Short-term Synaptic Plasticity;
Implementer(s): Via, Guillem; Baravalle, Roman;
Search NeuronDB for information about:  GabaA; I Na,t; I K; Gaba;
# Code adapted by Guillem Via from the original code implementing the Brunel (2000) network model from the exercise repository for
# the book by Wulfram Gerstner et al. (2014)

# Single model neuron equations are extracted in parts by Erisir et al, (1999), Golomb et al. (2007) and Economo and White (2012).
# Parameters are calibrated to reproduce passive properties and f/I curve heterogeneity from a dataset of 11 parvalbumin positive
# fast-spiking interneurons recorded in mouse dorsal Medial Entorhinal Cortex slices.

# Original description of the source code is as follows:

# This file is part of the exercise code repository accompanying
# the book: Neuronal Dynamics (see
# located at

# This free software: you can redistribute it and/or modify it under
# the terms of the GNU General Public License 2.0 as published by the
# Free Software Foundation. You should have received a copy of the
# GNU General Public License along with the repository. If not,
# see

# Should you reuse and publish the code for your own purposes,
# please cite the book or point to the webpage

# Wulfram Gerstner, Werner M. Kistler, Richard Naud, and Liam Paninski.
# Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition.
# Cambridge University Press, 2014.

import brian2 as b2
from brian2 import NeuronGroup, Synapses, PoissonInput, PoissonGroup, SpikeGeneratorGroup, Equations, TimedArray
from brian2.monitors import StateMonitor, SpikeMonitor, PopulationRateMonitor
from brian2 import Hz, nS, nF, Mohm, pA, mV, ms
import matplotlib.pyplot as plt
from time import time
import numpy as np
from numpy import random, exp, sqrt, log, pi
from random import shuffle, sample
import os
TigerFish = True
if TigerFish: 	# For Tigerfish use only
	cache_dir = os.path.expanduser('~/.cython/brian-pid-{}'.format(os.getpid()))
	b2.prefs.codegen.runtime.cython.cache_dir = cache_dir
	b2.prefs.codegen.runtime.cython.multiprocess_safe = False

def gewnet(
    eps=1.e-9, # Small parameter to avoid divide by zero in gating variable activation function
    connectivity_matrix=0, # To load connectivity matrices for chemical and electrical synapses from different files
    sim_time=1500., # Total simulated time (in ms)
    N=100, # Number of neurons
    Nrec=10, # Number of neurons to record state variables from (only if return_state==True)
    v_thresh=-30., # Threshold membrane potential whose crossing counts a spike and triggers a postsynaptic IPSC
    epsilon=.36, # Connection probability for chemical synapses. Used only if connectivity matrix is not read from file
    I=0., # External input current in the constant current protocol, i.e. if const_current==True (in pA)
    gext=7., # In case external drive is applied as a constant conductance, multiplied by a driving force with reversal Ese (in nS)
    Ips=0., # Value or array of values for the Ip applied to each neuron. If array, then need to be numpy.array() type.
    gsin=7., # In case external drive is applied as a conductance (times a driving force). Conductance oscillates between 0 and 2*gsin.
    gms=1.65, # Value for peak conductances if not heterogeneous, 1.65 by defaut, i.e. mean of the log-normal used for most simulations (in nS)
    gmmean=0., # if we assume Esyn!=-75, gms need to be rescaled by a factor a, i.e. lognormal dist mu->mu+log(a)
    gmstd=1., # If lognormal gm, its std
    gm_scale=1., # If scale-gm==True, scale peak conductance values gms by this factor gm_scale
    tau_fall=2.0, # Synaptic decay time constant. For the bi-exponential IPSC conductance
    tau_rise=.3, # Synaptic rise time constant. For the bi-exponential IPSC conductance
    Es=-75., # Reversal potential for the chemical synapses
    Ese=0.,  # Synaptic reversal for excitation. In case external through conductance.

    st=time(); # Start timer to compute simulation running time

    b2.defaultclock.dt = dt * ms # Set brian simulator default clock


    # Set fixed and homogeneous neuron parameters

    ENa=+50.*mV; EK=-90.*mV;

    km2=.1; kh1=.012; kh2=.2; kn2=FactorScaleKV3*.001; ka2=.02;  kn1 = FactorScaleKV3*1.;

    sigm1=4.*mV; sigh2=3.5*mV; sign1=12.*mV; sigm2=-13.*mV; sigh1=-20.*mV; sign2=-8.5*mV; siga1=12.*mV; siga2=-80.*mV;

    # Tuned parameters with heterogeneity and properties as observed in electrophysiological recordings

    # Load neuron parameters that can be heterogeneous (need to set folder and file names)
    params_act = np.transpose(np.loadtxt( "params/intrinsic/active/params_actRoman.dat" )); # Heterogeneous active parameters: Voltage-gated current peak conductances and thetas
    #params_act = np.loadtxt( "params/intrinsic/active/params_act.dat" ); # Heterogeneous active parameters: Voltage-gated current peak conductances and thetas
    if mod_gL==True:
        folder = "mod_gL/sigma04_higg12/gLmin%d/" % int(np.round(10.*gLmin));
        gLs = np.loadtxt( "params/intrinsic/passive/gLs/" + folder + "gLs%d.dat" % connectivity_matrix )[:N]; # Leakage conductances
        ELs = np.loadtxt( "params/intrinsic/passive/ELs/" + folder + "ELs%d.dat" % connectivity_matrix )[:N]; # Leakage reversal potentials
        folder = "original/";
        gLs = np.loadtxt( "params/intrinsic/passive/gLs/" + folder + "gLsRoman.dat" )[:N]; # Leakage conductances
        ELs = np.loadtxt( "params/intrinsic/passive/ELs/" + folder + "ELsRoman.dat" )[:N]; # Leakage reversal potentials
    taums = factorTau*np.loadtxt( "params/intrinsic/passive/taums/taumsRomanUniform.dat" )[:N]; # Membrane time constants

    Rins = 1.0e+3/gLs;
    Caps = taums/Rins;

    gNas, gKv3s, gKv1s, thm1s, thh2s, thn1s, tha1s = params_act[:,:N];

    #gLs = np.asarray([14.705184211981187])
    #ELs = np.asarray([-71.97508650437244])
    #taums = np.asarray([5.222681687884308310])#np.asarray([5.22])
    #Rins = 1.0e+3/gLs;
    #Caps = taums/Rins;
    #gNas, gKv3s, gKv1s, thm1s, thh2s, thn1s, tha1s = np.asarray([1.68049678e+04]), np.asarray([6.31700393e+02]), np.asarray([5.90431313e+01]),  np.asarray([-5.29951244e+01]), np.asarray([-5.57111801e+01]), np.asarray([5.90077045e+00]), np.asarray([5.13553887e+01])


    # If considering a homogeneous network, take passive parameters from a single model neuron
    if homo_neurons==True: gLs, ELs, taums, Rins, Caps = gLs[kneu]*np.ones(N), ELs[kneu]*np.ones(N), taums[kneu]*np.ones(N), Rins[kneu]*np.ones(N), Caps[kneu]*np.ones(N);

    # If considering a network with homogeneous intrinsic neuronal active parameters (from kneu set in the function call)
    if homo_neurons==True or homo_act==True: gNas, gKv3s, gKv1s, thm1s, thh2s, thn1s, tha1s = gNas[kneu]*np.ones(N), gKv3s[kneu]*np.ones(N), gKv1s[kneu]*np.ones(N), thm1s[kneu]*np.ones(N), thh2s[kneu]*np.ones(N), thn1s[kneu]*np.ones(N), tha1s[kneu]*np.ones(N);
    #print(gLs[kneu], ELs[kneu], taums[kneu], Rins[kneu], Caps[kneu],gNas[kneu], gKv3s[kneu], gKv1s[kneu], thm1s[kneu], thh2s[kneu], thn1s[kneu], tha1s[kneu]); 

    # Import synaptic parameters

    # Chemical synapses
    if act_syns==True:

        if read_syns==True: syns = np.loadtxt( "params/chem_syns/conn_mat/p036/syns%d.dat" % connectivity_matrix, dtype=type(1))[:N];

        if read_delays==True: delays = np.loadtxt( "params/chem_syns/delays/p036/delays%d.dat" % connectivity_matrix)[:len(syns[0])];

        if read_gms==True: gms = np.loadtxt( "params/chem_syns/gms/p036/" + gm_distribution + "/gms%d.dat" % connectivity_matrix)[:len(syns[0])]; # Read peak conductances

    # Gap junctions
    if gjs==True:
        if read_gjs==True: synsgj = np.loadtxt( "params/gjs/conn_mat/" + folder + "synsgj%d.dat" % connectivity_matrix).astype(int);
        if read_ggaps==True: ggaps = np.loadtxt( "params/gjs/ggaps/" + folder + "ggaps%d.dat" % connectivity_matrix)[:len(synsgj[0])];

    # Set units

    I = I*pA; Ips = Ips*pA; Isin = Isin*pA; # External input (depending on protocol)

    gext = gext*nS; gsin = gsin*nS; Ese = Ese*mV; # External conductance (depending on protocol - emulates Chr optogenetic input - reversal of 0 mV)

    gNas = gNas*nS; gKv3s = gKv3s*nS; gKv1s = gKv1s*nS; # Peak conductances for voltage-gated currents

    thm1s = thm1s*mV; thh2s = thh2s*mV; thn1s = thn1s*mV; tha1s = tha1s*mV; eps = eps*mV;

    taums = taums*ms; gLs = gLs*nS; Caps = Caps*nF; Rins = Rins*Mohm; ELs = ELs*mV; # Passive properties
    v_thresh = v_thresh*mV;
    #print("gLs,ELs,taums,gNas, gKv3s, gKv1s, thm1s, thh2s, thn1s, tha1s Caps")
    #print(gLs[kneu],ELs[kneu],taums[kneu],gNas[kneu], gKv3s[kneu], gKv1s[kneu], thm1s[kneu], thh2s[kneu], thn1s[kneu], tha1s[kneu],Caps[kneu])
    # Synaptic parameters
    if act_syns==True:
        tau_fall = tau_fall*ms; tau_rise = tau_rise*ms; Es = Es*mV; gms = gms*nS;
        if read_delays==True: delays = delays*ms;

    ggaps = ggaps*nS; # Gap junction peak conductances

    # Theta drive frequency
    fsin = fsin*Hz;

    # Short-term-depression recorvery time constant
    taudm = taudm*ms;

    # Normalizing factor for the IPSC conductance height
    c_fall = 1./tau_fall; c_rise = 1./tau_rise;
    norm_syn = 1./( exp(-c_fall*log(c_rise/c_fall)/(c_rise-c_fall)) - exp(-c_rise*log(c_rise/c_fall)/(c_rise-c_fall)) );


    # Applied step current (activated at simulation middle time)
    if step_current==True: intervals = [0.,1.]; stimulus = TimedArray(intervals,dt=sim_time/2.*;


    # Membrane-potential dynamics equations

    gat_vars =    """
                  alpham = -((v-thm1-eps)/sigm1)/(exp(-(v-thm1-eps)/sigm1)-1.)/ms : Hz
                  betam = km2*exp(v/sigm2)/ms : Hz
                  dm/dt = alpham*(1.-m) - betam*m : 1
                  alphah = kh1*exp(v/sigh1)/ms : Hz
                  betah = -kh2*(v-thh2)/(exp(-(v-thh2)/sigh2)-1.)/(ms*mV) : Hz
                  dh/dt = alphah*(1.-h) - betah*h : 1
                  alphan = -kn1*(v-thn1)/(exp(-(v-thn1)/sign1)-1.)/(ms*mV) : Hz
                  betan = kn2*exp(v/sign2)/ms : Hz
                  dn/dt = alphan*(1.-n) - betan*n : 1
                  alphaa = -(v-tha1)/(exp(-(v-tha1)/siga1)-1.)/(ms*mV) : Hz
                  betaa = ka2*exp(v/siga2)/ms : Hz
                  da/dt = alphaa*(1.-a) - betaa*a : 1
                  INa = - gNa*m*m*m*h*(v-ENa) : ampere
                  IKv3 = - gKv3*n*n*n*n*(v-EK) : ampere
                  IKv1 = - gKv1*a*a*a*a*(v-EK) : ampere
                  IL = - gL*(v-EL) : ampere

    # Heterogeneous parameters

    param_eqs =   """
                  Ip : ampere
                  EL : volt
                  taum : second
                  gL : siemens
                  Cap : farad
                  gNa : siemens
                  gKv3 : siemens
                  gKv1 : siemens
                  thm1 : volt
                  thh2 : volt
                  thn1 : volt
                  tha1 : volt

    # Gap-junction currents

    if gjs == True:
        param_eqs += """
                  Igap : ampere

    # Chemical synapse currents

    if act_syns == True:
        syn_eqs = """
                     Isyn = - (v-Es)*norm_syn*(g_fall-g_rise) : ampere
                     dg_fall/dt = -g_fall/tau_fall : siemens
                     dg_rise/dt = -g_rise/tau_rise : siemens
    else: syn_eqs = "";

    # External currents (conditional on the protocol: constant, step or theta-modulated; conductance or current-based)

    ext_eqs = "Iext = "

    if const_cond==True: ext_eqs += " - gext*(v-Ese) "
    if const_current==True: ext_eqs += " + I "
    if step_current==True: ext_eqs += " + Ip*stimulus(t) "
    if sin_cond==True: ext_eqs += " - gsin*(1.+sin(2.*pi*fsin*t-pi/2.))*(v-Ese) "
    if sin_current==True: ext_eqs += " + Isin*(1.+sin(2.*pi*fsin*t-pi/2.)) "	

    ext_eqs +=   " : ampere \n "

    # Master equation adding all above currents

    currs = " INa + IKv1 + IKv3 + IL + Iext "
    if act_syns == True: currs += " + Isyn "
    if gjs == True: currs += " + Igap "

    volt_eqs = "\n dv/dt = ( " + currs + " )/Cap : volt \n"

    equations = gat_vars + ext_eqs + syn_eqs + param_eqs + volt_eqs;

    print("Neuron equations are:\n %s\n\n" % ( equations ) );

    # Define neuron population. Use a threshold v_thresho to count a spike and trigger postsynaptic responses.
    # Since the spike has non-zero width a refractory period prevents the code from counting more than one spike while v>v_thresb
    net = NeuronGroup(N, equations, threshold="v>v_thresh", refractory=refractory, method=method)

    net.Ip = Ips;

    net.EL, net.gL, net.taum, net.Cap = ELs, gLs, taums, Caps;
    net.gNa, net.gKv3, net.gKv1, net.thm1, net.thh2, net.thn1, net.tha1 = gNas, gKv3s, gKv1s, thm1s, thh2s, thn1s, tha1s;

    # Initialize state variable at random values
    if randics==True:
        net.v = ELs+3.*mV*np.random.randn(len(ELs)); # Normally distributed membrane potential
        net.m = .01+.1*np.random.rand(len(ELs)); # Uniformly distributed gating variables
        net.h = .99-.1*np.random.rand(len(ELs));
        net.n = .01+.1*np.random.rand(len(ELs));
        net.a = .01+.1*np.random.rand(len(ELs));
        net.v = -30.*mV;
        net.m = .465391;
        net.h = .184172;
        net.n = .582671;
        net.a = .399927;
    if one_perturbed==True:
        kpert = 0;
        print("The perturbed neuron is kneu=%d.\n\n" % kpert);
        net.v[kpert] = -58.592274*mV;
        net.m[kpert] = .470296;
        net.h[kpert] = .035708;
        net.n[kpert] = .910829;
        net.a[kpert] = .452683;


    # If chemical synapses are activated
    if act_syns==True:

        # Definition of synapses
        print("Setting synapses...")

        # Include STD at the chem syns
        if std == True:

            synapses_eqs = '''
            # Fraction of synaptic neurotransmitter resources available:
            dxs/dt = (1 - xs)/tau_d : 1 (event-driven)
            gmax : siemens
            U_SE : 1
            tau_d : second
            synapses_action = '''
            rs = U_SE * xs
            xs -= rs

        # If STD is not included
        else: synapses_eqs = "gmax : siemens"; synapses_action = "g_fall_post+=gmax;g_rise_post+=gmax";

        # Print the used synaptic equations
        print("Synapse equations are:\n %s\n\n" % synapses_eqs )
        print("Synapse actions are:\n %s\n\n" % synapses_action )

        # Define the equation population
        S = Synapses(net,net,model=synapses_eqs,on_pre=synapses_action,method=method)

        # Connectivity matrix can be read from file. Otherwise Erdos-Renyi (connectivity="ER") or fixed number of presynaptic partners, i.e. in-degree (connectivity="FID")
        if read_syns==True: S.connect(i=syns[0], j=syns[1])
        elif connectivity=="ER": S.connect(p=epsilon, condition="i!=j")
        elif connectivity=="FID":
            C = int(np.round(epsilon*N));
            for k in range(N): pres = [ kp for kp in range(N) if kp!=k ]; shuffle(pres); S.connect(i=pres[:C],j=k);
        else: print("ERROR: Unknown connectivity type"); return "ERROR: Unknown connectivity type";

        if read_syns==True and read_delays==True: S.delay = delays;
        else: S.delay = (dmin+(dmax-dmin)*np.random.rand(len(S)))*ms;
        if std:
            S.xs = 1.; # Initialize auxiliary variable for STD
            S.U_SE = USEm; S.tau_d = taudm;
        if dist_syns==True:
            if read_syns==False or read_gms==False: gms = np.random.lognormal(gmmean,gmstd,size=len(S))*nS; # Set them to lognormal if not read
            if std==True:
                if scale_gm==False: S.gmax = gms/S.U_SE; # In order to see the effects of upscaling gms preserving their distribution 
                else: S.gmax = gm_scale*gms/S.U_SE; # If STD is added, gm is divided by U_SE. Then results with short tau_d and without STD agree.
                if scale_gm==False: S.gmax = gms;
                else: S.gmax = gm_scale*gms;
            if std: S.gmax = gms/S.U_SE;
            else: S.gmax = gms;


    if gjs==True:

        gap_eqs = """
                 g_gap : siemens # gap junction conductance
                 Igap_post = g_gap * (v_pre - v_post) : ampere (summed)

        print("Gap-junction equations are:\n %s\n\n" % gap_eqs );

        Sgj = Synapses(net, net, gap_eqs, method=method);

        if read_gjs==True: Sgj.connect(i=synsgj[0], j=synsgj[1])
        else: Sgj.connect(p=pgj,condition="i!=j"); # WRONG! THESE ARE NOT BIDIRECTIONAL, SYMMETRIC AND NON-RECTIFYING

        Sgj.g_gap = ggaps; # They can be read from file or set on function call. Can be an array of length synsgj or a single float value

    print("Setting monitors...")
    rMon = PopulationRateMonitor(net); # Population rate monitor to return population-averaged (and possibly time-smoothed) firing rates
    spMon = SpikeMonitor(net, record=True) # Monitor for spike times

    # Monitor for state variables. If return_state==True. If gap junctions are not activated, return membrane potential, synaptic current and external current.
    # If they are activated return also net gap junction current onto the neuron.
    # Record and return state variables for the first Nrec neurons
    if return_state==True:
        state_vars = ["v","Iext"];#["v","m","n","h","a","Iext"]#
        if gjs==True: state_vars.append("Igap");
        if act_syns==True: state_vars.append("Isyn");
        stMon = StateMonitor(net, state_vars, record=range(Nrec), dt=dt_rec*

    # RUN
    print("Running the simulation...")*

    rates = [rMon.t/,rMon.smooth_rate(window=windowType,width=window_width)/b2.Hz] # Save and return time-smoothed population-averaged firing rates
    # Return also time- and population-averaged spike-count rate. Averaged across the whole simulation.
    nu = np.sum(rates[1][-40000:])/len(rates[1][-40000:]); print( "\nnu = %.2lf Hz" % nu );

    # END
    print("Done (took %.2lf mins)" % ((time()-st)/60.) )
    if return_state==True: return rates, spMon, stMon
    elif return_syns==True: return rates, spMon, S
    else: return rates, spMon