Add BBP catalogue (#1163)

Add a catalogue of mechanisms used by the BBP models.

arbor/include/arbor/mechcat.hpp

@@ -122,5 +122,6 @@ private:
 
 const mechanism_catalogue& global_default_catalogue();
 const mechanism_catalogue& global_allen_catalogue();
+const mechanism_catalogue& global_bbp_catalogue();
 
 } // namespace arb

mechanisms/CMakeLists.txt

@@ -9,6 +9,46 @@ add_custom_target(build_all_mods)
 
 set(mech_sources "")
 
+# BBP
+set(bbp_mechanisms CaDynamics_E2 Ca_HVA Ca_LVAst Ih Im K_Pst K_Tst Nap_Et2 NaTa_t NaTs2_t SK_E2 SKv3_1)
+set(bbp_mod_srcdir "${CMAKE_CURRENT_SOURCE_DIR}/bbp")
+set(mech_dir "${CMAKE_CURRENT_BINARY_DIR}/generated/bbp")
+file(MAKE_DIRECTORY "${mech_dir}")
+
+build_modules(
+        ${bbp_mechanisms}
+        SOURCE_DIR "${bbp_mod_srcdir}"
+        DEST_DIR "${mech_dir}"
+        ${external_modcc}
+        MODCC_FLAGS -t cpu -t gpu ${ARB_MODCC_FLAGS}
+        GENERATES .hpp _cpu.cpp _gpu.cpp _gpu.cu
+        TARGET build_catalogue_bbp_mods)
+
+set(bbp_catalogue_source ${CMAKE_CURRENT_BINARY_DIR}/bbp_catalogue.cpp)
+set(bbp_catalogue_options -A arbor -I ${mech_dir} -o ${bbp_catalogue_source} -B multicore -C bbp)
+if(ARB_WITH_GPU)
+    list(APPEND bbp_catalogue_options -B gpu)
+endif()
+
+add_custom_command(
+        OUTPUT ${bbp_catalogue_source}
+        COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/generate_catalogue ${bbp_catalogue_options} ${bbp_mechanisms}
+        DEPENDS generate_catalogue
+)
+
+add_custom_target(bbp_catalogue_cpp_target DEPENDS ${bbp_catalogue_source})
+add_dependencies(build_catalogue_bbp_mods bbp_catalogue_cpp_target)
+add_dependencies(build_all_mods build_catalogue_bbp_mods)
+
+list(APPEND mech_sources ${bbp_catalogue_source})
+foreach(mech ${bbp_mechanisms})
+    list(APPEND mech_sources ${mech_dir}/${mech}_cpu.cpp)
+    if(ARB_WITH_GPU)
+        list(APPEND mech_sources ${mech_dir}/${mech}_gpu.cpp)
+        list(APPEND mech_sources ${mech_dir}/${mech}_gpu.cu)
+    endif()
+endforeach()
+
 # ALLEN
 set(allen_mechanisms CaDynamics Ca_HVA Ca_LVA Ih Im Im_v2 K_P K_T Kd Kv2like Kv3_1 NaTa NaTs NaV Nap SK)
 set(allen_mod_srcdir "${CMAKE_CURRENT_SOURCE_DIR}/allen")

mechanisms/bbp/CaDynamics_E2.mod

+: Dynamics that track inside calcium concentration
+: modified from Destexhe et al. 1994
+
+NEURON {
+   SUFFIX CaDynamics_E2
+   USEION ca READ ica WRITE cai
+   RANGE decay, gamma, minCai, depth, initCai
+}
+
+UNITS {
+   (mV)    = (millivolt)
+   (mA)    = (milliamp)
+   (molar) = (1/liter)
+   (mM)    = (millimolar)
+   (um)    = (micron)
+}
+
+PARAMETER {
+   F      = 96485.3321233100184 (coulomb/mole) : Faraday's constant
+   gamma  = 0.05                               : percent of free calcium (not buffered)
+   decay  = 80                  (ms)           : rate of removal of calcium
+   depth  = 0.1                 (um)           : depth of shell
+   minCai = 1e-4                (mM)
+   initCai = 5e-5
+}
+
+STATE {
+   cai
+}
+
+INITIAL {
+   cai = initCai
+}
+
+BREAKPOINT {
+   SOLVE states METHOD cnexp
+}
+
+DERIVATIVE states {
+   cai' = -5000*ica*gamma/(F*depth) - (cai - minCai)/decay
+}

mechanisms/bbp/Ca_HVA.mod

+: Reference: Reuveni, Friedman, Amitai, and Gutnick, J.Neurosci. 1993
+
+NEURON {
+   SUFFIX Ca_HVA
+   USEION ca READ eca WRITE ica
+   RANGE gCa_HVAbar
+}
+
+UNITS {
+   (S)  = (siemens)
+   (mV) = (millivolt)
+   (mA) = (milliamp)
+}
+
+PARAMETER {
+   gCa_HVAbar = 0.00001 (S/cm2)
+}
+
+STATE {
+   m
+   h
+}
+
+BREAKPOINT {
+   SOLVE states METHOD cnexp
+   ica = gCa_HVAbar*m*m*h*(v - eca)
+}
+
+DERIVATIVE states {
+   LOCAL mAlpha, mBeta, hAlpha, hBeta
+
+   mAlpha = m_alpha(v)
+   mBeta  = m_beta(v)
+   hAlpha = h_alpha(v)
+   hBeta  = h_beta(v)
+
+   m' = mAlpha - m*(mAlpha + mBeta)
+   h' = hAlpha - h*(hAlpha + hBeta)
+}
+
+INITIAL {
+   LOCAL mAlpha, mBeta, hAlpha, hBeta
+
+   mAlpha = m_alpha(v)
+   mBeta  = m_beta(v)
+   hAlpha = h_alpha(v)
+   hBeta  = h_beta(v)
+
+   m = mAlpha/(mAlpha + mBeta)
+   h = hAlpha/(hAlpha + hBeta)
+}
+
+FUNCTION m_alpha(v) {
+   m_alpha = 0.055*3.8*exprelr(-(27 + v)/3.8)
+}
+FUNCTION m_beta(v) {
+   m_beta  = 0.94*exp(-(v + 75)/17)
+}
+FUNCTION h_alpha(v) {
+   h_alpha = 0.000457*exp(-(v + 13)/50)
+}
+FUNCTION h_beta(v) {
+   h_beta  = 0.0065/(exp(-(v + 15)/28) + 1)
+}

mechanisms/bbp/Ca_LVAst.mod

+:Comment   LVA ca channel. Note: mtau is an approximation from the plots
+:Reference Avery and Johnston 1996, tau from Randall 1997
+:Comment   shifted by -10 mv to correct for junction potential
+:Comment   corrected rates using q10 = 2.3, target temperature 34, orginal 21
+
+NEURON {
+   SUFFIX Ca_LVAst
+   USEION ca READ eca WRITE ica
+   RANGE gCa_LVAstbar
+}
+
+UNITS {
+   (S) = (siemens)
+   (mV) = (millivolt)
+   (mA) = (milliamp)
+}
+
+PARAMETER {
+   gCa_LVAstbar = 0.00001 (S/cm2)
+}
+
+STATE {
+   m
+   h
+}
+
+BREAKPOINT {
+   SOLVE states METHOD cnexp
+   ica = gCa_LVAstbar*m*m*h*(v - eca)
+}
+
+DERIVATIVE states {
+    LOCAL qt, mInf, mTau, hInf, hTau
+
+    qt = 2.3^((34 - 21)/10)
+
+    mInf = m_inf(v)
+    hInf = h_inf(v)
+    mTau =  5 + 20/(1 + exp((v + 35)/5))
+    hTau = 20 + 50/(1 + exp((v + 50)/7))
+
+    m' = (mInf - m)*qt/mTau
+    h' = (hInf - h)*qt/hTau
+}
+
+INITIAL {
+    m = m_inf(v)
+    h = h_inf(v)
+}
+
+FUNCTION m_inf(v) {
+    m_inf = 1/(1 + exp(-(v + 40)/6.0))
+}
+FUNCTION h_inf(v) {
+    h_inf = 1/(1 + exp( (v + 90)/6.4))
+}
+

mechanisms/bbp/Ih.mod

+: Reference:  Kole,Hallermann,and Stuart, J. Neurosci. 2006
+
+NEURON {
+    SUFFIX Ih
+    NONSPECIFIC_CURRENT ihcn
+    RANGE gIhbar
+}
+
+UNITS {
+    (S) = (siemens)
+    (mV) = (millivolt)
+    (mA) = (milliamp)
+}
+
+PARAMETER {
+    gIhbar =   0.00001 (S/cm2)
+    ehcn   = -45.0     (mV)
+}
+
+STATE {
+    m
+}
+
+BREAKPOINT {
+    SOLVE states METHOD cnexp
+    ihcn = gIhbar*m*(v - ehcn)
+}
+
+DERIVATIVE states {
+    LOCAL mAlpha, mBeta
+
+    mAlpha = m_alpha(v)
+    mBeta  = m_beta(v)
+
+    m' = mAlpha - m*(mAlpha + mBeta)
+}
+
+INITIAL {
+    LOCAL mAlpha, mBeta
+
+    mAlpha = m_alpha(v)
+    mBeta  = m_beta(v)
+
+    m = mAlpha/(mAlpha + mBeta)
+}
+
+FUNCTION m_alpha(v) {
+    m_alpha = 0.001*6.43*11.9*exprelr((v + 154.9)/11.9)
+}
+FUNCTION m_beta(v) {
+    m_beta  = 0.001*193*exp(v/33.1)
+}
+

mechanisms/bbp/Im.mod

+: Reference: Adams et al. 1982 - M-currents and other potassium currents in bullfrog sympathetic neurones
+: Comment:   corrected rates using q10 = 2.3, target temperature 34, orginal 21
+
+NEURON {
+    SUFFIX Im
+    USEION k READ ek WRITE ik
+    RANGE gImbar
+}
+
+UNITS {
+    (S) = (siemens)
+    (mV) = (millivolt)
+    (mA) = (milliamp)
+}
+
+PARAMETER {
+    gImbar = 0.00001 (S/cm2)
+}
+
+STATE {
+    m
+}
+
+BREAKPOINT {
+    SOLVE states METHOD cnexp
+    ik = gImbar*m*(v - ek)
+}
+
+DERIVATIVE states {
+    LOCAL qt, mAlpha, mBeta
+
+    qt     = 2.3^((34-21)/10)
+    mAlpha = m_alpha(v)
+    mBeta  = m_beta(v)
+
+    m'     = qt*(mAlpha - m*(mAlpha + mBeta))
+}
+
+INITIAL {
+    LOCAL mAlpha, mBeta
+
+    mAlpha = m_alpha(v)
+    mBeta  = m_beta(v)
+
+    m = mAlpha/(mAlpha + mBeta)
+}
+
+FUNCTION m_alpha(v) {
+    m_alpha = 3.3e-3*exp( 2.5*0.04*(v + 35))
+}
+FUNCTION m_beta(v) {
+    m_beta  = 3.3e-3*exp(-2.5*0.04*(v + 35))
+}
+

mechanisms/bbp/K_Pst.mod

+:Comment    The persistent component of the K current
+:Reference  Voltage-gated K+ channels in layer 5 neocortical pyramidal neurones from young rats:subtypes and gradients,Korngreen and Sakmann, J. Physiology, 2000
+:Comment    shifted -10 mv to correct for junction potential
+:Comment    corrected rates using q10 = 2.3, target temperature 34, orginal 21
+
+
+NEURON {
+    SUFFIX K_Pst
+    USEION k READ ek WRITE ik
+    RANGE gK_Pstbar
+}
+
+UNITS {
+    (S) = (siemens)
+    (mV) = (millivolt)
+    (mA) = (milliamp)
+}
+
+PARAMETER {
+    gK_Pstbar = 0.00001 (S/cm2)
+}
+
+STATE {
+    m
+    h
+}
+
+BREAKPOINT {
+    SOLVE states METHOD cnexp
+    ik = gK_Pstbar*m*m*h*(v-ek)
+}
+
+DERIVATIVE states {
+    LOCAL qt, mTau, mInf, hTau, hInf
+    qt = 2.3^((34-21)/10)
+
+    mInf = m_inf(v)
+    if (v < -60){
+        mTau = 1.25 + 175.03*exp( (v + 10)*0.026)
+    } else {
+        mTau = 1.25 +  13.00*exp(-(v + 10)*0.026)
+    }
+    hInf = h_inf(v)
+    hTau = 360 + (1010 + 24*(v + 65))*exp(-((v + 85)/48)^2)
+
+    m' = qt*(mInf - m)/mTau
+    h' = qt*(hInf - h)/hTau
+}
+
+INITIAL {
+    m = m_inf(v)
+    h = h_inf(v)
+}
+
+FUNCTION m_inf(v) {
+    m_inf = 1/(1 + exp(-(v + 11)/12))
+}
+FUNCTION h_inf(v) {
+    h_inf = 1/(1 + exp( (v + 64)/11))
+}

mechanisms/bbp/K_Tst.mod

+: Comment    The transient component of the K current
+: Reference  Voltage-gated K+ channels in layer 5 neocortical pyramidal neurones from young rats:subtypes and gradients,Korngreen and Sakmann, J. Physiology, 2000
+: Comment:   shifted -10 mv to correct for junction potential
+: Comment    corrected rates using q10 = 2.3, target temperature 34, orginal 21
+
+NEURON {
+    SUFFIX K_Tst
+    USEION k READ ek WRITE ik
+    RANGE gK_Tstbar
+}
+
+UNITS {
+    (S)  = (siemens)
+    (mV) = (millivolt)
+    (mA) = (milliamp)
+}
+
+PARAMETER {
+    gK_Tstbar = 0.00001 (S/cm2)
+}
+
+STATE {
+    m
+    h
+}
+
+BREAKPOINT {
+    SOLVE states METHOD cnexp
+    ik = gK_Tstbar*(m^4)*h*(v-ek)
+}
+
+DERIVATIVE states {
+    LOCAL qt, mInf, hInf, mTau, hTau
+
+    qt   = 2.3^((34 - 21)/10)
+    mInf = m_inf(v)
+    mTau = 0.34 + 0.92*exp(-((v + 81)/59)^2)
+    hInf = h_inf(v)
+    hTau = 8 + 49*exp(-((v + 83)/23)^2)
+
+    m' = qt*(mInf - m)/mTau
+    h' = qt*(hInf - h)/hTau
+}
+
+INITIAL {
+    m = m_inf(v)
+    h = h_inf(v)
+}
+
+FUNCTION m_inf(v) {
+   m_inf = 1/(1 + exp(-(v + 10)/19))
+}
+FUNCTION h_inf(v) {
+   h_inf = 1/(1 + exp(-(v + 76)/-10))
+}

mechanisms/bbp/NaTa_t.mod

+NEURON  {
+    SUFFIX NaTa_t
+    USEION na READ ena WRITE ina
+    RANGE gNaTa_tbar
+}
+
+UNITS {
+    (S)  = (siemens)
+    (mV) = (millivolt)
+    (mA) = (milliamp)
+}
+
+PARAMETER {
+    gNaTa_tbar = 0.00001 (S/cm2)
+}
+
+STATE {
+    m
+    h
+}
+
+BREAKPOINT {
+    SOLVE states METHOD cnexp
+    ina = gNaTa_tbar*m*m*m*h*(v-ena)
+}
+
+DERIVATIVE states {
+    LOCAL qt, mAlpha, mBeta, mRate, hAlpha, hBeta, hRate
+    qt = 2.3^((34-21)/10)
+
+    mAlpha = m_alpha(v)
+    mBeta  = m_beta(v)
+    mRate  = mAlpha + mBeta
+
+    hAlpha = h_alpha(v)
+    hBeta  = h_beta(v)
+    hRate  = hAlpha + hBeta
+
+    m' = qt*(mAlpha - m*mRate)
+    h' = qt*(hAlpha - h*hRate)
+}
+
+INITIAL {
+    LOCAL mAlpha, mBeta, hAlpha, hBeta
+
+    mAlpha = m_alpha(v)
+    mBeta  = m_beta(v)
+    m      = mAlpha/(mAlpha + mBeta)
+
+    hAlpha = h_alpha(v)
+    hBeta  = h_beta(v)
+    h      = hAlpha/(hAlpha + hBeta)
+}
+
+FUNCTION m_alpha(v) {
+    m_alpha = 0.182*6*exprelr(-(v + 38)/6)
+}
+
+FUNCTION m_beta(v) {
+    m_beta = 0.124*6*exprelr( (v + 38)/6)
+}
+
+FUNCTION h_alpha(v) {
+    h_alpha = 0.015*6*exprelr( (v + 66)/6)
+}
+
+FUNCTION h_beta(v) {
+    h_beta = 0.015*6*exprelr(-(v + 66)/6)
+}

mechanisms/bbp/NaTs2_t.mod

+: Reference Colbert and Pan 2002
+: Comment   took the NaTa and shifted both activation/inactivation by 6 mv
+
+NEURON {
+    SUFFIX NaTs2_t
+    USEION na READ ena WRITE ina
+    RANGE gNaTs2_tbar
+}
+
+UNITS {
+    (S) = (siemens)
+    (mV) = (millivolt)
+    (mA) = (milliamp)
+}
+
+PARAMETER {
+    gNaTs2_tbar = 0.00001 (S/cm2)
+}
+
+STATE {
+    m
+    h
+}
+
+BREAKPOINT {
+    SOLVE states METHOD cnexp
+    ina = gNaTs2_tbar*m*m*m*h*(v - ena)
+}
+
+DERIVATIVE states {
+    LOCAL qt, mAlpha, mBeta, mRate, hAlpha, hBeta, hRate
+
+    qt = 2.3^((34-21)/10)
+
+    mAlpha = m_alpha(v)
+    mBeta  = m_beta(v)
+    mRate  = mAlpha + mBeta
+
+    hAlpha = h_alpha(v)
+    hBeta  = h_beta(v)
+    hRate  = hAlpha + hBeta
+
+    m' = qt*(mAlpha - m*mRate)
+    h' = qt*(hAlpha - h*hRate)
+}
+
+INITIAL {
+    LOCAL mAlpha, mBeta, hAlpha, hBeta
+
+    mAlpha = m_alpha(v)
+    mBeta  = m_beta(v)
+    m      = mAlpha/(mAlpha + mBeta)
+
+    hAlpha = h_alpha(v)
+    hBeta  = h_beta(v)
+    h      = hAlpha/(hAlpha + hBeta)
+}
+
+FUNCTION m_alpha(v) {
+    m_alpha = 0.182*6*exprelr(-(v + 32)/6)
+}
+FUNCTION m_beta(v) {
+    m_beta  = 0.124*6*exprelr( (v + 32)/6)
+}
+FUNCTION h_alpha(v) {
+    h_alpha = 0.015*6*exprelr( (v + 60)/6)
+}
+FUNCTION h_beta(v) {
+    h_beta  = 0.015*6*exprelr(-(v + 60)/6)
+}
+

mechanisms/bbp/Nap_Et2.mod

+NEURON {
+    SUFFIX Nap_Et2
+    USEION na READ ena WRITE ina
+    RANGE gNap_Et2bar
+}
+
+UNITS {
+    (S)  = (siemens)
+    (mV) = (millivolt)
+    (mA) = (milliamp)
+}
+
+PARAMETER {
+    gNap_Et2bar = 0.00001 (S/cm2)
+}
+
+STATE {
+    m
+    h
+}
+
+BREAKPOINT {
+    SOLVE states METHOD cnexp
+    ina = gNap_Et2bar*m*m*m*h*(v - ena)
+}
+
+DERIVATIVE states {
+    LOCAL qt, mInf, mAlpha, mBeta, mRho, hInf, hAlpha, hBeta, hRho
+
+    qt = 2.3^((34 - 21)/10)
+
+    mInf = m_inf(v)
+    mAlpha = 0.182*6*exprelr(-(v + 38)/6)
+    mBeta  = 0.124*6*exprelr( (v + 38)/6)
+    mRho   = mAlpha + mBeta
+
+    hInf = h_inf(v)
+    hAlpha = 2.88e-6*4.63*exprelr( (v + 17.0)/4.63)
+    hBeta  = 6.94e-6*2.63*exprelr(-(v + 64.4)/2.63)
+    hRho   = hAlpha + hBeta
+
+    m' = qt*(mInf - m)*mRho/6.0     : equivalent to mTau = 6.0/mRho; m' = qt*(mInf - m)/mTau
+    h' = qt*(hInf - h)*hRho         : equivalent to hTau = 1.0/hRho; h' = qt*(hInf - h)/hTau
+}
+
+INITIAL {
+    m = m_inf(v)
+    h = h_inf(v)
+}
+
+FUNCTION m_inf(v) {
+    m_inf = 1.0/(1 + exp(-(v + 52.6)/4.6))
+}
+
+FUNCTION h_inf(v) {
+    h_inf = 1.0/(1 + exp( (v + 48.8)/10))
+}

mechanisms/bbp/SK_E2.mod

+: SK-type calcium-activated potassium current
+: Reference : Kohler et al. 1996
+
+NEURON {
+    SUFFIX SK_E2
+    USEION k READ ek WRITE ik
+    USEION ca READ cai
+    RANGE gSK_E2bar
+}
+
+UNITS {
+    (mV) = (millivolt)
+    (mA) = (milliamp)
+    (mM) = (milli/liter)
+}
+
+PARAMETER {
+    gSK_E2bar = 0.000001 (mho/cm2)
+    zTau      = 1        (ms)
+}
+
+STATE {
+    z FROM 0 TO 1
+}
+
+BREAKPOINT {
+    SOLVE states METHOD cnexp
+    ik = gSK_E2bar*z*(v - ek)
+}
+
+DERIVATIVE states {
+    z' = (z_inf(cai) - z)/zTau
+}
+
+INITIAL {
+    z = z_inf(cai)
+}
+
+FUNCTION z_inf(ca) {
+    if (ca < 1e-7) {
+        z_inf = 0
+    } else {
+        z_inf =  1/(1 + (0.00043/ca)^4.8)
+    }
+}

mechanisms/bbp/SKv3_1.mod

+:Reference  Characterization of a Shaw-related potassium channel family in rat brain, The EMBO Journal, vol.11, no.7,2473-2486 (1992)
+
+NEURON {
+    SUFFIX SKv3_1
+    USEION k READ ek WRITE ik
+    RANGE gSKv3_1bar
+}
+
+UNITS    {
+    (S)  = (siemens)
+    (mV) = (millivolt)
+    (mA) = (milliamp)
+}
+
+PARAMETER {
+    gSKv3_1bar = 0.00001 (S/cm2)
+}
+
+STATE {
+    m
+}
+
+BREAKPOINT {
+    SOLVE states METHOD cnexp
+    ik = gSKv3_1bar*m*(v - ek)
+}
+
+DERIVATIVE states {
+    LOCAL mInf, mRho
+
+    mInf = m_inf(v)
+    mRho = 0.25*(1 + exp((v + 46.56)/(-44.140)))
+
+    m' = (mInf - m)*mRho
+}
+
+INITIAL {
+    m = m_inf(v)
+}
+
+FUNCTION m_inf(v) {
+    m_inf = 1/(1 + exp((18.7 - v)/9.7))
+}

python/mechanism.cpp

@@ -136,6 +136,7 @@ void register_mechanisms(pybind11::module& m) {
 
     m.def("default_catalogue", [](){return arb::global_default_catalogue();});
     m.def("allen_catalogue", [](){return arb::global_allen_catalogue();});
+    m.def("bbp_catalogue", [](){return arb::global_bbp_catalogue();});
 
     // arb::mechanism_desc
     // For specifying a mechanism in the cable_cell interface.