.. currentmodule:: brian

.. index::
   pair: example usage; NeuronGroup
   pair: example usage; run
   pair: example usage; StateMonitor

.. _example-frompapers_Rothman_Manis_2003:

Example: Rothman_Manis_2003 (frompapers)
========================================

Cochlear neuron model of Rothman & Manis
----------------------------------------
Rothman JS, Manis PB (2003) The roles potassium currents play in
regulating the electrical activity of ventral cochlear nucleus neurons.
J Neurophysiol 89:3097-113.

All model types differ only by the maximal conductances.

Adapted from their Neuron implementation by Romain Brette

::

    from brian import *
    
    #defaultclock.dt=0.025*ms # for better precision
    
    '''
    Simulation parameters: choose current amplitude and neuron type
    (from type1c, type1t, type12, type 21, type2, type2o)
    '''
    neuron_type = 'type1c'
    Ipulse = 250 * pA
    
    C = 12 * pF
    Eh = -43 * mV
    EK = -70 * mV # -77*mV in mod file
    El = -65 * mV
    ENa = 50 * mV
    nf = 0.85 # proportion of n vs p kinetics
    zss = 0.5 # steady state inactivation of glt
    celsius = 22. # temperature
    q10 = 3. ** ((celsius - 22) / 10.)
    # hcno current (octopus cell)
    frac = 0.0
    qt = 4.5 ** ((celsius - 33.) / 10.)
    
    # Maximal conductances of different cell types in nS
    maximal_conductances = dict(
    type1c=(1000, 150, 0, 0, 0.5, 0, 2),
    type1t=(1000, 80, 0, 65, 0.5, 0, 2),
    type12=(1000, 150, 20, 0, 2, 0, 2),
    type21=(1000, 150, 35, 0, 3.5, 0, 2),
    type2=(1000, 150, 200, 0, 20, 0, 2),
    type2o=(1000, 150, 600, 0, 0, 40, 2) # octopus cell
    )
    gnabar, gkhtbar, gkltbar, gkabar, ghbar, gbarno, gl = [x * nS for x in maximal_conductances[neuron_type]]
    
    # Classical Na channel
    eqs_na = """
    ina = gnabar*m**3*h*(ENa-v) : amp
    dm/dt=q10*(minf-m)/mtau : 1
    dh/dt=q10*(hinf-h)/htau : 1
    minf = 1./(1+exp(-(vu + 38.) / 7.)) : 1
    hinf = 1./(1+exp((vu + 65.) / 6.)) : 1
    mtau =  ((10. / (5*exp((vu+60.) / 18.) + 36.*exp(-(vu+60.) / 25.))) + 0.04)*ms : ms
    htau =  ((100. / (7*exp((vu+60.) / 11.) + 10.*exp(-(vu+60.) / 25.))) + 0.6)*ms : ms
    """
    
    # KHT channel (delayed-rectifier K+)
    eqs_kht = """
    ikht = gkhtbar*(nf*n**2 + (1-nf)*p)*(EK-v) : amp
    dn/dt=q10*(ninf-n)/ntau : 1
    dp/dt=q10*(pinf-p)/ptau : 1
    ninf =   (1 + exp(-(vu + 15) / 5.))**-0.5 : 1
    pinf =  1. / (1 + exp(-(vu + 23) / 6.)) : 1
    ntau =  ((100. / (11*exp((vu+60) / 24.) + 21*exp(-(vu+60) / 23.))) + 0.7)*ms : ms
    ptau = ((100. / (4*exp((vu+60) / 32.) + 5*exp(-(vu+60) / 22.))) + 5)*ms : ms
    """
    
    # Ih channel (subthreshold adaptive, non-inactivating)
    eqs_ih = """
    ih = ghbar*r*(Eh-v) : amp
    dr/dt=q10*(rinf-r)/rtau : 1
    rinf = 1. / (1+exp((vu + 76.) / 7.)) : 1
    rtau = ((100000. / (237.*exp((vu+60.) / 12.) + 17.*exp(-(vu+60.) / 14.))) + 25.)*ms : ms
    """
    
    # KLT channel (low threshold K+)
    eqs_klt = """
    iklt = gkltbar*w**4*z*(EK-v) : amp
    dw/dt=q10*(winf-w)/wtau : 1
    dz/dt=q10*(zinf-z)/wtau : 1
    winf = (1. / (1 + exp(-(vu + 48.) / 6.)))**0.25 : 1
    zinf = zss + ((1.-zss) / (1 + exp((vu + 71.) / 10.))) : 1
    wtau = ((100. / (6.*exp((vu+60.) / 6.) + 16.*exp(-(vu+60.) / 45.))) + 1.5)*ms : ms
    ztau = ((1000. / (exp((vu+60.) / 20.) + exp(-(vu+60.) / 8.))) + 50)*ms : ms
    """
    
    # Ka channel (transient K+)
    eqs_ka = """
    ika = gkabar*a**4*b*c*(EK-v): amp
    da/dt=q10*(ainf-a)/atau : 1
    db/dt=q10*(binf-b)/btau : 1
    dc/dt=q10*(cinf-c)/ctau : 1
    ainf = (1. / (1 + exp(-(vu + 31) / 6.)))**0.25 : 1
    binf = 1. / (1 + exp((vu + 66) / 7.))**0.5 : 1
    cinf = 1. / (1 + exp((vu + 66) / 7.))**0.5 : 1
    atau =  ((100. / (7*exp((vu+60) / 14.) + 29*exp(-(vu+60) / 24.))) + 0.1)*ms : ms
    btau =  ((1000. / (14*exp((vu+60) / 27.) + 29*exp(-(vu+60) / 24.))) + 1)*ms : ms
    ctau = ((90. / (1 + exp((-66-vu) / 17.))) + 10)*ms : ms
    """
    
    # Leak
    eqs_leak = """
    ileak = gl*(El-v) : amp
    """
    
    # h current for octopus cells
    eqs_hcno = """
    ihcno = gbarno*(h1*frac + h2*(1-frac))*(Eh-v) : amp
    dh1/dt=(hinfno-h1)/tau1 : 1
    dh2/dt=(hinfno-h2)/tau2 : 1
    hinfno = 1./(1+exp((vu+66.)/7.)) : 1
    tau1 = bet1/(qt*0.008*(1+alp1))*ms : ms
    tau2 = bet2/(qt*0.0029*(1+alp2))*ms : ms
    alp1 = exp(1e-3*3*(vu+50)*9.648e4/(8.315*(273.16+celsius))) : 1
    bet1 = exp(1e-3*3*0.3*(vu+50)*9.648e4/(8.315*(273.16+celsius))) : 1 
    alp2 = exp(1e-3*3*(vu+84)*9.648e4/(8.315*(273.16+celsius))) : 1
    bet2 = exp(1e-3*3*0.6*(vu+84)*9.648e4/(8.315*(273.16+celsius))) : 1
    """
    
    eqs = """
    dv/dt=(ileak+ina+ikht+iklt+ika+ih+ihcno+I)/C : volt
    vu = v/mV : 1 # unitless v
    I : amp
    """
    eqs += eqs_leak + eqs_ka + eqs_na + eqs_ih + eqs_klt + eqs_kht + eqs_hcno
    
    neuron = NeuronGroup(1, eqs, implicit=True)
    neuron.v = El
    
    run(50 * ms) # Go to rest
    
    M = StateMonitor(neuron, 'v', record=0)
    neuron.I = Ipulse
    
    run(100 * ms, report='text')
    
    plot(M.times / ms, M[0] / mV)
    show()