single_file_custom_catalogue_arbor.py

import os.path
import numpy as np
import re
import subprocess
import tempfile
import arbor
from math import sqrt
import matplotlib.pyplot as plt

ARBOR_BUILD_CATALOGUE = 'arbor-build-catalogue'

mechanism_source = '''
NEURON {
    SUFFIX noise
    USEION na READ ena WRITE ina
    USEION k  READ ek  WRITE ik
    NONSPECIFIC_CURRENT il
    RANGE gnabar, gkbar, gl, el, q10
}

UNITS {
    (mV) = (millivolt)
    (S) = (siemens)
}

PARAMETER {
    gnabar =   0.12   (S/cm2)
    gkbar  =   0.036  (S/cm2)
    gl     =   0.0003 (S/cm2)
    el     = -54.3    (mV)
    celsius           (degC)
}

STATE { m h n }

ASSIGNED { q10 }

BREAKPOINT {
    SOLVE states METHOD cnexp
    LOCAL gk, gna, n2

    n2 = n*n
    gk = gkbar*n2*n2
    gna = gnabar*m*m*m*h
    ina = gna*(v - ena)
    ik  = gk*(v - ek)
    il  = gl*(v - el)
}

INITIAL {
    LOCAL alpha, beta

    q10 = 3^((celsius - 6.3)/10.0)

    : sodium activation system
    alpha = m_alpha(v)
    beta  = m_beta(v)
    m     = alpha/(alpha + beta)

    : sodium inactivation system
    alpha = h_alpha(v)
    beta  = h_beta(v)
    h     = alpha/(alpha + beta)

    : potassium activation system
    alpha = n_alpha(v)
    beta  = n_beta(v)

    n     = alpha/(alpha + beta)
}

DERIVATIVE states {
    LOCAL alpha, beta, sum

    : sodium activation system
    alpha = m_alpha(v)
    beta  = m_beta(v)
    sum   = alpha + beta
    m'    = (alpha - m*sum)*q10

    : sodium inactivation system
    alpha = h_alpha(v)
    beta  = h_beta(v)
    sum   = alpha + beta
    h'    = (alpha - h*sum)*q10

    : potassium activation system
    alpha = n_alpha(v)
    beta  = n_beta(v)
    sum   = alpha + beta
    n'    = (alpha - n*sum)*q10
}

FUNCTION vtrap(x,y) { vtrap   = y*exprelr(x/y) }

FUNCTION m_alpha(v) { m_alpha = 0.1*vtrap(-(v + 40.0), 10.0) }
FUNCTION h_alpha(v) { h_alpha = 0.07*exp(-(v + 65.0)/20.0) }
FUNCTION n_alpha(v) { n_alpha = 0.01*vtrap(-(v + 55.0), 10.0) }

FUNCTION m_beta(v)  { m_beta  = 4.0*exp(-(v + 65.0)/18.0) }
FUNCTION h_beta(v)  { h_beta  = 1.0/(exp(-(v + 35.0)/10.0) + 1.0) }
FUNCTION n_beta(v)  { n_beta  = 0.125*exp(-(v + 65.0)/80.0) }
'''

def build_noise_catalogue():
    with tempfile.TemporaryDirectory() as tmpdir:
        with open(os.path.join(tmpdir, 'noise.mod'), 'w') as f:
            print(mechanism_source, file=f)
        out = subprocess.getoutput(f'{ARBOR_BUILD_CATALOGUE} noise {tmpdir}')
        m = re.search('[^ ]*\.so', out)
        if m is None:
            print(out)
            exit(1)
        return arbor.load_catalogue(m.group(0))


def make_cable_cell(gid):
    tree = arbor.segment_tree()
    soma = tree.append(arbor.mnpos, arbor.mpoint(-12, 0, 0, 6), arbor.mpoint(0, 0, 0, 6), tag=1)
    labels = arbor.label_dict(dict(soma="(tag 1)", root= "(root)"))
    decor = arbor.decor()
    decor.paint('"soma"', arbor.density("noise"))
    return arbor.cable_cell(tree, labels, decor)

class Recipe(arbor.recipe):
    def __init__(self, ncells=100):
        arbor.recipe.__init__(self)
        self.ncells = ncells
        self.props = arbor.neuron_cable_properties()
        self.props.catalogue.extend(build_noise_catalogue(), '')
    def num_cells(self): return self.ncells
    def cell_description(self, gid): return make_cable_cell(gid)
    def cell_kind(self, gid): return arbor.cell_kind.cable
    def connections_on(self, gid): return []
    def event_generators(self, gid): return []
    def probes(self, gid): return [arbor.cable_probe_membrane_voltage('"root"')]
    def global_properties(self, kind): return self.props

# (11) Instantiate recipe
recipe = Recipe(ncells=100)
ctx = arbor.context('avail_threads')
sim = arbor.simulation(recipe, ctx)
handles = [sim.sample((gid, 0), arbor.regular_schedule(1)) for gid in range(recipe.num_cells())]
sim.run(100)
voltages = np.array([sim.samples(handles[gid])[0][0][:, 1] for gid in range(recipe.num_cells())])

plt.plot(voltages.T)
plt.show()