.. currentmodule:: brian2

.. COBAHH:

Example: COBAHH
===============


        .. only:: html

            .. |launchbinder| image:: file:///usr/share/doc/python-brian-doc/docs/badge.svg
            .. _launchbinder: https://mybinder.org/v2/gh/brian-team/brian2-binder/master?filepath=examples/COBAHH.ipynb

            .. note::
               You can launch an interactive, editable version of this
               example without installing any local files
               using the Binder service (although note that at some times this
               may be slow or fail to open): |launchbinder|_

        

This is an implementation of a benchmark described
in the following review paper:

Simulation of networks of spiking neurons: A review of tools and strategies (2006).
Brette, Rudolph, Carnevale, Hines, Beeman, Bower, Diesmann, Goodman, Harris, Zirpe,
Natschläger, Pecevski, Ermentrout, Djurfeldt, Lansner, Rochel, Vibert, Alvarez, Muller,
Davison, El Boustani and Destexhe.
Journal of Computational Neuroscience

Benchmark 3: random network of HH neurons with exponential synaptic conductances

Clock-driven implementation
(no spike time interpolation)

R. Brette - Dec 2007

::

    
    from brian2 import *
    
    # Parameters
    area = 20000*umetre**2
    Cm = (1*ufarad*cm**-2) * area
    gl = (5e-5*siemens*cm**-2) * area
    
    El = -60*mV
    EK = -90*mV
    ENa = 50*mV
    g_na = (100*msiemens*cm**-2) * area
    g_kd = (30*msiemens*cm**-2) * area
    VT = -63*mV
    # Time constants
    taue = 5*ms
    taui = 10*ms
    # Reversal potentials
    Ee = 0*mV
    Ei = -80*mV
    we = 6*nS  # excitatory synaptic weight
    wi = 67*nS  # inhibitory synaptic weight
    
    # The model
    eqs = Equations('''
    dv/dt = (gl*(El-v)+ge*(Ee-v)+gi*(Ei-v)-
             g_na*(m*m*m)*h*(v-ENa)-
             g_kd*(n*n*n*n)*(v-EK))/Cm : volt
    dm/dt = alpha_m*(1-m)-beta_m*m : 1
    dn/dt = alpha_n*(1-n)-beta_n*n : 1
    dh/dt = alpha_h*(1-h)-beta_h*h : 1
    dge/dt = -ge*(1./taue) : siemens
    dgi/dt = -gi*(1./taui) : siemens
    alpha_m = 0.32*(mV**-1)*4*mV/exprel((13*mV-v+VT)/(4*mV))/ms : Hz
    beta_m = 0.28*(mV**-1)*5*mV/exprel((v-VT-40*mV)/(5*mV))/ms : Hz
    alpha_h = 0.128*exp((17*mV-v+VT)/(18*mV))/ms : Hz
    beta_h = 4./(1+exp((40*mV-v+VT)/(5*mV)))/ms : Hz
    alpha_n = 0.032*(mV**-1)*5*mV/exprel((15*mV-v+VT)/(5*mV))/ms : Hz
    beta_n = .5*exp((10*mV-v+VT)/(40*mV))/ms : Hz
    ''')
    
    P = NeuronGroup(4000, model=eqs, threshold='v>-20*mV', refractory=3*ms,
                    method='exponential_euler')
    Pe = P[:3200]
    Pi = P[3200:]
    Ce = Synapses(Pe, P, on_pre='ge+=we')
    Ci = Synapses(Pi, P, on_pre='gi+=wi')
    Ce.connect(p=0.02)
    Ci.connect(p=0.02)
    
    # Initialization
    P.v = 'El + (randn() * 5 - 5)*mV'
    P.ge = '(randn() * 1.5 + 4) * 10.*nS'
    P.gi = '(randn() * 12 + 20) * 10.*nS'
    
    # Record a few traces
    trace = StateMonitor(P, 'v', record=[1, 10, 100])
    run(1 * second, report='text')
    plot(trace.t/ms, trace[1].v/mV)
    plot(trace.t/ms, trace[10].v/mV)
    plot(trace.t/ms, trace[100].v/mV)
    xlabel('t (ms)')
    ylabel('v (mV)')
    show()
    

