From 083aaf936514567cf22c626e3258eaacf6a05d9e Mon Sep 17 00:00:00 2001 From: Andrew Davison Date: Wed, 8 Jul 2026 15:00:20 +0200 Subject: [PATCH 1/6] Add point-neuron support to the Arbor backend (phase 1: IF_curr_delta) Implement IF_curr_delta by mapping it onto Arbor's native leaky integrate-and-fire cell (arbor.lif_cell, cell_kind.lif). Note that native lif cells cannot take an injected current, so non-zero `i_offset` and `DCSource` injection raise NotImplementedError. Verified on Arbor 0.10.0 and 0.12.2. --- pyNN/arbor/populations.py | 29 ++++++++++++++++- pyNN/arbor/projections.py | 32 ++++++++++++++++-- pyNN/arbor/recording.py | 6 +++- pyNN/arbor/standardmodels.py | 41 ++++++++++++++++++++++++ test/system/scenarios/test_cell_types.py | 25 ++++++++++++++- test/unittests/test_arbor.py | 29 +++++++++++++++++ 6 files changed, 156 insertions(+), 6 deletions(-) diff --git a/pyNN/arbor/populations.py b/pyNN/arbor/populations.py index 02237271..5b762c0b 100644 --- a/pyNN/arbor/populations.py +++ b/pyNN/arbor/populations.py @@ -5,6 +5,7 @@ from warnings import warn import numpy as np import arbor +from arbor import units as U from .. import common, errors from ..standardmodels import StandardCellType @@ -100,6 +101,25 @@ def arbor_cell_description(self, gid): schedule_params[key] = value * unit schedule = self.celltype.arbor_schedule(**schedule_params) return arbor.spike_source_cell("spike-source", schedule) + elif self.celltype.arbor_cell_kind == arbor.cell_kind.lif: + cell_descr = self._arbor_cell_description + if not cell_descr._evaluated: + cell_descr.evaluate() + params = list(cell_descr)[index] + if params.get("i_offset", 0.0) != 0.0: + raise NotImplementedError( + "Arbor's native lif_cell (IF_curr_delta) cannot inject a " + "constant current, so i_offset must be 0. Use the cable-cell " + "IF models for current injection.") + # The source label ("detector") matches what projections.py uses for + # injectable presynaptic cells; "syn" is the single delta-synapse target. + cell = arbor.lif_cell("detector", "syn") + for name, unit in self.celltype.lif_param_units.items(): + setattr(cell, name, float(params[name]) * unit) + initial_v = getattr(self, "_lif_initial_v", None) + if initial_v is not None: + cell.V_m = float(initial_v._partially_evaluate(index, simplify=True)) * U.mV + return cell else: args = self._arbor_cell_description[index] return _compat.make_cable_cell( @@ -117,7 +137,8 @@ def _create_cells(self): parameter_space.shape = (self.size,) - if self.celltype.arbor_cell_kind == arbor.cell_kind.spike_source: + if self.celltype.arbor_cell_kind in (arbor.cell_kind.spike_source, + arbor.cell_kind.lif): self._arbor_cell_description = parameter_space else: self._arbor_cell_description = parameter_space["cell_description"] @@ -147,6 +168,12 @@ def _create_cells(self): self._mask_local = np.ones_like(id_range, dtype=bool) def _set_initial_value_array(self, variable, initial_values): + if self.celltype.arbor_cell_kind == arbor.cell_kind.lif: + # Native lif cells have no decor; the initial v is applied as V_m when + # the cell is built in arbor_cell_description. + if variable == "v": + self._lif_initial_v = initial_values + return if variable != "v": warn("todo: handle initial values for ion channel states") # may have to handle this at the same time as setting parameters diff --git a/pyNN/arbor/projections.py b/pyNN/arbor/projections.py index 9a3e995d..c74d9ccf 100644 --- a/pyNN/arbor/projections.py +++ b/pyNN/arbor/projections.py @@ -67,6 +67,19 @@ def _convergent_connect(self, presynaptic_indices, postsynaptic_index, ConnectionGroup(pre_idx, postsynaptic_index, self.receptor_type, location_selector, **other_attributes) ) + def _lif_post_cm_pF(self, gid): + """The C_m (in pF) of the postsynaptic native lif cell with the given gid. + + The native C_m is stored in nF (its PyNN unit); the delta-weight charge + relation ΔV[mV] = weight / C_m[pF] needs it in pF. + """ + post_pop = self.post.parent if hasattr(self.post, "parent") else self.post + native = post_pop._arbor_cell_description + if not native._evaluated: + native.evaluate() + cm_nF = float(list(native)[post_pop.id_to_index(gid)]["C_m"]) + return 1000.0 * cm_nF + def arbor_connections(self, gid): """Return a list of incoming connections to the cell with the given gid""" try: @@ -80,16 +93,29 @@ def arbor_connections(self, gid): source = "spike-source" connections = [] - all_labels = list(self.post._arbor_cell_description[postsynaptic_index]["labels"]) + is_lif_post = self.post.celltype.arbor_cell_kind == arbor.cell_kind.lif + if is_lif_post: + # A lif_cell's delta synapse adds weight/C_m [mV] to V_m, i.e. the + # weight is a charge. PyNN's IF_curr_delta weight is a voltage step + # (mV), so scale by the post cell's C_m [pF] to recover it. + weight_scale = self._lif_post_cm_pF(gid) + else: + weight_scale = 1.0 + all_labels = list(self.post._arbor_cell_description[postsynaptic_index]["labels"]) for cg in self.connections[postsynaptic_index]: if cg.location_selector in (None, "all"): - target_labels = [lbl for lbl in all_labels if lbl.startswith(cg.receptor_type)] + if is_lif_post: + # A native lif_cell has a single built-in delta synapse; + # excitatory vs inhibitory is set by the sign of the weight. + target_labels = ["syn"] + else: + target_labels = [lbl for lbl in all_labels if lbl.startswith(cg.receptor_type)] for target in target_labels: connections.append( arbor.connection( (self.pre[cg.presynaptic_index], source), target, - cg.weight, + cg.weight * weight_scale, cg.delay * U.ms ) ) diff --git a/pyNN/arbor/recording.py b/pyNN/arbor/recording.py index 92f31cb2..ac1565f5 100644 --- a/pyNN/arbor/recording.py +++ b/pyNN/arbor/recording.py @@ -90,7 +90,11 @@ def _get_arbor_probes(self, gid): tag = str(probe_index) probe_index += 1 if variable.name == "v": - probe = arbor.cable_probe_membrane_voltage(locset, tag) + if self.population.celltype.arbor_cell_kind == arbor.cell_kind.lif: + # Native lif cells have no morphology/locset. + probe = arbor.lif_probe_voltage(tag) + else: + probe = arbor.cable_probe_membrane_voltage(locset, tag) else: mech_name, state_name = variable.name.split(".") arbor_model = mech_name # to do: find_arbor_model(mech_name) diff --git a/pyNN/arbor/standardmodels.py b/pyNN/arbor/standardmodels.py index 94c22e34..5c31c9f2 100644 --- a/pyNN/arbor/standardmodels.py +++ b/pyNN/arbor/standardmodels.py @@ -48,6 +48,35 @@ class SpikeSourceArray(cells.SpikeSourceArray): arbor_schedule_units = {"times": U.ms} +class IF_curr_delta(cells.IF_curr_delta): + __doc__ = cells.IF_curr_delta.__doc__ + + # Maps onto Arbor's native leaky integrate-and-fire cell (arbor.lif_cell, + # cell_kind.lif). Its synapses are delta, but an incoming event adds + # weight/C_m to V_m (the event weight is a charge, not a voltage), so + # IF_curr_delta's mV voltage-step weight is recovered by scaling the + # connection weight by C_m (see Projection._lif_post_cm_pF). + translations = build_translations( + ('v_rest', 'E_L'), + ('v_reset', 'E_R'), + ('v_thresh', 'V_th'), + ('tau_refrac', 't_ref'), + ('tau_m', 'tau_m'), + ('cm', 'C_m'), + # A native lif_cell has no way to inject a constant current, so i_offset + # is carried through untranslated and rejected at cell-build time if + # non-zero (see Population.arbor_cell_description). + ('i_offset', 'i_offset'), + ) + arbor_cell_kind = arbor.cell_kind.lif + # Units for the native lif_cell attributes (Arbor requires unit-typed values, + # and handles the conversion to its internal units itself). + lif_param_units = { + 'E_L': U.mV, 'E_R': U.mV, 'V_th': U.mV, + 't_ref': U.ms, 'tau_m': U.ms, 'C_m': U.nF, + } + + class BaseCurrentSource(object): pass @@ -62,6 +91,18 @@ class DCSource(BaseCurrentSource, electrodes.DCSource): ) def inject_into(self, cells, location=None): # rename to `locations` ? + # Native lif cells (IF_curr_delta) have no decor and cannot take an + # i_clamp, so current injection is impossible for them. + if hasattr(cells, "parent"): + target_pop = cells.parent + elif hasattr(cells, "_arbor_cell_description"): + target_pop = cells + else: + target_pop = cells[0].parent + if target_pop.celltype.arbor_cell_kind == arbor.cell_kind.lif: + raise NotImplementedError( + "Current injection into Arbor's native lif_cell (IF_curr_delta) " + "is not supported; use the cable-cell IF models instead.") if hasattr(cells, "parent"): cell_descr = cells.parent._arbor_cell_description.base_value index = cells.parent.id_to_index(cells.all_cells.astype(int)) diff --git a/test/system/scenarios/test_cell_types.py b/test/system/scenarios/test_cell_types.py index 96e3cb8e..b635fa97 100644 --- a/test/system/scenarios/test_cell_types.py +++ b/test/system/scenarios/test_cell_types.py @@ -312,7 +312,7 @@ def test_SpikeSourceArray_delivers_spike_times(sim): sim.end() -@run_with_simulators("nest", "brian2") +@run_with_simulators("arbor", "nest", "brian2") def test_IF_curr_delta_voltage_step(sim): """A single delta-synapse input steps V by the weight (mV). @@ -338,6 +338,29 @@ def test_IF_curr_delta_voltage_step(sim): sim.end() +@run_with_simulators("arbor", "brian2") +def test_IF_curr_delta_fires_and_resets(sim): + """A supra-threshold delta train makes the wired cell fire; the other is silent.""" + sim.setup(timestep=0.1) + source = sim.Population(1, sim.SpikeSourceArray( + spike_times=[10.0, 11.0, 12.0, 13.0])) + cells = sim.Population(2, sim.IF_curr_delta( + v_rest=-65.0, v_reset=-65.0, v_thresh=-50.0, + tau_m=20.0, cm=1.0, tau_refrac=5.0), + initial_values={'v': -65.0}) + sim.Projection(source, cells[0:1], sim.AllToAllConnector(), + sim.StaticSynapse(weight=8.0, delay=1.0), + receptor_type="excitatory") + cells.record('spikes') + sim.run(60.0) + spiketrains = cells.get_data().segments[0].spiketrains + assert len(spiketrains) == 2 + counts = sorted(len(st) for st in spiketrains) + assert counts[0] == 0 # unconnected cell stays silent + assert counts[1] >= 1 # wired cell fires + sim.end() + + @run_with_simulators("nest", "neuron", "brian2") def test_composed_neuron_model_homogeneous_receptors(sim, plot_figure=False): sim.setup() diff --git a/test/unittests/test_arbor.py b/test/unittests/test_arbor.py index 5bc81e05..453b1ec7 100644 --- a/test/unittests/test_arbor.py +++ b/test/unittests/test_arbor.py @@ -93,6 +93,15 @@ def test_electrode_param_units(self): self.assertIs(units["duration"], U.ms) self.assertIs(units["current"], U.nA) + def test_if_curr_delta_lif_param_units(self): + units = arbor_standardmodels.IF_curr_delta.lif_param_units + self.assertIs(units["E_L"], U.mV) + self.assertIs(units["E_R"], U.mV) + self.assertIs(units["V_th"], U.mV) + self.assertIs(units["t_ref"], U.ms) + self.assertIs(units["tau_m"], U.ms) + self.assertIs(units["C_m"], U.nF) + @unittest.skipUnless(have_arbor, "Requires Arbor") class TestTranslations(unittest.TestCase): @@ -125,6 +134,26 @@ def test_multicompartment_schema(self): for key in ("morphology", "cm", "Ra", "ionic_species"): self.assertIn(key, schema) + def test_if_curr_delta_is_native_lif(self): + self.assertIs(arbor_standardmodels.IF_curr_delta.arbor_cell_kind, + arbor.cell_kind.lif) + + def test_if_curr_delta_translation(self): + m = arbor_standardmodels.IF_curr_delta( + v_rest=-65.0, cm=1.0, tau_m=20.0, tau_refrac=2.0, + v_reset=-70.0, v_thresh=-50.0, i_offset=0.0) + native = m.native_parameters + native.shape = (1,) + native.evaluate(simplify=True) + d = native.as_dict() + self.assertAlmostEqual(d["E_L"], -65.0) + self.assertAlmostEqual(d["E_R"], -70.0) + self.assertAlmostEqual(d["V_th"], -50.0) + self.assertAlmostEqual(d["t_ref"], 2.0) + self.assertAlmostEqual(d["tau_m"], 20.0) + self.assertAlmostEqual(d["C_m"], 1.0) # cm passes through in nF + self.assertAlmostEqual(d["i_offset"], 0.0) + @unittest.skipUnless(have_arbor, "Requires Arbor") class TestMorphologyLocsets(unittest.TestCase): From 9a4cceb961e13ffb1af52324f82ac32f99b24ce4 Mon Sep 17 00:00:00 2001 From: Andrew Davison Date: Thu, 9 Jul 2026 11:15:29 +0200 Subject: [PATCH 2/6] Add cable-cell point neurons to the Arbor backend (phase 2) Implement point neurons (IF_cond_exp, IF_curr_exp and the composable PointNeuron) as single-compartment Arbor cable cells, complementing the native lif_cell used for IF_curr_delta in phase 1. - pas leak with specific cm/g_leak derived from the whole-cell cm/tau_m so the absolute C_m and g_leak match PyNN's values; - new lif.mod reset mechanism: the threshold_detector emits the network spike and, via POST_EVENT, reloads a refractory countdown that clamps V to v_reset (Arbor NMODL has neither absolute time nor WATCH); - new expsyn_curr.mod current-based synapse; conductance synapses reuse expsyn; i_offset and DCSource become iclamp stimuli. IF_cond_exp/IF_curr_exp are flat standard models whose translate() calls the base class translate() and wraps the result; the composable PointNeuron assembles the same builder input from its components. Cross-checked against NEURON (exact) and Brian2/NEST on f-I and EPSP amplitude; verified on Arbor 0.10.0 and 0.12.2. Also enable the arbor backend on the scenario tests it now passes (simulation control, DCSource injection, SpikeSourcePoisson parameters, procedural+OO recording). --- pyNN/arbor/cells.py | 187 ++++++++++++++++-- pyNN/arbor/nmodl/expsyn_curr.mod | 40 ++++ pyNN/arbor/nmodl/lif.mod | 58 ++++++ pyNN/arbor/populations.py | 15 +- pyNN/arbor/recording.py | 4 +- pyNN/arbor/standardmodels.py | 140 ++++++++++++- .../scenarios/test__simulation_control.py | 10 +- test/system/scenarios/test_cell_types.py | 72 +++++++ test/system/scenarios/test_electrodes.py | 2 +- .../scenarios/test_parameter_handling.py | 2 +- test/system/scenarios/test_recording.py | 2 +- test/system/scenarios/test_scenario2.py | 2 +- test/unittests/test_arbor.py | 88 +++++++++ 13 files changed, 590 insertions(+), 32 deletions(-) create mode 100644 pyNN/arbor/nmodl/expsyn_curr.mod create mode 100644 pyNN/arbor/nmodl/lif.mod diff --git a/pyNN/arbor/cells.py b/pyNN/arbor/cells.py index 6cbbfb33..13504990 100644 --- a/pyNN/arbor/cells.py +++ b/pyNN/arbor/cells.py @@ -1,5 +1,6 @@ from collections import defaultdict +import numpy as np from lazyarray import larray import arbor from arbor import units as U @@ -33,9 +34,39 @@ def region_name_to_tag(name): return map.get(name, -1) -class CellDescriptionBuilder: +class BaseCellDescriptionBuilder: + """Shared protocol for the Arbor cable-cell builders consumed by + :mod:`populations.py`. + + A builder is a callable ``builder(i)`` returning the ``tree``/``decor``/ + ``labels``/``discretization`` needed to construct the ``i``-th cell, plus + ``set_shape``, ``set_initial_values`` and ``add_current_source``. This base + holds only the parts common to :class:`CellDescriptionBuilder` (multicompartment) + and :class:`PointCellDescriptionBuilder` (point neuron); ``set_shape``, + ``_build_tree``, ``_build_decor`` and ``__call__`` are builder-specific. + """ + + def __init__(self): + self.initial_values = {} + self.current_sources = defaultdict(list) + + def set_initial_values(self, variable, initial_values): + assert isinstance(initial_values, larray) + self.initial_values[variable] = initial_values + + def add_current_source(self, model_name, location_generator, index, parameters): + for i in index: + self.current_sources[i].append({ + "model_name": model_name, + "parameters": parameters, + "location_generator": location_generator + }) + + +class CellDescriptionBuilder(BaseCellDescriptionBuilder): def __init__(self, parameters, ion_channels, post_synaptic_entities=None): + super().__init__() assert isinstance(parameters, ParameterSpace) self.parameters = parameters self.ion_channels = ion_channels @@ -50,8 +81,6 @@ def __init__(self, parameters, ion_channels, post_synaptic_entities=None): "root": "(root)", "mid-dend": "(location 0 0.5)" } - self.initial_values = {} - self.current_sources = defaultdict(list) def _build_tree(self, i): self.parameters["morphology"].dtype = Morphology @@ -152,18 +181,6 @@ def _build_decor(self, i): return decor - def set_initial_values(self, variable, initial_values): - assert isinstance(initial_values, larray) - self.initial_values[variable] = initial_values - - def add_current_source(self, model_name, location_generator, index, parameters): - for i in index: - self.current_sources[i].append({ - "model_name": model_name, - "parameters": parameters, - "location_generator": location_generator - }) - def set_shape(self, value): self.parameters.shape = value for ion_channel in self.ion_channels.values(): @@ -183,6 +200,146 @@ def __call__(self, i): } +# Geometry of the synthetic single compartment used for point neurons. The +# absolute surface area cancels out of the membrane dynamics (it is divided back +# out when deriving the specific cm/leak below), so it doesn't affect results; we +# match the pyNN.neuron backend's SingleCompartmentNeuron (L = 100 um, +# diam = 1000/pi um) so the two backends use an identical synthetic cell. Both +# NEURON and Arbor take the cylinder area to be the lateral surface only (no end +# caps), giving area = pi*L*diam = exactly 1e-3 cm2 -- a round value chosen so the +# derived specific capacitance equals the whole-cell cm numerically (1 nF -> 1 +# uF/cm2), avoiding an irrational specific cm. +POINT_CELL_LENGTH_UM = 100.0 +POINT_CELL_DIAMETER_UM = 1000.0 / np.pi +_POINT_CELL_AREA_CM2 = np.pi * POINT_CELL_DIAMETER_UM * POINT_CELL_LENGTH_UM * 1e-8 +# A clamp conductance large enough that the post-spike reset settles within one +# timestep (tau_clamp = C_m / LIF_RESET_CONDUCTANCE << dt for any realistic C_m). +LIF_RESET_CONDUCTANCE = 1000.0 # uS + + +class PointCellDescriptionBuilder(BaseCellDescriptionBuilder): + """Builds an Arbor cable cell that behaves as a PyNN point neuron. + + A point neuron (LIF, and its IF_cond_exp/IF_curr_exp specialisations) is + realised as a single-compartment cable cell: + + * the sub-threshold leak is a ``pas`` density mechanism, with the specific + capacitance and leak conductance derived from the whole-cell ``cm`` [nF] and + ``tau_m`` [ms] so that the *absolute* C_m and g_leak match PyNN's values; + * the network spike is emitted by the cell's ``threshold_detector`` (at + ``v_thresh``); + * the post-spike reset and refractory clamp are provided by the ``lif`` + point mechanism (see nmodl/lif.mod), driven by the detector via POST_EVENT; + * one synapse per receptor (``expsyn`` for conductance-based, + ``expsyn_curr`` for current-based) is placed at the soma; + * ``i_offset`` and injected current sources become ``iclamp`` stimuli. + + Its ``__call__(i)`` / ``set_shape`` / ``set_initial_values`` / + ``add_current_source`` interface mirrors :class:`CellDescriptionBuilder`, so a + point-neuron Population reuses the backend's existing cable-cell machinery. + + ``neuron_parameters`` is a native :class:`ParameterSpace` with the LIF keys + (E_L, E_R, V_th, t_ref, tau_m, C_m, i_offset). ``post_synaptic_receptors`` maps + each receptor label to an ``(arbor_synapse_model, native_synapse_parameters)`` + pair. This is the form produced both by the composable + :class:`~pyNN.arbor.standardmodels.PointNeuron` (from its components) and by the + flat IF_cond_exp/IF_curr_exp standard models (from their translations). + """ + + def __init__(self, neuron_parameters, post_synaptic_receptors): + super().__init__() + self.neuron_parameters = neuron_parameters + self.post_synaptic_receptors = post_synaptic_receptors + self.shape = None + + def set_shape(self, value): + self.shape = value + self.neuron_parameters.shape = value + for (_model, synapse_parameters) in self.post_synaptic_receptors.values(): + synapse_parameters.shape = value + + def _specific_properties(self, cm_nF, tau_m_ms): + """Specific membrane capacitance [uF/cm2] and leak conductance [S/cm2] + reproducing the whole-cell C_m [nF] and g_leak = C_m/tau_m [uS].""" + c_spec = cm_nF * 1e-3 / _POINT_CELL_AREA_CM2 + g_spec = (cm_nF / tau_m_ms) * 1e-6 / _POINT_CELL_AREA_CM2 + return c_spec, g_spec + + def _build_tree(self, i): + d = POINT_CELL_DIAMETER_UM + tree = arbor.segment_tree() + tree.append(arbor.mnpos, arbor.mpoint(0, 0, 0, d / 2), + arbor.mpoint(POINT_CELL_LENGTH_UM, 0, 0, d / 2), tag=1) + return tree + + def _build_decor(self, i): + p = self.neuron_parameters + + E_L = p["E_L"][i] + c_spec, g_spec = self._specific_properties(p["C_m"][i], p["tau_m"][i]) + labels = { + "all": "(all)", "soma": "(tag 1)", "root": "(root)", + } + decor = arbor.decor() + decor.set_property( + Vm=self.initial_values["v"][i] * U.mV, + cm=c_spec * U.uF / U.cm2, + ) + # sub-threshold leak (e is a GLOBAL parameter of pas, set via the mechanism name) + decor.paint("(all)", arbor.density(f"pas/e={E_L}", {"g": g_spec})) + # integrate-and-fire reset + refractory clamp + decor.place( + "(root)", + arbor.synapse("lif", {"v_reset": p["E_R"][i], "t_ref": p["t_ref"][i], + "g_reset": LIF_RESET_CONDUCTANCE}), + "lif_reset", + ) + # network spike source + decor.place("(root)", arbor.threshold_detector(p["V_th"][i] * U.mV), "detector") + + # one synapse per receptor, labelled by receptor name so that + # Projection.arbor_connections can match it by receptor_type prefix + for label, (model, synapse_parameters) in self.post_synaptic_receptors.items(): + params = {key: synapse_parameters[key][i] + for key in synapse_parameters.keys() + if key != "locations"} # a point neuron has a single location + decor.place("(root)", arbor.synapse(model, params), label) + labels[label] = "(root)" + + # constant offset current + i_offset = p["i_offset"][i] if "i_offset" in p.keys() else 0.0 + if i_offset != 0.0: + self._place_iclamp(decor, "(root)", + tstart=0.0, duration=1e12, current=i_offset, + label="i_offset") + # injected current sources (e.g. DCSource) + for source in self.current_sources[i]: + location_generator = source["location_generator"] + params = source["parameters"].evaluate(simplify=True).as_dict() + for locset, label in location_generator.generate_locations(None, label=source["model_name"]): + self._place_iclamp( + decor, locset, + tstart=params["tstart"], duration=params["duration"], + current=params["current"], label=label) + + return decor, arbor.label_dict(labels) + + def _place_iclamp(self, decor, locset, tstart, duration, current, label): + mechanism = _compat.get_electrode_mechanism("iclamp") + mech = mechanism(tstart * U.ms, duration * U.ms, current * U.nA) + _compat.place_current_source(decor, locset, mech, label) + + def __call__(self, i): + tree = self._build_tree(i) + decor, labels = self._build_decor(i) + return { + "tree": tree, + "decor": decor, + "labels": labels, + "discretization": arbor.cv_policy_single(), + } + + class NativeCellType(BaseCellType): arbor_cell_kind = arbor.cell_kind.cable units = {"v": "mV"} diff --git a/pyNN/arbor/nmodl/expsyn_curr.mod b/pyNN/arbor/nmodl/expsyn_curr.mod new file mode 100644 index 00000000..4c58a134 --- /dev/null +++ b/pyNN/arbor/nmodl/expsyn_curr.mod @@ -0,0 +1,40 @@ +: Current-based synapse with exponential decay, for IF_curr_exp / current-based +: point neurons. An incoming event steps the synaptic current by `weight` (nA), +: which then decays exponentially with time constant `tau`. A positive weight +: injects depolarising (inward) current, matching PyNN's isyn convention, so the +: contributed membrane current is i = -isyn. +NEURON { + POINT_PROCESS expsyn_curr + RANGE tau + NONSPECIFIC_CURRENT i +} + +UNITS { + (nA) = (nanoamp) + (ms) = (millisecond) +} + +PARAMETER { + tau = 5.0 (ms) +} + +STATE { + isyn (nA) +} + +INITIAL { + isyn = 0 +} + +BREAKPOINT { + SOLVE state METHOD cnexp + i = -isyn +} + +DERIVATIVE state { + isyn' = -isyn/tau +} + +NET_RECEIVE(weight) { + isyn = isyn + weight +} diff --git a/pyNN/arbor/nmodl/lif.mod b/pyNN/arbor/nmodl/lif.mod new file mode 100644 index 00000000..be07c4c5 --- /dev/null +++ b/pyNN/arbor/nmodl/lif.mod @@ -0,0 +1,58 @@ +: Integrate-and-fire reset mechanism, for building point neurons as Arbor cable +: cells (a single-compartment cell whose sub-threshold leak is a separate `pas` +: mechanism and whose network spike is emitted by the cell's threshold_detector). +: +: This mechanism performs only the post-spike reset and refractory clamp: when the +: cell spikes, Arbor delivers a POST_EVENT here, which reloads a refractory +: countdown `refrac` (in ms). While the countdown is positive the mechanism injects +: a strong current g_reset*(v - v_reset) that holds the membrane at v_reset; when it +: is not, it contributes nothing. Arbor NMODL exposes neither absolute time nor +: WATCH (NEURON's adexp.mod approach), hence the countdown-plus-POST_EVENT design. +: +: The clamp is gated by a voltage-independent multiplier `clamp` (0/1) rather than +: writing the piecewise current directly, so that the automatically-derived +: conductance di/dv is exactly clamp*g_reset (0 when not refractory); guarding the +: current expression with a plain `if (refrac > 0)` instead makes modcc leak a +: spurious conductance even when the current is zero. +NEURON { + POINT_PROCESS lif + RANGE v_reset, t_ref, g_reset + NONSPECIFIC_CURRENT i +} + +UNITS { + (mV) = (millivolt) + (nA) = (nanoamp) + (uS) = (microsiemens) + (ms) = (millisecond) +} + +PARAMETER { + v_reset = -65 (mV) + t_ref = 0.1 (ms) + g_reset = 1000 (uS) : large clamp conductance (tau_clamp = C_m/g_reset << dt) +} + +STATE { refrac (ms) } + +ASSIGNED { clamp } + +INITIAL { + refrac = 0 + clamp = 0 +} + +BREAKPOINT { + SOLVE states METHOD cnexp + clamp = 0 + if (refrac > 0) { clamp = 1 } + i = clamp * g_reset * (v - v_reset) +} + +DERIVATIVE states { + refrac' = -1 +} + +POST_EVENT(time) { + refrac = t_ref +} diff --git a/pyNN/arbor/populations.py b/pyNN/arbor/populations.py index 5b762c0b..0936a469 100644 --- a/pyNN/arbor/populations.py +++ b/pyNN/arbor/populations.py @@ -175,11 +175,16 @@ def _set_initial_value_array(self, variable, initial_values): self._lif_initial_v = initial_values return if variable != "v": - warn("todo: handle initial values for ion channel states") - # may have to handle this at the same time as setting parameters - # it is not clear to me if Arbor supports updating decors - # after their creation, other than by set_property - # maybe keep a reference to the return values of arbor.paint? + # Receptor/synapse state variables (e.g. "excitatory.gsyn") and ion + # channel states are initialised to zero by the mechanism's INITIAL + # block; only a non-default request needs handling, which is not yet + # supported. + if "." not in variable: + warn("todo: handle initial values for ion channel states") + # may have to handle this at the same time as setting parameters + # it is not clear to me if Arbor supports updating decors + # after their creation, other than by set_property + # maybe keep a reference to the return values of arbor.paint? return self._arbor_cell_description.base_value.set_initial_values(variable, initial_values) diff --git a/pyNN/arbor/recording.py b/pyNN/arbor/recording.py index ac1565f5..2f3d9014 100644 --- a/pyNN/arbor/recording.py +++ b/pyNN/arbor/recording.py @@ -79,7 +79,9 @@ def _get_arbor_probes(self, gid): probe_index = 0 for variable in self.recorded: if variable.location is None: - pass + # Point neurons (single-compartment cable cells) are recorded + # without an explicit location; default to the soma. + locset = "(root)" else: locset = variable.location diff --git a/pyNN/arbor/standardmodels.py b/pyNN/arbor/standardmodels.py index 5c31c9f2..800bdd09 100644 --- a/pyNN/arbor/standardmodels.py +++ b/pyNN/arbor/standardmodels.py @@ -15,7 +15,7 @@ from ..standardmodels import cells, ion_channels, synapses, electrodes, receptors, build_translations from ..parameters import ParameterSpace, IonicSpecies from ..morphology import Morphology, NeuriteDistribution, LocationGenerator -from .cells import CellDescriptionBuilder +from .cells import CellDescriptionBuilder, PointCellDescriptionBuilder from .simulator import state from .morphology import LabelledLocations @@ -117,7 +117,8 @@ def inject_into(self, cells, location=None): # rename to `locations` ? self.parameter_space.shape = (1,) if location is None: - raise NotImplementedError + # Point neurons (and, by default, any cell) inject at the soma. + location = LabelledLocations("soma") elif isinstance(location, str): location = LabelledLocations(location) elif isinstance(location, LocationGenerator): @@ -361,3 +362,138 @@ class CondExpPostSynapticResponse(receptors.CondExpPostSynapticResponse): model = "expsyn" recordable = ["gsyn"] variable_map = {"gsyn": "g"} + + +class CurrExpPostSynapticResponse(receptors.CurrExpPostSynapticResponse): + + translations = build_translations( + ('locations', 'locations'), + ('tau_syn', 'tau') + ) + model = "expsyn_curr" + recordable = ["isyn"] + variable_map = {"isyn": "isyn"} + + +class LIF(cells.LIF): + __doc__ = cells.LIF.__doc__ + + translations = build_translations( + ('v_rest', 'E_L'), + ('v_reset', 'E_R'), + ('v_thresh', 'V_th'), + ('tau_refrac', 't_ref'), + ('tau_m', 'tau_m'), + ('cm', 'C_m'), + ('i_offset', 'i_offset'), + ) + variable_map = {"v": "v"} + + +class PointNeuron(cells.PointNeuron): + """Composable point neuron, realised as a single-compartment Arbor cable cell. + + Combines a leaky integrate-and-fire ``neuron`` (an :class:`LIF` instance) with + one or more post-synaptic receptors (:class:`CondExpPostSynapticResponse` or + :class:`CurrExpPostSynapticResponse`). See :class:`PointCellDescriptionBuilder` + for how the cable cell is assembled. + """ + + arbor_cell_kind = arbor.cell_kind.cable + + def translate(self, parameters, copy=True): + """Build the Arbor cable-cell description for this point neuron. + + ``parameters`` (the composable parameter space) is not consumed directly: + the neuron and receptor components carry their own (translatable) parameter + spaces, which are assembled into the form the builder expects. + """ + neuron_parameters = self.neuron.native_parameters + post_synaptic_receptors = { + name: (psr.model, psr.native_parameters) + for name, psr in self.post_synaptic_receptors.items() + } + builder = PointCellDescriptionBuilder(neuron_parameters, post_synaptic_receptors) + return ParameterSpace({"cell_description": builder}, schema=None, shape=parameters.shape) + + def reverse_translate(self, native_parameters): + raise NotImplementedError + + def can_record(self, variable, location=None): + return True # todo: implement this properly + + +# The native (LIF) parameter names carried through to PointCellDescriptionBuilder; +# the remaining native names produced by the classic IF models below describe their +# synapses. +_LIF_NATIVE_NAMES = ("E_L", "E_R", "V_th", "t_ref", "tau_m", "C_m", "i_offset") + + +def _point_cell_description(native, receptor_specs, shape): + """Wrap a flat native parameter space (from a classic IF model's base + ``translate()``) into a point-neuron ``cell_description`` ParameterSpace. + + ``receptor_specs`` maps each receptor label to + ``(arbor_synapse_model, {arbor_synapse_param: native_name})``. + """ + neuron_parameters = ParameterSpace( + {name: native[name] for name in _LIF_NATIVE_NAMES}, shape=shape) + post_synaptic_receptors = { + label: (model, + ParameterSpace({arbor_param: native[native_name] + for arbor_param, native_name in param_map.items()}, + shape=shape)) + for label, (model, param_map) in receptor_specs.items() + } + builder = PointCellDescriptionBuilder(neuron_parameters, post_synaptic_receptors) + return ParameterSpace({"cell_description": builder}, schema=None, shape=shape) + + +class IF_curr_exp(cells.IF_curr_exp): + __doc__ = cells.IF_curr_exp.__doc__ + + translations = build_translations( + ('v_rest', 'E_L'), + ('v_reset', 'E_R'), + ('v_thresh', 'V_th'), + ('tau_refrac', 't_ref'), + ('tau_m', 'tau_m'), + ('cm', 'C_m'), + ('i_offset', 'i_offset'), + ('tau_syn_E', 'tau_syn_E'), + ('tau_syn_I', 'tau_syn_I'), + ) + arbor_cell_kind = arbor.cell_kind.cable + + def translate(self, parameters, copy=True): + native = super().translate(parameters, copy) + return _point_cell_description(native, { + "excitatory": ("expsyn_curr", {"tau": "tau_syn_E"}), + "inhibitory": ("expsyn_curr", {"tau": "tau_syn_I"}), + }, parameters.shape) + + +class IF_cond_exp(cells.IF_cond_exp): + __doc__ = cells.IF_cond_exp.__doc__ + + translations = build_translations( + ('v_rest', 'E_L'), + ('v_reset', 'E_R'), + ('v_thresh', 'V_th'), + ('tau_refrac', 't_ref'), + ('tau_m', 'tau_m'), + ('cm', 'C_m'), + ('i_offset', 'i_offset'), + ('tau_syn_E', 'tau_syn_E'), + ('tau_syn_I', 'tau_syn_I'), + ('e_rev_E', 'e_rev_E'), + ('e_rev_I', 'e_rev_I'), + ) + arbor_cell_kind = arbor.cell_kind.cable + + def translate(self, parameters, copy=True): + native = super().translate(parameters, copy) + return _point_cell_description(native, { + "excitatory": ("expsyn", {"tau": "tau_syn_E", "e": "e_rev_E"}), + "inhibitory": ("expsyn", {"tau": "tau_syn_I", "e": "e_rev_I"}), + }, parameters.shape) diff --git a/test/system/scenarios/test__simulation_control.py b/test/system/scenarios/test__simulation_control.py index 1520decd..bc163243 100644 --- a/test/system/scenarios/test__simulation_control.py +++ b/test/system/scenarios/test__simulation_control.py @@ -4,7 +4,7 @@ import pytest -@run_with_simulators("nest", "neuron", "brian2") +@run_with_simulators("nest", "neuron", "brian2", "arbor") def test_reset(sim): """ Run the same simulation n times without recreating the network, @@ -28,7 +28,7 @@ def test_reset(sim): data.segments[0].analogsignals[0], 10) -@run_with_simulators("nest", "neuron", "brian2") +@run_with_simulators("nest", "neuron", "brian2", "arbor") def test_reset_with_clear(sim): """ Run the same simulation n times without recreating the network, @@ -54,7 +54,7 @@ def test_reset_with_clear(sim): data[0].segments[0].analogsignals[0].magnitude, 1e-11) -@run_with_simulators("nest", "neuron", "brian2") +@run_with_simulators("nest", "neuron", "brian2", "arbor") def test_reset_with_spikes(sim): """ Run the same simulation n times without recreating the network, @@ -84,7 +84,7 @@ def test_reset_with_spikes(sim): data.segments[0].analogsignals[0], 10) -@run_with_simulators("nest", "neuron", "brian2") +@run_with_simulators("nest", "neuron", "brian2", "arbor") def test_setup(sim): """ Run the same simulation n times, recreating the network each time, @@ -111,7 +111,7 @@ def test_setup(sim): assert_array_equal(signals[0], data[0].segments[0].analogsignals[0]) -@run_with_simulators("nest", "neuron", "brian2") +@run_with_simulators("nest", "neuron", "brian2", "arbor") def test_run_until(sim): sim.setup(timestep=0.1) diff --git a/test/system/scenarios/test_cell_types.py b/test/system/scenarios/test_cell_types.py index b635fa97..5f6ceded 100644 --- a/test/system/scenarios/test_cell_types.py +++ b/test/system/scenarios/test_cell_types.py @@ -361,6 +361,78 @@ def test_IF_curr_delta_fires_and_resets(sim): sim.end() +def _lif_isi_theory(v_rest, v_reset, v_thresh, tau_m, cm, tau_refrac, i_offset): + """Analytic steady-state ISI of a leaky integrate-and-fire neuron driven by a + constant current i_offset (independent of simulator backend).""" + v_inf = v_rest + i_offset * tau_m / cm + return tau_refrac + tau_m * np.log((v_inf - v_reset) / (v_inf - v_thresh)) + + +@run_with_simulators("arbor", "neuron", "nest", "brian2") +def test_IF_exp_point_neuron_fI(sim): + """IF_cond_exp / IF_curr_exp fire regularly under a constant current, with the + inter-spike interval matching the analytic LIF prediction on every backend. + """ + params = dict(v_rest=-65.0, v_reset=-65.0, v_thresh=-50.0, + tau_m=20.0, cm=1.0, tau_refrac=5.0) + isi_theory = _lif_isi_theory(i_offset=1.0, **params) + for cell_class in (sim.IF_curr_exp, sim.IF_cond_exp): + sim.setup(timestep=0.025) + cell = sim.Population(1, cell_class(i_offset=1.0, tau_syn_E=5.0, tau_syn_I=5.0, + **params), + initial_values={'v': -65.0}) + cell.record('spikes') + sim.run(500.0) + spikes = np.array(cell.get_data().segments[0].spiketrains[0]) + assert len(spikes) >= 10, (cell_class.__name__, len(spikes)) + isi = np.diff(spikes)[1:].mean() # skip the first (from-rest) interval + assert abs(isi - isi_theory) < 1.0, (cell_class.__name__, isi, isi_theory) + sim.end() + + +@run_with_simulators("arbor", "neuron", "nest", "brian2") +def test_IF_curr_exp_EPSP(sim): + """A single presynaptic spike into IF_curr_exp produces an EPSP whose peak + matches the analytic current-based-synapse prediction on every backend.""" + v_rest, cm, tau_m, tau_syn, weight = -65.0, 1.0, 20.0, 5.0, 0.5 + sim.setup(timestep=0.025) + source = sim.Population(1, sim.SpikeSourceArray(spike_times=[20.0])) + cell = sim.Population(1, sim.IF_curr_exp( + v_rest=v_rest, v_reset=v_rest, v_thresh=-50.0, tau_m=tau_m, cm=cm, + tau_refrac=5.0, tau_syn_E=tau_syn, i_offset=0.0), + initial_values={'v': v_rest}) + sim.Projection(source, cell, sim.AllToAllConnector(), + sim.StaticSynapse(weight=weight, delay=1.0), + receptor_type="excitatory") + cell.record('v') + sim.run(100.0) + v = cell.get_data().segments[0].filter(name='v')[0].magnitude[:, 0] + epsp = v.max() - v_rest + # analytic peak of a current-based exponential synapse + t_peak = np.log(tau_m / tau_syn) / (1 / tau_syn - 1 / tau_m) + prefactor = (weight / cm) * (tau_m * tau_syn / (tau_m - tau_syn)) + epsp_theory = prefactor * (np.exp(-t_peak / tau_m) - np.exp(-t_peak / tau_syn)) + assert abs(epsp - epsp_theory) < 0.05, (epsp, epsp_theory) + sim.end() + + +@run_with_simulators("arbor", "neuron", "nest", "brian2") +def test_IF_point_neuron_heterogeneous_current(sim): + """Per-cell i_offset values give per-cell firing rates (guards Arbor's + per-cell cable-cell construction / scalar coercion).""" + sim.setup(timestep=0.025) + cells = sim.Population(3, sim.IF_cond_exp( + v_rest=-65.0, v_reset=-65.0, v_thresh=-50.0, tau_m=20.0, cm=1.0, + tau_refrac=5.0, i_offset=[0.5, 1.0, 1.5]), + initial_values={'v': -65.0}) + cells.record('spikes') + sim.run(500.0) + counts = [len(st) for st in cells.get_data().segments[0].spiketrains] + assert counts[0] == 0, counts # 0.5 nA is sub-threshold + assert 0 < counts[1] < counts[2], counts # firing rate rises with current + sim.end() + + @run_with_simulators("nest", "neuron", "brian2") def test_composed_neuron_model_homogeneous_receptors(sim, plot_figure=False): sim.setup() diff --git a/test/system/scenarios/test_electrodes.py b/test/system/scenarios/test_electrodes.py index 6b130b98..ad1bc07f 100644 --- a/test/system/scenarios/test_electrodes.py +++ b/test/system/scenarios/test_electrodes.py @@ -64,7 +64,7 @@ def test_ticket226(sim): assert v_10p1 > -59.99, v_10p1 -@run_with_simulators("nest", "neuron", "brian2") +@run_with_simulators("nest", "neuron", "brian2", "arbor") def test_issue165(sim): """Ensure that anonymous current sources are not lost.""" sim.setup(timestep=0.1) diff --git a/test/system/scenarios/test_parameter_handling.py b/test/system/scenarios/test_parameter_handling.py index c944d2d6..ea3a3c9a 100644 --- a/test/system/scenarios/test_parameter_handling.py +++ b/test/system/scenarios/test_parameter_handling.py @@ -5,7 +5,7 @@ from .fixtures import run_with_simulators -@run_with_simulators("nest", "neuron", "brian2") +@run_with_simulators("nest", "neuron", "brian2", "arbor") def test_issue241(sim): # "Nest SpikeSourcePoisson populations require all parameters to be passed to constructor" sim.setup() diff --git a/test/system/scenarios/test_recording.py b/test/system/scenarios/test_recording.py index cb9079da..e16113b5 100644 --- a/test/system/scenarios/test_recording.py +++ b/test/system/scenarios/test_recording.py @@ -134,7 +134,7 @@ def test_sampling_interval(sim): sim.end() -@run_with_simulators("nest", "neuron", "brian2") +@run_with_simulators("nest", "neuron", "brian2", "arbor") def test_mix_procedural_and_oo(sim): # cf Issues #217, #234 fn_proc = "test_write_procedural.pkl" diff --git a/test/system/scenarios/test_scenario2.py b/test/system/scenarios/test_scenario2.py index af9823af..2d8ed5ed 100644 --- a/test/system/scenarios/test_scenario2.py +++ b/test/system/scenarios/test_scenario2.py @@ -3,7 +3,7 @@ from .fixtures import run_with_simulators -@run_with_simulators("nest", "neuron", "brian2") +@run_with_simulators("nest", "neuron", "brian2", "arbor") def test_scenario2(sim): """ Array of neurons, each injected with a different current. diff --git a/test/unittests/test_arbor.py b/test/unittests/test_arbor.py index 453b1ec7..2d175f85 100644 --- a/test/unittests/test_arbor.py +++ b/test/unittests/test_arbor.py @@ -155,6 +155,94 @@ def test_if_curr_delta_translation(self): self.assertAlmostEqual(d["i_offset"], 0.0) +@unittest.skipUnless(have_arbor, "Requires Arbor") +class TestPointNeurons(unittest.TestCase): + """Point neurons realised as single-compartment cable cells (Phase 2).""" + + def test_lif_translation(self): + m = arbor_standardmodels.LIF( + v_rest=-65.0, cm=1.0, tau_m=20.0, tau_refrac=2.0, + v_reset=-70.0, v_thresh=-50.0, i_offset=0.1) + native = m.native_parameters + native.shape = (1,) + native.evaluate(simplify=True) + d = native.as_dict() + self.assertAlmostEqual(d["E_L"], -65.0) + self.assertAlmostEqual(d["E_R"], -70.0) + self.assertAlmostEqual(d["V_th"], -50.0) + self.assertAlmostEqual(d["t_ref"], 2.0) + self.assertAlmostEqual(d["tau_m"], 20.0) + self.assertAlmostEqual(d["C_m"], 1.0) # cm passes through in nF + self.assertAlmostEqual(d["i_offset"], 0.1) + + def test_curr_exp_psr(self): + psr = arbor_standardmodels.CurrExpPostSynapticResponse(tau_syn=3.0) + self.assertEqual(psr.model, "expsyn_curr") + self.assertFalse(psr.conductance_based) + native = psr.native_parameters + native.shape = (1,) + native.evaluate(simplify=True) + self.assertAlmostEqual(native.as_dict()["tau"], 3.0) + + def test_cond_exp_psr(self): + psr = arbor_standardmodels.CondExpPostSynapticResponse(tau_syn=3.0, e_syn=-10.0) + self.assertEqual(psr.model, "expsyn") + self.assertTrue(psr.conductance_based) + native = psr.native_parameters + native.shape = (1,) + native.evaluate(simplify=True) + d = native.as_dict() + self.assertAlmostEqual(d["tau"], 3.0) + self.assertAlmostEqual(d["e"], -10.0) + + def test_point_neuron_is_cable_cell(self): + ct = arbor_standardmodels.PointNeuron( + arbor_standardmodels.LIF(), + excitatory=arbor_standardmodels.CondExpPostSynapticResponse(), + inhibitory=arbor_standardmodels.CondExpPostSynapticResponse(e_syn=-70.0), + ) + self.assertIs(ct.arbor_cell_kind, arbor.cell_kind.cable) + self.assertEqual(ct.receptor_types, ["excitatory", "inhibitory"]) + self.assertTrue(ct.conductance_based) + + def test_point_neuron_mixed_receptor_kinds_rejected(self): + ct = arbor_standardmodels.PointNeuron( + arbor_standardmodels.LIF(), + excitatory=arbor_standardmodels.CondExpPostSynapticResponse(), + inhibitory=arbor_standardmodels.CurrExpPostSynapticResponse(), + ) + with self.assertRaises(Exception): + ct.conductance_based + + def test_if_cond_exp_is_cable_cell(self): + ct = arbor_standardmodels.IF_cond_exp(tau_syn_E=1.5, e_rev_E=0.0) + self.assertIs(ct.arbor_cell_kind, arbor.cell_kind.cable) + self.assertTrue(ct.conductance_based) + self.assertEqual(tuple(ct.receptor_types), ("excitatory", "inhibitory")) + + def test_if_curr_exp_is_cable_cell(self): + ct = arbor_standardmodels.IF_curr_exp(tau_syn_E=1.5) + self.assertIs(ct.arbor_cell_kind, arbor.cell_kind.cable) + self.assertFalse(ct.conductance_based) + + def test_if_curr_exp_translation(self): + # the classic model uses the standard flat translate(); check a couple of + # translated native names appear in the resulting synapse parameters + ct = arbor_standardmodels.IF_curr_exp(tau_syn_E=1.5, tau_syn_I=2.5, cm=0.5) + builder = ct.native_parameters["cell_description"].base_value + builder.set_shape((1,)) + exc_model, exc_params = builder.post_synaptic_receptors["excitatory"] + self.assertEqual(exc_model, "expsyn_curr") + self.assertAlmostEqual(exc_params["tau"][0], 1.5) + self.assertAlmostEqual(builder.neuron_parameters["C_m"][0], 0.5) + + def test_reset_and_current_synapse_mechanisms_in_catalogue(self): + cat = arbor.load_catalogue(arbor_simulator.catalogue_path()) + mechs = list(cat) + self.assertIn("lif", mechs) + self.assertIn("expsyn_curr", mechs) + + @unittest.skipUnless(have_arbor, "Requires Arbor") class TestMorphologyLocsets(unittest.TestCase): """Locset generation for recording/placement locations.""" From 0d0daae4df0b1d7e1a1035d606d6d4b0d55258db Mon Sep 17 00:00:00 2001 From: Andrew Davison Date: Thu, 9 Jul 2026 12:01:16 +0200 Subject: [PATCH 3/6] Implement Step/AC/Noisy current sources for the Arbor backend Only DCSource injected current on the Arbor backend; StepCurrentSource, ACSource and NoisyCurrentSource inherited the base inject_into that raises NotImplementedError. This implements all three so every standard current source works on Arbor cable cells (point neurons and multicompartment). Each source is now realised as one or more Arbor iclamp "components", a (envelope, frequency_Hz, phase_deg) tuple where the envelope is a list of (time_ms, amplitude_nA) points. This unifies the two previously divergent current-source loops in the cell-description builders into one shared placement path, and moves all waveform/shape logic onto the source classes (as the NEURON backend does). --- pyNN/arbor/cells.py | 66 +++++----- pyNN/arbor/standardmodels.py | 150 ++++++++++++++++++----- test/system/scenarios/test_electrodes.py | 81 ++++++++++++ test/unittests/test_arbor.py | 83 ++++++++++++- 4 files changed, 314 insertions(+), 66 deletions(-) diff --git a/pyNN/arbor/cells.py b/pyNN/arbor/cells.py index 13504990..d31b4aab 100644 --- a/pyNN/arbor/cells.py +++ b/pyNN/arbor/cells.py @@ -11,14 +11,6 @@ from ..parameters import ParameterSpace -# Units for the parameters of current-source mechanisms (e.g. iclamp). -ELECTRODE_PARAM_UNITS = { - "tstart": U.ms, - "duration": U.ms, - "current": U.nA, -} - - def convert_point(p3d: Point3DWithDiam) -> arbor.mpoint: return arbor.mpoint(p3d.x, p3d.y, p3d.z, p3d.diameter/2) @@ -54,14 +46,42 @@ def set_initial_values(self, variable, initial_values): assert isinstance(initial_values, larray) self.initial_values[variable] = initial_values - def add_current_source(self, model_name, location_generator, index, parameters): + def add_current_source(self, components, location_generator, index): + # ``components`` is a list of (envelope, frequency_Hz, phase_deg) tuples, + # where ``envelope`` is a list of (time_ms, amplitude_nA) points. Each + # standard current source (DC/Step/AC/Noisy) produces its own components; + # see standardmodels.BaseCurrentSource. for i in index: self.current_sources[i].append({ - "model_name": model_name, - "parameters": parameters, + "components": components, "location_generator": location_generator }) + def _make_iclamp(self, component): + """Build an Arbor iclamp mechanism from one (envelope, frequency_Hz, + phase_deg) component. A zero frequency gives a plain (non-sinusoidal) + clamp; a non-zero frequency amplitude-modulates a sine of that envelope.""" + envelope, frequency, phase = component + env = [(t * U.ms, a * U.nA) for (t, a) in envelope] + kwargs = {} + if frequency: + kwargs["frequency"] = frequency * U.Hz + kwargs["phase"] = phase * U.deg + return _compat.get_electrode_mechanism("iclamp")(env, **kwargs) + + def _place_current_sources(self, decor, i, morph): + """Place every current source registered for cell ``i`` onto ``decor``. + ``morph`` is the (evaluated) morphology for locating the injection site, + or ``None`` for point neurons (which inject at the single soma location).""" + for source in self.current_sources[i]: + location_generator = source["location_generator"] + for locset, label in location_generator.generate_locations(morph, label="current_source"): + components = source["components"] + for k, component in enumerate(components): + place_label = label if len(components) == 1 else f"{label}_{k}" + _compat.place_current_source( + decor, locset, self._make_iclamp(component), place_label) + class CellDescriptionBuilder(BaseCellDescriptionBuilder): @@ -163,17 +183,7 @@ def _build_decor(self, i): decor.place(locset, arbor.synapse(pse.model, pse_parameters.as_dict()), label) # insert current sources - for current_source in self.current_sources[i]: - location_generator = current_source["location_generator"] - mechanism = _compat.get_electrode_mechanism(current_source["model_name"]) - for locset, label in location_generator.generate_locations(morph, label=f"{current_source['model_name']}_label"): - params = current_source["parameters"].evaluate(simplify=True).as_dict() - params = { - name: value * ELECTRODE_PARAM_UNITS[name] if name in ELECTRODE_PARAM_UNITS else value - for name, value in params.items() - } - mech = mechanism(**params) - _compat.place_current_source(decor, locset, mech, label) + self._place_current_sources(decor, i, morph) # add spike source decor.place('"root"', arbor.threshold_detector(-10 * U.mV), "detector") @@ -312,15 +322,9 @@ def _build_decor(self, i): self._place_iclamp(decor, "(root)", tstart=0.0, duration=1e12, current=i_offset, label="i_offset") - # injected current sources (e.g. DCSource) - for source in self.current_sources[i]: - location_generator = source["location_generator"] - params = source["parameters"].evaluate(simplify=True).as_dict() - for locset, label in location_generator.generate_locations(None, label=source["model_name"]): - self._place_iclamp( - decor, locset, - tstart=params["tstart"], duration=params["duration"], - current=params["current"], label=label) + # injected current sources (DCSource, StepCurrentSource, ACSource, + # NoisyCurrentSource); a point neuron injects at its single soma location. + self._place_current_sources(decor, i, morph=None) return decor, arbor.label_dict(labels) diff --git a/pyNN/arbor/standardmodels.py b/pyNN/arbor/standardmodels.py index 800bdd09..7ec3e0b0 100644 --- a/pyNN/arbor/standardmodels.py +++ b/pyNN/arbor/standardmodels.py @@ -13,7 +13,7 @@ from arbor import units as U from ..standardmodels import cells, ion_channels, synapses, electrodes, receptors, build_translations -from ..parameters import ParameterSpace, IonicSpecies +from ..parameters import ParameterSpace, IonicSpecies, Sequence from ..morphology import Morphology, NeuriteDistribution, LocationGenerator from .cells import CellDescriptionBuilder, PointCellDescriptionBuilder from .simulator import state @@ -77,44 +77,81 @@ class IF_curr_delta(cells.IF_curr_delta): } -class BaseCurrentSource(object): - pass +# --- current sources --------------------------------------------------------- +# +# Every standard current source is realised as one or more Arbor iclamp +# "components", each a (envelope, frequency_Hz, phase_deg) tuple where ``envelope`` +# is a list of (time_ms, amplitude_nA) points (see cells.BaseCellDescriptionBuilder). +# Arbor iclamp envelope semantics: the current is 0 before the first point, held at +# the last amplitude after the last point, linearly interpolated between points, and +# steps discontinuously where two points share a time. With a non-zero frequency the +# envelope amplitude-modulates a sine, sin(2*pi*f*t + phase), referenced to t=0. +_CURRENT_STOP_SENTINEL = 1e11 # PyNN's "run to the end" stop default is 1e12 -class DCSource(BaseCurrentSource, electrodes.DCSource): - __doc__ = electrodes.DCSource.__doc__ - translations = build_translations( - ('amplitude', 'current'), - ('start', 'tstart'), - ('stop', 'duration', "stop - start", "tstart + duration") - ) +def _as_array(value): + """Coerce a (possibly Sequence-wrapped) parameter to a float ndarray.""" + return np.asarray(value.value if isinstance(value, Sequence) else value, dtype=float) + + +def _box_envelope(start, stop, amplitude): + """A rectangular pulse of ``amplitude`` over [start, stop], zero outside.""" + return [(start, amplitude), (stop, amplitude), (stop, 0.0)] + + +def _staircase_envelope(times, amplitudes): + """A piecewise-constant staircase: ``amplitudes[k]`` is held from ``times[k]`` + until ``times[k+1]`` (and the last value until the end of the run). The current + is zero before ``times[0]``. Duplicated timestamps make the steps discontinuous.""" + envelope = [] + for k in range(len(times)): + if k > 0: + envelope.append((times[k], amplitudes[k - 1])) + envelope.append((times[k], amplitudes[k])) + return envelope + + +def _check_step_times(times): + """Validate StepCurrentSource times (mirrors the checks in the NEURON backend).""" + if not (times >= 0.0).all(): + raise ValueError("Step current cannot accept negative timestamps.") + if not (np.diff(times) > 0.0).all(): + raise ValueError("Step current timestamps should be monotonically increasing.") + + +class BaseCurrentSource(object): + """Base class for the Arbor current sources. + + Subclasses implement :meth:`_iclamp_components`, returning the list of + (envelope, frequency_Hz, phase_deg) components for the injected current; this + base handles resolving the injection target and registering the components on + the target cells' description builder. + """ def inject_into(self, cells, location=None): # rename to `locations` ? - # Native lif cells (IF_curr_delta) have no decor and cannot take an - # i_clamp, so current injection is impossible for them. if hasattr(cells, "parent"): target_pop = cells.parent + cell_descr = target_pop._arbor_cell_description.base_value + index = target_pop.id_to_index(cells.all_cells.astype(int)) elif hasattr(cells, "_arbor_cell_description"): target_pop = cells - else: - target_pop = cells[0].parent - if target_pop.celltype.arbor_cell_kind == arbor.cell_kind.lif: - raise NotImplementedError( - "Current injection into Arbor's native lif_cell (IF_curr_delta) " - "is not supported; use the cable-cell IF models instead.") - if hasattr(cells, "parent"): - cell_descr = cells.parent._arbor_cell_description.base_value - index = cells.parent.id_to_index(cells.all_cells.astype(int)) - elif hasattr(cells, "_arbor_cell_description"): cell_descr = cells._arbor_cell_description.base_value index = cells.id_to_index(cells.all_cells.astype(int)) else: assert isinstance(cells, (list, tuple)) # we're assuming all cells have the same parent here - cell_descr = cells[0].parent._arbor_cell_description.base_value + target_pop = cells[0].parent + cell_descr = target_pop._arbor_cell_description.base_value index = np.array(cells, dtype=int) + # Native lif cells (IF_curr_delta) have no decor and cannot take an + # i_clamp, so current injection is impossible for them. + if target_pop.celltype.arbor_cell_kind == arbor.cell_kind.lif: + raise NotImplementedError( + "Current injection into Arbor's native lif_cell (IF_curr_delta) " + "is not supported; use the cable-cell IF models instead.") + self.parameter_space.shape = (1,) if location is None: # Point neurons (and, by default, any cell) inject at the soma. @@ -122,20 +159,41 @@ def inject_into(self, cells, location=None): # rename to `locations` ? elif isinstance(location, str): location = LabelledLocations(location) elif isinstance(location, LocationGenerator): - # morphology = cells._arbor_cell_description.base_value.parameters["morphology"].base_value # todo: evaluate lazyarray - # locations = location.generate_locations(morphology, label="dc_current_source") - # assert len(locations) == 1 - # locset = locations[0] pass else: raise TypeError("location must be a string or a LocationGenerator") + cell_descr.add_current_source( - model_name="iclamp", + components=self._iclamp_components(), location_generator=location, index=index, - parameters=self.native_parameters ) + def _native_parameters(self): + """The source's native parameters as a plain {name: scalar} dict.""" + native = self.native_parameters + native.shape = (1,) + native.evaluate(simplify=True) + return native.as_dict() + + def _iclamp_components(self): + raise NotImplementedError("Should be redefined in the individual current sources") + + +class DCSource(BaseCurrentSource, electrodes.DCSource): + __doc__ = electrodes.DCSource.__doc__ + + translations = build_translations( + ('amplitude', 'current'), + ('start', 'tstart'), + ('stop', 'duration', "stop - start", "tstart + duration") + ) + + def _iclamp_components(self): + p = self._native_parameters() + start = p["tstart"] + return [(_box_envelope(start, start + p["duration"], p["current"]), 0.0, 0.0)] + class StepCurrentSource(BaseCurrentSource, electrodes.StepCurrentSource): __doc__ = electrodes.StepCurrentSource.__doc__ @@ -145,6 +203,13 @@ class StepCurrentSource(BaseCurrentSource, electrodes.StepCurrentSource): ('times', 'times') ) + def _iclamp_components(self): + p = self._native_parameters() + times = _as_array(p["times"]) + amplitudes = _as_array(p["amplitudes"]) + _check_step_times(times) + return [(_staircase_envelope(times, amplitudes), 0.0, 0.0)] + class ACSource(BaseCurrentSource, electrodes.ACSource): __doc__ = electrodes.ACSource.__doc__ @@ -158,8 +223,21 @@ class ACSource(BaseCurrentSource, electrodes.ACSource): ('phase', 'phase') ) + def _iclamp_components(self): + p = self._native_parameters() + start, stop = p["start"], p["stop"] + # Arbor references the sine to t=0, so shift the phase to make PyNN's + # ``phase`` hold at ``start`` instead. + phase = p["phase"] - 360.0 * p["frequency"] * start / 1000.0 + components = [(_box_envelope(start, stop, p["amplitude"]), p["frequency"], phase)] + if p["offset"] != 0.0: + # Arbor sums co-located clamps, so the DC offset is a separate clamp. + components.append((_box_envelope(start, stop, p["offset"]), 0.0, 0.0)) + return components + class NoisyCurrentSource(BaseCurrentSource, electrodes.NoisyCurrentSource): + __doc__ = electrodes.NoisyCurrentSource.__doc__ translations = build_translations( ('mean', 'mean'), @@ -169,6 +247,20 @@ class NoisyCurrentSource(BaseCurrentSource, electrodes.NoisyCurrentSource): ('dt', 'dt') ) + def _iclamp_components(self): + p = self._native_parameters() + start, stop = p["start"], p["stop"] + if stop >= _CURRENT_STOP_SENTINEL: + raise ValueError( + "NoisyCurrentSource on the Arbor backend must be given a finite `stop` " + "(the noise is precomputed as a per-sample current envelope).") + dt = max(p["dt"], state.dt) + n = int(round((stop - start) / dt)) + times = np.append(start + dt * np.arange(n), stop) + amplitudes = p["mean"] + p["stdev"] * np.random.randn(len(times)) + amplitudes[-1] = 0.0 # switch the current off at `stop` + return [(_staircase_envelope(times, amplitudes), 0.0, 0.0)] + class StaticSynapse(synapses.StaticSynapse): __doc__ = synapses.StaticSynapse.__doc__ diff --git a/test/system/scenarios/test_electrodes.py b/test/system/scenarios/test_electrodes.py index ad1bc07f..d290d408 100644 --- a/test/system/scenarios/test_electrodes.py +++ b/test/system/scenarios/test_electrodes.py @@ -78,6 +78,84 @@ def test_issue165(sim): assert data[150, 0] > -65.0 +@run_with_simulators("arbor", "nest", "neuron", "brian2") +def test_step_current_source(sim): + """A subthreshold StepCurrentSource drives IF_curr_exp: v is flat before the + first step, depolarises while the step is on, and decays once it returns to 0.""" + dt = 0.1 + sim.setup(timestep=dt, min_delay=dt) + v_rest = -60.0 + cell = sim.Population(1, sim.IF_curr_exp(v_rest=v_rest, v_thresh=-50.0, + tau_refrac=5.0, tau_m=20.0, cm=1.0)) + cell.initialize(v=v_rest) + cell.inject(sim.StepCurrentSource(times=[20.0, 50.0], amplitudes=[0.2, 0.0])) + cell.record('v') + sim.run(80.0) + v = cell.get_data().segments[0].filter(name="v")[0] + sim.end() + t = v.times.magnitude + vm = v.magnitude[:, 0] + # no current well before the first step time -> v stays at rest + assert np.allclose(vm[t < 19.0], v_rest, atol=1e-6) + # while the 0.2 nA step is on, v depolarises clearly above rest + assert vm[int(45.0 / dt)] > v_rest + 0.5 + # after the step returns to zero, v decays back towards rest + assert vm[int(75.0 / dt)] < vm[int(50.0 / dt)] + + +@run_with_simulators("arbor", "nest", "neuron", "brian2") +def test_ac_current_source(sim): + """A subthreshold ACSource drives IF_curr_exp: v is flat before start, is + perturbed during injection, and decays monotonically after stop.""" + dt = 0.1 + sim.setup(timestep=dt, min_delay=dt) + v_rest = -60.0 + cell = sim.Population(1, sim.IF_curr_exp(v_rest=v_rest, v_thresh=-50.0, + tau_refrac=5.0, tau_m=20.0, cm=1.0)) + cell.initialize(v=v_rest) + cell.inject(sim.ACSource(start=20.0, stop=60.0, amplitude=0.5, offset=0.2, + frequency=50.0, phase=0.0)) + cell.record('v') + sim.run(90.0) + v = cell.get_data().segments[0].filter(name="v")[0] + sim.end() + t = v.times.magnitude + vm = v.magnitude[:, 0] + # no current well before start -> v stays at rest + assert np.allclose(vm[t < 19.0], v_rest, atol=1e-6) + # during injection v is perturbed away from rest + assert abs(vm[int(40.0 / dt)] - v_rest) > 0.5 + # after stop, v decays monotonically back towards rest + post = vm[int(60.0 / dt):] + assert all(a >= b - 1e-9 for a, b in zip(post, post[1:])) + + +@run_with_simulators("arbor", "nest", "neuron", "brian2") +def test_noisy_current_source(sim): + """A NoisyCurrentSource drives IF_curr_exp: v is flat before start, is + perturbed during injection, and decays monotonically after stop.""" + dt = 0.1 + sim.setup(timestep=dt, min_delay=dt) + v_rest = -60.0 + cell = sim.Population(1, sim.IF_curr_exp(v_rest=v_rest, v_thresh=-50.0, + tau_refrac=5.0, tau_m=20.0, cm=1.0)) + cell.initialize(v=v_rest) + cell.inject(sim.NoisyCurrentSource(mean=0.5, stdev=0.1, start=20.0, stop=60.0, dt=dt)) + cell.record('v') + sim.run(90.0) + v = cell.get_data().segments[0].filter(name="v")[0] + sim.end() + t = v.times.magnitude + vm = v.magnitude[:, 0] + # no current well before start -> v stays at rest + assert np.allclose(vm[t < 19.0], v_rest, atol=1e-6) + # during injection v is perturbed away from rest + assert abs(vm[int(50.0 / dt)] - v_rest) > 0.5 + # after stop, v decays monotonically back towards rest + post = vm[int(60.0 / dt):] + assert all(a >= b - 1e-9 for a, b in zip(post, post[1:])) + + @run_with_simulators("nest", "neuron", "brian2") def test_issue321(sim): """Check that non-zero currents at t=0 are taken into account.""" @@ -695,6 +773,9 @@ def test_issue759(sim): test_changing_electrode(sim) test_ticket226(sim) test_issue165(sim) + test_step_current_source(sim) + test_ac_current_source(sim) + test_noisy_current_source(sim) test_issue321(sim) test_issue437(sim) test_issue442(sim) diff --git a/test/unittests/test_arbor.py b/test/unittests/test_arbor.py index 2d175f85..ea527761 100644 --- a/test/unittests/test_arbor.py +++ b/test/unittests/test_arbor.py @@ -87,12 +87,6 @@ def test_spike_source_array_units(self): self.assertEqual(list(units.keys()), ["times"]) self.assertIs(units["times"], U.ms) - def test_electrode_param_units(self): - units = arbor_cells.ELECTRODE_PARAM_UNITS - self.assertIs(units["tstart"], U.ms) - self.assertIs(units["duration"], U.ms) - self.assertIs(units["current"], U.nA) - def test_if_curr_delta_lif_param_units(self): units = arbor_standardmodels.IF_curr_delta.lif_param_units self.assertIs(units["E_L"], U.mV) @@ -243,6 +237,83 @@ def test_reset_and_current_synapse_mechanisms_in_catalogue(self): self.assertIn("expsyn_curr", mechs) +@unittest.skipUnless(have_arbor, "Requires Arbor") +class TestCurrentSources(unittest.TestCase): + """Each standard current source is realised as one or more Arbor iclamp + (envelope, frequency_Hz, phase_deg) components; check they are built correctly.""" + + @staticmethod + def _components(source): + source.parameter_space.shape = (1,) + return source._iclamp_components() + + def test_dcsource_is_a_box(self): + components = self._components( + arbor_standardmodels.DCSource(amplitude=0.5, start=10.0, stop=20.0)) + self.assertEqual(len(components), 1) + envelope, frequency, phase = components[0] + self.assertEqual(frequency, 0.0) + # rectangular pulse: on at start, off at stop + self.assertEqual(envelope, [(10.0, 0.5), (20.0, 0.5), (20.0, 0.0)]) + + def test_stepcurrentsource_is_a_staircase(self): + components = self._components(arbor_standardmodels.StepCurrentSource( + times=[10.0, 15.0, 20.0], amplitudes=[0.1, 0.2, 0.3])) + self.assertEqual(len(components), 1) + envelope, frequency, _ = components[0] + self.assertEqual(frequency, 0.0) + # piecewise-constant with duplicated breakpoints; holds the last value + self.assertEqual( + [(round(t, 6), round(a, 6)) for (t, a) in envelope], + [(10.0, 0.1), (15.0, 0.1), (15.0, 0.2), (20.0, 0.2), (20.0, 0.3)]) + + def test_stepcurrentsource_rejects_bad_times(self): + with self.assertRaises(ValueError): + self._components(arbor_standardmodels.StepCurrentSource( + times=[10.0, -5.0], amplitudes=[0.1, 0.2])) + with self.assertRaises(ValueError): + self._components(arbor_standardmodels.StepCurrentSource( + times=[10.0, 5.0], amplitudes=[0.1, 0.2])) + + def test_acsource_sine_plus_offset(self): + components = self._components(arbor_standardmodels.ACSource( + start=10.0, stop=20.0, amplitude=0.5, offset=0.1, + frequency=100.0, phase=30.0)) + self.assertEqual(len(components), 2) + (sine_env, freq, phase), (offset_env, offset_freq, _) = components + self.assertEqual(freq, 100.0) + # phase shifted so PyNN's 30 deg holds at start=10 (f=100 Hz -> 360 deg/10 ms) + self.assertAlmostEqual(phase, 30.0 - 360.0 * 100.0 * 10.0 / 1000.0) + self.assertEqual(sine_env, [(10.0, 0.5), (20.0, 0.5), (20.0, 0.0)]) + self.assertEqual(offset_freq, 0.0) + self.assertEqual(offset_env, [(10.0, 0.1), (20.0, 0.1), (20.0, 0.0)]) + + def test_acsource_omits_zero_offset(self): + components = self._components(arbor_standardmodels.ACSource( + start=10.0, stop=20.0, amplitude=0.5, offset=0.0, + frequency=100.0, phase=0.0)) + self.assertEqual(len(components), 1) + + def test_noisycurrentsource_samples_and_zeroes_at_stop(self): + import pyNN.arbor as sim + sim.setup(timestep=0.1, min_delay=0.1) + components = self._components(arbor_standardmodels.NoisyCurrentSource( + mean=0.5, stdev=0.05, start=10.0, stop=12.0, dt=0.5)) + self.assertEqual(len(components), 1) + envelope, frequency, _ = components[0] + self.assertEqual(frequency, 0.0) + self.assertAlmostEqual(envelope[0][0], 10.0) # starts at `start` + self.assertAlmostEqual(envelope[-1][0], 12.0) # ends at `stop` + self.assertEqual(envelope[-1][1], 0.0) # current off at `stop` + + def test_noisycurrentsource_requires_finite_stop(self): + import pyNN.arbor as sim + sim.setup(timestep=0.1, min_delay=0.1) + with self.assertRaises(ValueError): + self._components(arbor_standardmodels.NoisyCurrentSource( + mean=0.5, stdev=0.05, start=10.0)) + + @unittest.skipUnless(have_arbor, "Requires Arbor") class TestMorphologyLocsets(unittest.TestCase): """Locset generation for recording/placement locations.""" From d1c6ee50e0cbda037d65f6132040385a9b3b1539 Mon Sep 17 00:00:00 2001 From: Andrew Davison Date: Thu, 9 Jul 2026 14:29:41 +0200 Subject: [PATCH 4/6] Generalise the Arbor point-neuron builder with a dynamics abstraction Prepare for the remaining standard point-neuron cell types, which differ from LIF in the membrane properties, the leak/channel densities, the reset/dynamics mechanism and the spike-detector threshold. --- pyNN/arbor/cells.py | 135 ++++++++++++++++++++++++++--------- pyNN/arbor/standardmodels.py | 25 ++++--- test/unittests/test_arbor.py | 2 +- 3 files changed, 114 insertions(+), 48 deletions(-) diff --git a/pyNN/arbor/cells.py b/pyNN/arbor/cells.py index d31b4aab..15a70186 100644 --- a/pyNN/arbor/cells.py +++ b/pyNN/arbor/cells.py @@ -227,6 +227,90 @@ def __call__(self, i): LIF_RESET_CONDUCTANCE = 1000.0 # uS +class PointNeuronDynamics: + """Neuron-specific decoration of the single-compartment cable cell that + :class:`PointCellDescriptionBuilder` uses to realise a PyNN point neuron. + + A dynamics object is built from the neuron component's native + :class:`ParameterSpace` and knows how to decorate a ``decor`` with the membrane + properties, the leak / channel densities, and the reset-or-dynamics point + process, plus the voltage at which the cell's ``threshold_detector`` fires. + This lets one builder serve every point-neuron kind (LIF, AdExp, Izhikevich, + HH, ...): the kind-specific slots live here, the geometry, receptor placement + and current-source handling stay in the builder. + + ``native_names`` lists the keys a subclass reads from a flat native parameter + space (used by the flat IF_* models via ``_point_cell_description``). + """ + + native_names = () + + def __init__(self, neuron_parameters): + self.neuron_parameters = neuron_parameters + + def set_shape(self, value): + self.neuron_parameters.shape = value + + def _specific_properties(self, cm_nF, tau_m_ms): + """Specific membrane capacitance [uF/cm2] and leak conductance [S/cm2] + reproducing the whole-cell C_m [nF] and g_leak = C_m/tau_m [uS].""" + c_spec = cm_nF * 1e-3 / _POINT_CELL_AREA_CM2 + g_spec = (cm_nF / tau_m_ms) * 1e-6 / _POINT_CELL_AREA_CM2 + return c_spec, g_spec + + def specific_cm(self, i): + """Specific membrane capacitance [uF/cm2] for cell ``i``.""" + raise NotImplementedError + + def set_property(self, decor, i, initial_v): + decor.set_property(Vm=initial_v * U.mV, cm=self.specific_cm(i) * U.uF / U.cm2) + + def paint(self, decor, i): + """Paint density mechanisms (leak / ion channels). May be a no-op.""" + + def place_reset(self, decor, i): + """Place the reset / dynamics point process at ``(root)``. May be a no-op.""" + + def detector_threshold(self, i): + """Membrane voltage [mV] at which the network-spike detector fires.""" + raise NotImplementedError + + def i_offset(self, i): + p = self.neuron_parameters + return p["i_offset"][i] if "i_offset" in p.keys() else 0.0 + + +class LIFDynamics(PointNeuronDynamics): + """Leaky integrate-and-fire dynamics: a ``pas`` leak whose specific cm/g + reproduce the whole-cell C_m/tau_m, and the ``lif`` reset point process + (nmodl/lif.mod) driven by the threshold_detector via POST_EVENT.""" + + native_names = ("E_L", "E_R", "V_th", "t_ref", "tau_m", "C_m", "i_offset") + + def specific_cm(self, i): + c_spec, _ = self._specific_properties( + self.neuron_parameters["C_m"][i], self.neuron_parameters["tau_m"][i]) + return c_spec + + def paint(self, decor, i): + p = self.neuron_parameters + _c_spec, g_spec = self._specific_properties(p["C_m"][i], p["tau_m"][i]) + # e is a GLOBAL parameter of pas, set via the mechanism name + decor.paint("(all)", arbor.density(f"pas/e={p['E_L'][i]}", {"g": g_spec})) + + def place_reset(self, decor, i): + p = self.neuron_parameters + decor.place( + "(root)", + arbor.synapse("lif", {"v_reset": p["E_R"][i], "t_ref": p["t_ref"][i], + "g_reset": LIF_RESET_CONDUCTANCE}), + "lif_reset", + ) + + def detector_threshold(self, i): + return self.neuron_parameters["V_th"][i] + + class PointCellDescriptionBuilder(BaseCellDescriptionBuilder): """Builds an Arbor cable cell that behaves as a PyNN point neuron. @@ -248,33 +332,27 @@ class PointCellDescriptionBuilder(BaseCellDescriptionBuilder): ``add_current_source`` interface mirrors :class:`CellDescriptionBuilder`, so a point-neuron Population reuses the backend's existing cable-cell machinery. - ``neuron_parameters`` is a native :class:`ParameterSpace` with the LIF keys - (E_L, E_R, V_th, t_ref, tau_m, C_m, i_offset). ``post_synaptic_receptors`` maps - each receptor label to an ``(arbor_synapse_model, native_synapse_parameters)`` - pair. This is the form produced both by the composable + ``dynamics`` is a :class:`PointNeuronDynamics` supplying the neuron-specific + decoration (membrane properties, leak/channels, reset process, detector + threshold). ``post_synaptic_receptors`` maps each receptor label to an + ``(arbor_synapse_model, native_synapse_parameters)`` pair. This is the form + produced both by the composable :class:`~pyNN.arbor.standardmodels.PointNeuron` (from its components) and by the - flat IF_cond_exp/IF_curr_exp standard models (from their translations). + flat IF_* standard models (from their translations). """ - def __init__(self, neuron_parameters, post_synaptic_receptors): + def __init__(self, dynamics, post_synaptic_receptors): super().__init__() - self.neuron_parameters = neuron_parameters + self.dynamics = dynamics self.post_synaptic_receptors = post_synaptic_receptors self.shape = None def set_shape(self, value): self.shape = value - self.neuron_parameters.shape = value + self.dynamics.set_shape(value) for (_model, synapse_parameters) in self.post_synaptic_receptors.values(): synapse_parameters.shape = value - def _specific_properties(self, cm_nF, tau_m_ms): - """Specific membrane capacitance [uF/cm2] and leak conductance [S/cm2] - reproducing the whole-cell C_m [nF] and g_leak = C_m/tau_m [uS].""" - c_spec = cm_nF * 1e-3 / _POINT_CELL_AREA_CM2 - g_spec = (cm_nF / tau_m_ms) * 1e-6 / _POINT_CELL_AREA_CM2 - return c_spec, g_spec - def _build_tree(self, i): d = POINT_CELL_DIAMETER_UM tree = arbor.segment_tree() @@ -283,29 +361,18 @@ def _build_tree(self, i): return tree def _build_decor(self, i): - p = self.neuron_parameters - - E_L = p["E_L"][i] - c_spec, g_spec = self._specific_properties(p["C_m"][i], p["tau_m"][i]) labels = { "all": "(all)", "soma": "(tag 1)", "root": "(root)", } decor = arbor.decor() - decor.set_property( - Vm=self.initial_values["v"][i] * U.mV, - cm=c_spec * U.uF / U.cm2, - ) - # sub-threshold leak (e is a GLOBAL parameter of pas, set via the mechanism name) - decor.paint("(all)", arbor.density(f"pas/e={E_L}", {"g": g_spec})) - # integrate-and-fire reset + refractory clamp - decor.place( - "(root)", - arbor.synapse("lif", {"v_reset": p["E_R"][i], "t_ref": p["t_ref"][i], - "g_reset": LIF_RESET_CONDUCTANCE}), - "lif_reset", - ) + # neuron-specific decoration (membrane properties, leak/channels, reset) + self.dynamics.set_property(decor, i, self.initial_values["v"][i]) + self.dynamics.paint(decor, i) + self.dynamics.place_reset(decor, i) # network spike source - decor.place("(root)", arbor.threshold_detector(p["V_th"][i] * U.mV), "detector") + decor.place("(root)", + arbor.threshold_detector(self.dynamics.detector_threshold(i) * U.mV), + "detector") # one synapse per receptor, labelled by receptor name so that # Projection.arbor_connections can match it by receptor_type prefix @@ -317,7 +384,7 @@ def _build_decor(self, i): labels[label] = "(root)" # constant offset current - i_offset = p["i_offset"][i] if "i_offset" in p.keys() else 0.0 + i_offset = self.dynamics.i_offset(i) if i_offset != 0.0: self._place_iclamp(decor, "(root)", tstart=0.0, duration=1e12, current=i_offset, diff --git a/pyNN/arbor/standardmodels.py b/pyNN/arbor/standardmodels.py index 7ec3e0b0..23c6b17a 100644 --- a/pyNN/arbor/standardmodels.py +++ b/pyNN/arbor/standardmodels.py @@ -15,7 +15,7 @@ from ..standardmodels import cells, ion_channels, synapses, electrodes, receptors, build_translations from ..parameters import ParameterSpace, IonicSpecies, Sequence from ..morphology import Morphology, NeuriteDistribution, LocationGenerator -from .cells import CellDescriptionBuilder, PointCellDescriptionBuilder +from .cells import CellDescriptionBuilder, PointCellDescriptionBuilder, LIFDynamics from .simulator import state from .morphology import LabelledLocations @@ -480,6 +480,8 @@ class LIF(cells.LIF): ('i_offset', 'i_offset'), ) variable_map = {"v": "v"} + # the PointNeuronDynamics that realises this neuron component as a cable cell + dynamics_class = LIFDynamics class PointNeuron(cells.PointNeuron): @@ -500,12 +502,12 @@ def translate(self, parameters, copy=True): the neuron and receptor components carry their own (translatable) parameter spaces, which are assembled into the form the builder expects. """ - neuron_parameters = self.neuron.native_parameters + dynamics = self.neuron.dynamics_class(self.neuron.native_parameters) post_synaptic_receptors = { name: (psr.model, psr.native_parameters) for name, psr in self.post_synaptic_receptors.items() } - builder = PointCellDescriptionBuilder(neuron_parameters, post_synaptic_receptors) + builder = PointCellDescriptionBuilder(dynamics, post_synaptic_receptors) return ParameterSpace({"cell_description": builder}, schema=None, shape=parameters.shape) def reverse_translate(self, native_parameters): @@ -515,21 +517,18 @@ def can_record(self, variable, location=None): return True # todo: implement this properly -# The native (LIF) parameter names carried through to PointCellDescriptionBuilder; -# the remaining native names produced by the classic IF models below describe their -# synapses. -_LIF_NATIVE_NAMES = ("E_L", "E_R", "V_th", "t_ref", "tau_m", "C_m", "i_offset") - - -def _point_cell_description(native, receptor_specs, shape): +def _point_cell_description(native, receptor_specs, shape, dynamics_class=LIFDynamics): """Wrap a flat native parameter space (from a classic IF model's base ``translate()``) into a point-neuron ``cell_description`` ParameterSpace. - ``receptor_specs`` maps each receptor label to + ``dynamics_class`` is the :class:`~pyNN.arbor.cells.PointNeuronDynamics` for the + neuron; the neuron parameters it needs are taken from ``native`` by name (its + ``native_names``). ``receptor_specs`` maps each receptor label to ``(arbor_synapse_model, {arbor_synapse_param: native_name})``. """ neuron_parameters = ParameterSpace( - {name: native[name] for name in _LIF_NATIVE_NAMES}, shape=shape) + {name: native[name] for name in dynamics_class.native_names}, shape=shape) + dynamics = dynamics_class(neuron_parameters) post_synaptic_receptors = { label: (model, ParameterSpace({arbor_param: native[native_name] @@ -537,7 +536,7 @@ def _point_cell_description(native, receptor_specs, shape): shape=shape)) for label, (model, param_map) in receptor_specs.items() } - builder = PointCellDescriptionBuilder(neuron_parameters, post_synaptic_receptors) + builder = PointCellDescriptionBuilder(dynamics, post_synaptic_receptors) return ParameterSpace({"cell_description": builder}, schema=None, shape=shape) diff --git a/test/unittests/test_arbor.py b/test/unittests/test_arbor.py index ea527761..6eef0c75 100644 --- a/test/unittests/test_arbor.py +++ b/test/unittests/test_arbor.py @@ -228,7 +228,7 @@ def test_if_curr_exp_translation(self): exc_model, exc_params = builder.post_synaptic_receptors["excitatory"] self.assertEqual(exc_model, "expsyn_curr") self.assertAlmostEqual(exc_params["tau"][0], 1.5) - self.assertAlmostEqual(builder.neuron_parameters["C_m"][0], 0.5) + self.assertAlmostEqual(builder.dynamics.neuron_parameters["C_m"][0], 0.5) def test_reset_and_current_synapse_mechanisms_in_catalogue(self): cat = arbor.load_catalogue(arbor_simulator.catalogue_path()) From a4a7b11550e880ff1faee5bec1514f0a3b69a00f Mon Sep 17 00:00:00 2001 From: Andrew Davison Date: Thu, 9 Jul 2026 17:13:25 +0200 Subject: [PATCH 5/6] Add alpha-synapse and AdExp/EIF point neurons to the Arbor backend Extends the Arbor point-neuron support with five more standard cell types, all realised as single-compartment cable cells via the PointNeuronDynamics abstraction and cross-validated against the NEURON backend. --- pyNN/arbor/cells.py | 38 +++++++ pyNN/arbor/nmodl/adexp.mod | 79 ++++++++++++++ pyNN/arbor/nmodl/alphasyn.mod | 55 ++++++++++ pyNN/arbor/nmodl/alphasyn_curr.mod | 46 ++++++++ pyNN/arbor/standardmodels.py | 127 ++++++++++++++++++++++- test/system/scenarios/test_cell_types.py | 105 +++++++++++++++++++ test/unittests/test_arbor.py | 46 ++++++++ 7 files changed, 495 insertions(+), 1 deletion(-) create mode 100644 pyNN/arbor/nmodl/adexp.mod create mode 100644 pyNN/arbor/nmodl/alphasyn.mod create mode 100644 pyNN/arbor/nmodl/alphasyn_curr.mod diff --git a/pyNN/arbor/cells.py b/pyNN/arbor/cells.py index 15a70186..c6190fd2 100644 --- a/pyNN/arbor/cells.py +++ b/pyNN/arbor/cells.py @@ -311,6 +311,44 @@ def detector_threshold(self, i): return self.neuron_parameters["V_th"][i] +class AdExpDynamics(PointNeuronDynamics): + """Adaptive-exponential (Brette-Gerstner) dynamics: a ``pas`` leak (as for LIF) + plus the ``adexp`` point process (nmodl/adexp.mod), which adds the exponential + spike current and the adaptation current w and does the reset via POST_EVENT. + The detector fires at the spike-detection threshold ``V_spike`` (not the softer + exponential threshold ``V_th``).""" + + native_names = ("E_L", "E_R", "V_spike", "V_th", "delta", "t_ref", + "tau_m", "C_m", "a", "b", "tau_w", "i_offset") + + def specific_cm(self, i): + c_spec, _ = self._specific_properties( + self.neuron_parameters["C_m"][i], self.neuron_parameters["tau_m"][i]) + return c_spec + + def paint(self, decor, i): + p = self.neuron_parameters + _c_spec, g_spec = self._specific_properties(p["C_m"][i], p["tau_m"][i]) + decor.paint("(all)", arbor.density(f"pas/e={p['E_L'][i]}", {"g": g_spec})) + + def place_reset(self, decor, i): + p = self.neuron_parameters + g_leak = p["C_m"][i] / p["tau_m"][i] # whole-cell leak conductance [uS] + decor.place( + "(root)", + arbor.synapse("adexp", { + "v_reset": p["E_R"][i], "t_ref": p["t_ref"][i], + "g_reset": LIF_RESET_CONDUCTANCE, "GL": g_leak, + "delta": p["delta"][i], "vthresh": p["V_th"][i], + "a": p["a"][i], "b": p["b"][i], "tau_w": p["tau_w"][i], + "EL": p["E_L"][i]}), + "adexp_reset", + ) + + def detector_threshold(self, i): + return self.neuron_parameters["V_spike"][i] + + class PointCellDescriptionBuilder(BaseCellDescriptionBuilder): """Builds an Arbor cable cell that behaves as a PyNN point neuron. diff --git a/pyNN/arbor/nmodl/adexp.mod b/pyNN/arbor/nmodl/adexp.mod new file mode 100644 index 00000000..f1668013 --- /dev/null +++ b/pyNN/arbor/nmodl/adexp.mod @@ -0,0 +1,79 @@ +: Adaptive-exponential (Brette-Gerstner) dynamics for building AdExp / EIF point +: neurons as Arbor cable cells. Like lif.mod this is a POINT_PROCESS placed at the +: soma whose reset is driven by the cell's threshold_detector via POST_EVENT; it +: additionally contributes the exponential spike-generating current and the +: adaptation current w. The sub-threshold leak -GL*(v - EL) is a separate `pas` +: density mechanism (with the same GL), exactly as for lif.mod. +: +: Membrane dynamics realised (with C dv/dt = -sum of mechanism currents): +: C dv/dt = -GL(v-EL) + GL*delta*exp((v-vthresh)/delta) - w + I +: tau_w dw/dt = a*(v - EL) - w +: so this mechanism contributes i = GL*delta term (as iexp, negative/inward) + w. +: On a spike (detector crosses v_spike) POST_EVENT starts a refractory countdown and +: increments w by b; while refractory a strong clamp holds v at v_reset (and the +: exponential/adaptation currents are gated off), matching lif.mod. The clamp is +: gated by a voltage-independent 0/1 multiplier so modcc derives the conductance +: correctly (see lif.mod). +NEURON { + POINT_PROCESS adexp + RANGE v_reset, t_ref, g_reset, GL, delta, vthresh, a, b, tau_w, EL + NONSPECIFIC_CURRENT i +} + +UNITS { + (mV) = (millivolt) + (nA) = (nanoamp) + (uS) = (microsiemens) + (ms) = (millisecond) +} + +PARAMETER { + v_reset = -70.6 (mV) + t_ref = 0.1 (ms) + g_reset = 1000 (uS) + GL = 0.03 (uS) : leak conductance (= C_m/tau_m; matches the pas leak) + delta = 2 (mV) : steepness of the exponential + vthresh = -50.4 (mV) : exponential (soft) threshold V_T + a = 0.004 (uS) : subthreshold adaptation conductance + b = 0.0805 (nA) : spike-triggered adaptation increment + tau_w = 144 (ms) : adaptation time constant + EL = -70.6 (mV) : leak reversal (= v_rest) +} + +STATE { + refrac (ms) + w (nA) +} + +ASSIGNED { clamp iexp (nA) } + +INITIAL { + refrac = 0 + w = 0 + clamp = 0 +} + +BREAKPOINT { + SOLVE states METHOD cnexp + clamp = 0 + if (refrac > 0) { clamp = 1 } + iexp = exp_current(v) + i = clamp * g_reset * (v - v_reset) + (1 - clamp) * (iexp + w) +} + +DERIVATIVE states { + refrac' = -1 + w' = (a * (v - EL) - w) / tau_w +} + +FUNCTION exp_current(v (mV)) (nA) { + LOCAL arg + arg = (v - vthresh) / delta + if (arg > 50) { arg = 50 } : guard against overflow before the reset clamps v + exp_current = -GL * delta * exp(arg) +} + +POST_EVENT(time) { + refrac = t_ref + w = w + b +} diff --git a/pyNN/arbor/nmodl/alphasyn.mod b/pyNN/arbor/nmodl/alphasyn.mod new file mode 100644 index 00000000..13b4de1c --- /dev/null +++ b/pyNN/arbor/nmodl/alphasyn.mod @@ -0,0 +1,55 @@ +: Conductance-based synapse with an alpha-function time course, for IF_cond_alpha +: and conductance-based point neurons. An incoming event of weight `w` (uS) produces +: a conductance g(t) = w * (t/tau) * exp(1 - t/tau), peaking at `w` a time `tau` +: after the event; the membrane current is i = g*(v - e). +: +: The alpha shape is the impulse response of the coupled linear system +: g' = z - g/tau +: z' = -z/tau +: with each event adding w*exp(1)/tau to z (so that the peak of g equals w). Unlike +: the NEURON backend's alphasyn.mod this needs no spike queue. The system is a +: non-diagonalisable Jordan block, so `cnexp` cannot solve it; `sparse` (Arbor's +: implicit/backward-Euler solver) is used instead, which is first-order accurate -- +: the alpha peak is slightly under-estimated at large dt (a few % of the synaptic +: current for dt ~ tau), but the downstream membrane response matches the exact +: alpha to <1% for dt <= 0.05 ms. +NEURON { + POINT_PROCESS alphasyn + RANGE tau, e + NONSPECIFIC_CURRENT i +} + +UNITS { + (mV) = (millivolt) + (uS) = (microsiemens) + (ms) = (millisecond) +} + +PARAMETER { + tau = 5.0 (ms) + e = 0 (mV) +} + +STATE { + g (uS) + z (uS/ms) +} + +INITIAL { + g = 0 + z = 0 +} + +BREAKPOINT { + SOLVE state METHOD sparse + i = g * (v - e) +} + +DERIVATIVE state { + g' = z - g/tau + z' = -z/tau +} + +NET_RECEIVE(weight (uS)) { + z = z + weight * exp(1) / tau +} diff --git a/pyNN/arbor/nmodl/alphasyn_curr.mod b/pyNN/arbor/nmodl/alphasyn_curr.mod new file mode 100644 index 00000000..5e8396d7 --- /dev/null +++ b/pyNN/arbor/nmodl/alphasyn_curr.mod @@ -0,0 +1,46 @@ +: Current-based synapse with an alpha-function time course, for IF_curr_alpha and +: current-based point neurons. An incoming event of weight `w` (nA) produces a +: synaptic current isyn(t) = w * (t/tau) * exp(1 - t/tau), peaking at `w` a time +: `tau` after the event. A positive weight injects depolarising (inward) current, +: matching PyNN's convention, so the contributed membrane current is i = -isyn. +: +: Same coupled-ODE realisation of the alpha shape as alphasyn.mod (no spike queue); +: solved with `sparse` (implicit) because the Jordan-block system is not diagonal. +NEURON { + POINT_PROCESS alphasyn_curr + RANGE tau + NONSPECIFIC_CURRENT i +} + +UNITS { + (nA) = (nanoamp) + (ms) = (millisecond) +} + +PARAMETER { + tau = 5.0 (ms) +} + +STATE { + isyn (nA) + z (nA/ms) +} + +INITIAL { + isyn = 0 + z = 0 +} + +BREAKPOINT { + SOLVE state METHOD sparse + i = -isyn +} + +DERIVATIVE state { + isyn' = z - isyn/tau + z' = -z/tau +} + +NET_RECEIVE(weight (nA)) { + z = z + weight * exp(1) / tau +} diff --git a/pyNN/arbor/standardmodels.py b/pyNN/arbor/standardmodels.py index 23c6b17a..31d02f5a 100644 --- a/pyNN/arbor/standardmodels.py +++ b/pyNN/arbor/standardmodels.py @@ -15,7 +15,8 @@ from ..standardmodels import cells, ion_channels, synapses, electrodes, receptors, build_translations from ..parameters import ParameterSpace, IonicSpecies, Sequence from ..morphology import Morphology, NeuriteDistribution, LocationGenerator -from .cells import CellDescriptionBuilder, PointCellDescriptionBuilder, LIFDynamics +from .cells import (CellDescriptionBuilder, PointCellDescriptionBuilder, + LIFDynamics, AdExpDynamics) from .simulator import state from .morphology import LabelledLocations @@ -467,6 +468,18 @@ class CurrExpPostSynapticResponse(receptors.CurrExpPostSynapticResponse): variable_map = {"isyn": "isyn"} +class CondAlphaPostSynapticResponse(receptors.CondAlphaPostSynapticResponse): + + translations = build_translations( + ('locations', 'locations'), + ('e_syn', 'e'), + ('tau_syn', 'tau') + ) + model = "alphasyn" + recordable = ["gsyn"] + variable_map = {"gsyn": "g"} + + class LIF(cells.LIF): __doc__ = cells.LIF.__doc__ @@ -484,6 +497,31 @@ class LIF(cells.LIF): dynamics_class = LIFDynamics +# translations shared by the AdExp neuron component and the flat EIF_* cell types +_ADEXP_TRANSLATIONS = build_translations( + ('v_rest', 'E_L'), + ('v_reset', 'E_R'), + ('v_spike', 'V_spike'), + ('v_thresh', 'V_th'), # the (soft) exponential threshold V_T + ('delta_T', 'delta'), + ('tau_refrac', 't_ref'), + ('tau_m', 'tau_m'), + ('cm', 'C_m'), + ('a', 'a', "a*0.001", "a*1000"), # subthreshold adaptation, nS -> uS + ('b', 'b'), + ('tau_w', 'tau_w'), + ('i_offset', 'i_offset'), +) + + +class AdExp(cells.AdExp): + __doc__ = cells.AdExp.__doc__ + + translations = _ADEXP_TRANSLATIONS + variable_map = {"v": "v", "w": "w"} + dynamics_class = AdExpDynamics + + class PointNeuron(cells.PointNeuron): """Composable point neuron, realised as a single-compartment Arbor cable cell. @@ -588,3 +626,90 @@ def translate(self, parameters, copy=True): "excitatory": ("expsyn", {"tau": "tau_syn_E", "e": "e_rev_E"}), "inhibitory": ("expsyn", {"tau": "tau_syn_I", "e": "e_rev_I"}), }, parameters.shape) + + +class IF_curr_alpha(cells.IF_curr_alpha): + __doc__ = cells.IF_curr_alpha.__doc__ + + translations = build_translations( + ('v_rest', 'E_L'), + ('v_reset', 'E_R'), + ('v_thresh', 'V_th'), + ('tau_refrac', 't_ref'), + ('tau_m', 'tau_m'), + ('cm', 'C_m'), + ('i_offset', 'i_offset'), + ('tau_syn_E', 'tau_syn_E'), + ('tau_syn_I', 'tau_syn_I'), + ) + arbor_cell_kind = arbor.cell_kind.cable + + def translate(self, parameters, copy=True): + native = super().translate(parameters, copy) + return _point_cell_description(native, { + "excitatory": ("alphasyn_curr", {"tau": "tau_syn_E"}), + "inhibitory": ("alphasyn_curr", {"tau": "tau_syn_I"}), + }, parameters.shape) + + +class IF_cond_alpha(cells.IF_cond_alpha): + __doc__ = cells.IF_cond_alpha.__doc__ + + translations = build_translations( + ('v_rest', 'E_L'), + ('v_reset', 'E_R'), + ('v_thresh', 'V_th'), + ('tau_refrac', 't_ref'), + ('tau_m', 'tau_m'), + ('cm', 'C_m'), + ('i_offset', 'i_offset'), + ('tau_syn_E', 'tau_syn_E'), + ('tau_syn_I', 'tau_syn_I'), + ('e_rev_E', 'e_rev_E'), + ('e_rev_I', 'e_rev_I'), + ) + arbor_cell_kind = arbor.cell_kind.cable + + def translate(self, parameters, copy=True): + native = super().translate(parameters, copy) + return _point_cell_description(native, { + "excitatory": ("alphasyn", {"tau": "tau_syn_E", "e": "e_rev_E"}), + "inhibitory": ("alphasyn", {"tau": "tau_syn_I", "e": "e_rev_I"}), + }, parameters.shape) + + +# conductance-synapse translations shared by the EIF_cond_* cell types +_EIF_SYNAPSE_TRANSLATIONS = build_translations( + ('tau_syn_E', 'tau_syn_E'), + ('tau_syn_I', 'tau_syn_I'), + ('e_rev_E', 'e_rev_E'), + ('e_rev_I', 'e_rev_I'), +) + + +class EIF_cond_exp_isfa_ista(cells.EIF_cond_exp_isfa_ista): + __doc__ = cells.EIF_cond_exp_isfa_ista.__doc__ + + translations = {**_ADEXP_TRANSLATIONS, **_EIF_SYNAPSE_TRANSLATIONS} + arbor_cell_kind = arbor.cell_kind.cable + + def translate(self, parameters, copy=True): + native = super().translate(parameters, copy) + return _point_cell_description(native, { + "excitatory": ("expsyn", {"tau": "tau_syn_E", "e": "e_rev_E"}), + "inhibitory": ("expsyn", {"tau": "tau_syn_I", "e": "e_rev_I"}), + }, parameters.shape, dynamics_class=AdExpDynamics) + + +class EIF_cond_alpha_isfa_ista(cells.EIF_cond_alpha_isfa_ista): + __doc__ = cells.EIF_cond_alpha_isfa_ista.__doc__ + + translations = {**_ADEXP_TRANSLATIONS, **_EIF_SYNAPSE_TRANSLATIONS} + arbor_cell_kind = arbor.cell_kind.cable + + def translate(self, parameters, copy=True): + native = super().translate(parameters, copy) + return _point_cell_description(native, { + "excitatory": ("alphasyn", {"tau": "tau_syn_E", "e": "e_rev_E"}), + "inhibitory": ("alphasyn", {"tau": "tau_syn_I", "e": "e_rev_I"}), + }, parameters.shape, dynamics_class=AdExpDynamics) diff --git a/test/system/scenarios/test_cell_types.py b/test/system/scenarios/test_cell_types.py index 5f6ceded..fab2585c 100644 --- a/test/system/scenarios/test_cell_types.py +++ b/test/system/scenarios/test_cell_types.py @@ -416,6 +416,111 @@ def test_IF_curr_exp_EPSP(sim): sim.end() +@run_with_simulators("arbor", "nest", "neuron", "brian2") +def test_EIF_cond_alpha_isfa_ista_spike_times(sim): + """The AdExp dynamics of EIF_cond_alpha_isfa_ista driven by a constant current + reproduce the reference spike times (from the NEURON backend) to <1% on every + backend. (Same setup as test_EIF_cond_alpha_isfa_ista but recording only + spikes, so it can also run on Arbor, which does not yet record `w`.)""" + sim.setup(timestep=0.01, min_delay=0.1, max_delay=4.0) + ifcell = sim.Population(1, sim.EIF_cond_alpha_isfa_ista( + i_offset=1.0, tau_refrac=2.0, v_spike=-40)) + ifcell.initialize(v=-65, w=0) + ifcell.record('spikes') + sim.run(200.0) + spike_times = ifcell.get_data().segments[0].spiketrains[0].rescale(pq.ms).magnitude + sim.end() + expected_spike_times = np.array( + [10.015, 25.515, 43.168, 63.41, 86.649, 113.112, 142.663, 174.76]) + assert len(spike_times) == len(expected_spike_times), (spike_times, expected_spike_times) + diff = (spike_times - expected_spike_times) / expected_spike_times + assert abs(diff).max() < 0.01, abs(diff).max() + + +@run_with_simulators("arbor", "neuron", "nest", "brian2") +def test_EIF_cond_exp_isfa_ista_adaptation(sim): + """A constant supra-threshold current into EIF_cond_exp_isfa_ista produces a + regular spike train with spike-frequency adaptation (monotonically increasing + ISIs), matching the NEURON backend (cross-checked to <0.2% on mean ISI).""" + sim.setup(timestep=0.01, min_delay=0.1) + cell = sim.Population(1, sim.EIF_cond_exp_isfa_ista( + i_offset=1.0, tau_refrac=2.0, v_spike=-40.0)) + cell.initialize(v=-70.6, w=0.0) + cell.record('spikes') + sim.run(200.0) + st = cell.get_data().segments[0].spiketrains[0].magnitude + sim.end() + assert len(st) >= 6, len(st) + isis = np.diff(st) + # spike-frequency adaptation: each ISI is longer than the previous one + assert np.all(np.diff(isis) > 0), isis + # adaptation is substantial (last ISI at least 1.5x the first) + assert isis[-1] > 1.5 * isis[0], isis + + +@run_with_simulators("arbor", "neuron", "nest", "brian2") +def test_IF_curr_alpha_EPSP(sim): + """A single presynaptic spike into IF_curr_alpha produces an alpha-shaped EPSP + whose peak matches the analytic current-based-alpha-synapse prediction on every + backend (Arbor realises the alpha with an implicit solver; cross-checked to + <0.2% of the NEURON/analytic value at this timestep).""" + v_rest, cm, tau_m, tau_syn, weight = -65.0, 1.0, 20.0, 2.0, 0.2 + sim.setup(timestep=0.025) + source = sim.Population(1, sim.SpikeSourceArray(spike_times=[20.0])) + cell = sim.Population(1, sim.IF_curr_alpha( + v_rest=v_rest, v_reset=v_rest, v_thresh=-50.0, tau_m=tau_m, cm=cm, + tau_refrac=5.0, tau_syn_E=tau_syn, i_offset=0.0), + initial_values={'v': v_rest}) + sim.Projection(source, cell, sim.AllToAllConnector(), + sim.StaticSynapse(weight=weight, delay=1.0), + receptor_type="excitatory") + cell.record('v') + sim.run(120.0) + v = cell.get_data().segments[0].filter(name='v')[0].magnitude[:, 0] + epsp = v.max() - v_rest + # analytic peak of a current-based alpha synapse into an RC membrane: + # u(t) = (w*e)/(C*tau_syn) * exp(-t/tau_m) * (exp(k t)(k t - 1) + 1)/k^2 + t = np.arange(0, 300, 0.001) + k = 1 / tau_m - 1 / tau_syn + u = ((weight * np.e) / (cm * tau_syn) * np.exp(-t / tau_m) + * (np.exp(k * t) * (k * t - 1) + 1) / k ** 2) + epsp_theory = u.max() + assert abs(epsp - epsp_theory) < 0.01 * epsp_theory, (epsp, epsp_theory) + sim.end() + + +@run_with_simulators("arbor", "neuron", "nest", "brian2") +def test_IF_cond_alpha_EPSP(sim): + """A single presynaptic spike into IF_cond_alpha produces an alpha-shaped + (rise-then-decay) depolarising EPSP peaking a few ms after the input.""" + v_rest, tau_syn, weight = -65.0, 2.0, 0.05 + onset, delay = 20.0, 1.0 + sim.setup(timestep=0.025) + source = sim.Population(1, sim.SpikeSourceArray(spike_times=[onset])) + cell = sim.Population(1, sim.IF_cond_alpha( + v_rest=v_rest, v_reset=v_rest, v_thresh=-50.0, tau_m=20.0, cm=1.0, + tau_refrac=5.0, tau_syn_E=tau_syn, e_rev_E=0.0), + initial_values={'v': v_rest}) + sim.Projection(source, cell, sim.AllToAllConnector(), + sim.StaticSynapse(weight=weight, delay=delay), + receptor_type="excitatory") + cell.record('v') + sim.run(120.0) + signal = cell.get_data().segments[0].filter(name='v')[0] + v = signal.magnitude[:, 0] + t = signal.times.magnitude + # depolarising and alpha-shaped: quiescent before the input, single peak after + assert np.allclose(v[t < onset], v_rest, atol=1e-6) + i_peak = int(np.argmax(v)) + assert v[i_peak] - v_rest > 0.5 + # peak occurs after the synaptic peak (onset + delay + tau_syn) and before the + # membrane time constant washes it out + assert onset + delay + tau_syn < t[i_peak] < onset + delay + 20.0 + # decays monotonically back towards rest after the peak + assert v[-1] < v[i_peak] + sim.end() + + @run_with_simulators("arbor", "neuron", "nest", "brian2") def test_IF_point_neuron_heterogeneous_current(sim): """Per-cell i_offset values give per-cell firing rates (guards Arbor's diff --git a/test/unittests/test_arbor.py b/test/unittests/test_arbor.py index 6eef0c75..2be63188 100644 --- a/test/unittests/test_arbor.py +++ b/test/unittests/test_arbor.py @@ -230,11 +230,57 @@ def test_if_curr_exp_translation(self): self.assertAlmostEqual(exc_params["tau"][0], 1.5) self.assertAlmostEqual(builder.dynamics.neuron_parameters["C_m"][0], 0.5) + def test_if_curr_alpha_is_cable_cell(self): + ct = arbor_standardmodels.IF_curr_alpha(tau_syn_E=1.5) + self.assertIs(ct.arbor_cell_kind, arbor.cell_kind.cable) + self.assertFalse(ct.conductance_based) + + def test_if_curr_alpha_translation(self): + ct = arbor_standardmodels.IF_curr_alpha(tau_syn_E=1.5, tau_syn_I=2.5) + builder = ct.native_parameters["cell_description"].base_value + builder.set_shape((1,)) + exc_model, exc_params = builder.post_synaptic_receptors["excitatory"] + self.assertEqual(exc_model, "alphasyn_curr") + self.assertAlmostEqual(exc_params["tau"][0], 1.5) + + def test_if_cond_alpha_translation(self): + ct = arbor_standardmodels.IF_cond_alpha(tau_syn_E=1.5, e_rev_E=0.0) + self.assertIs(ct.arbor_cell_kind, arbor.cell_kind.cable) + self.assertTrue(ct.conductance_based) + builder = ct.native_parameters["cell_description"].base_value + builder.set_shape((1,)) + exc_model, exc_params = builder.post_synaptic_receptors["excitatory"] + self.assertEqual(exc_model, "alphasyn") + self.assertAlmostEqual(exc_params["tau"][0], 1.5) + self.assertAlmostEqual(exc_params["e"][0], 0.0) + + def test_eif_cond_exp_is_cable_cell(self): + ct = arbor_standardmodels.EIF_cond_exp_isfa_ista() + self.assertIs(ct.arbor_cell_kind, arbor.cell_kind.cable) + self.assertTrue(ct.conductance_based) + + def test_eif_cond_exp_translation(self): + # a (subthreshold adaptation) is translated nS -> uS; the dynamics is AdExp + ct = arbor_standardmodels.EIF_cond_exp_isfa_ista(a=4.0, tau_syn_E=1.5) + builder = ct.native_parameters["cell_description"].base_value + builder.set_shape((1,)) + self.assertIsInstance(builder.dynamics, arbor_cells.AdExpDynamics) + self.assertAlmostEqual(builder.dynamics.neuron_parameters["a"][0], 0.004) # 4 nS -> uS + exc_model, exc_params = builder.post_synaptic_receptors["excitatory"] + self.assertEqual(exc_model, "expsyn") + self.assertAlmostEqual(exc_params["tau"][0], 1.5) + + def test_adexp_component_uses_adexp_dynamics(self): + self.assertIs(arbor_standardmodels.AdExp.dynamics_class, arbor_cells.AdExpDynamics) + def test_reset_and_current_synapse_mechanisms_in_catalogue(self): cat = arbor.load_catalogue(arbor_simulator.catalogue_path()) mechs = list(cat) self.assertIn("lif", mechs) self.assertIn("expsyn_curr", mechs) + self.assertIn("alphasyn", mechs) + self.assertIn("alphasyn_curr", mechs) + self.assertIn("adexp", mechs) @unittest.skipUnless(have_arbor, "Requires Arbor") From e6f6a35e9ae156230da35e1a9519f4200cbc2515 Mon Sep 17 00:00:00 2001 From: Andrew Davison Date: Thu, 9 Jul 2026 17:14:49 +0200 Subject: [PATCH 6/6] Add gsfa_grr, Izhikevich and HH_cond_exp point neurons to the Arbor backend Completes the tractable missing standard cell types, all as single-compartment cable cells via PointNeuronDynamics, cross-validated against NEURON. --- pyNN/arbor/cells.py | 85 +++++++++++++++++++++++ pyNN/arbor/nmodl/hh_traub.mod | 64 +++++++++++++++++ pyNN/arbor/nmodl/izhikevich.mod | 75 ++++++++++++++++++++ pyNN/arbor/nmodl/lif_gsfa_grr.mod | 73 ++++++++++++++++++++ pyNN/arbor/standardmodels.py | 88 +++++++++++++++++++++++- test/system/scenarios/test_cell_types.py | 47 ++++++++++++- test/unittests/test_arbor.py | 46 +++++++++++++ 7 files changed, 476 insertions(+), 2 deletions(-) create mode 100644 pyNN/arbor/nmodl/hh_traub.mod create mode 100644 pyNN/arbor/nmodl/izhikevich.mod create mode 100644 pyNN/arbor/nmodl/lif_gsfa_grr.mod diff --git a/pyNN/arbor/cells.py b/pyNN/arbor/cells.py index c6190fd2..0e186783 100644 --- a/pyNN/arbor/cells.py +++ b/pyNN/arbor/cells.py @@ -349,6 +349,91 @@ def detector_threshold(self, i): return self.neuron_parameters["V_spike"][i] +class GsfaGrrDynamics(LIFDynamics): + """LIF dynamics plus conductance-based spike-frequency adaptation and a + relative-refractory mechanism (Muller 2007), via the ``lif_gsfa_grr`` point + process (nmodl/lif_gsfa_grr.mod). Everything else (pas leak, detector at V_th) + is inherited from :class:`LIFDynamics`.""" + + native_names = LIFDynamics.native_names + ( + "E_s", "E_r", "tau_s", "tau_r", "q_s", "q_r") + + def place_reset(self, decor, i): + p = self.neuron_parameters + decor.place( + "(root)", + arbor.synapse("lif_gsfa_grr", { + "v_reset": p["E_R"][i], "t_ref": p["t_ref"][i], + "g_reset": LIF_RESET_CONDUCTANCE, + "E_s": p["E_s"][i], "E_r": p["E_r"][i], + "tau_s": p["tau_s"][i], "tau_r": p["tau_r"][i], + "q_s": p["q_s"][i], "q_r": p["q_r"][i]}), + "gsfa_grr_reset", + ) + + +class HHDynamics(PointNeuronDynamics): + """Hodgkin-Huxley (Traub) dynamics: the ``hh_traub`` density mechanism + (nmodl/hh_traub.mod) with whole-cell Na/K/leak conductances converted to + specific S/cm2, and the Na/K reversal potentials set on the cell's ions. HH is + a genuine spiking model, so there is no reset process; the network-spike + detector fires at 10 mV (matching the NEURON backend's HH threshold).""" + + native_names = ("gnabar", "gkbar", "gl", "ena", "ek", "el", "vT", + "C_m", "i_offset") + + def specific_cm(self, i): + return self.neuron_parameters["C_m"][i] * 1e-3 / _POINT_CELL_AREA_CM2 + + def _specific_g(self, g_uS): + """Whole-cell conductance [uS] -> specific conductance [S/cm2].""" + return g_uS * 1e-6 / _POINT_CELL_AREA_CM2 + + def paint(self, decor, i): + p = self.neuron_parameters + decor.paint("(all)", arbor.density("hh_traub", { + "gnabar": self._specific_g(p["gnabar"][i]), + "gkbar": self._specific_g(p["gkbar"][i]), + "gl": self._specific_g(p["gl"][i]), + "el": p["el"][i], "vT": p["vT"][i]})) + decor.set_ion("na", rev_pot=p["ena"][i] * U.mV) + decor.set_ion("k", rev_pot=p["ek"][i] * U.mV) + + def detector_threshold(self, i): + return 10.0 + + +IZHIKEVICH_CM_NF = 0.001 # fixed membrane capacitance of the Izhikevich model [nF] +IZHIKEVICH_VTHRESH_MV = 30.0 # spike / reset threshold [mV] + + +class IzhikevichDynamics(PointNeuronDynamics): + """Izhikevich quadratic-integrate-and-fire dynamics (nmodl/izhikevich.mod): no + leak, a fixed capacitance Cm, and the reset (v -> c, u += d) done via the + detector + POST_EVENT one-step clamp. The membrane current is entirely supplied + by the izhikevich point process; the injected/offset current enters as an + iclamp (contributing I/Cm to dv/dt).""" + + native_names = ("a", "b", "c", "d", "i_offset") + + def specific_cm(self, i): + return IZHIKEVICH_CM_NF * 1e-3 / _POINT_CELL_AREA_CM2 + + def place_reset(self, decor, i): + p = self.neuron_parameters + decor.place( + "(root)", + arbor.synapse("izhikevich", { + "a": p["a"][i], "b": p["b"][i], "c": p["c"][i], "d": p["d"][i], + "Cm": IZHIKEVICH_CM_NF, "vthresh": IZHIKEVICH_VTHRESH_MV, + "g_reset": LIF_RESET_CONDUCTANCE}), + "izhikevich", + ) + + def detector_threshold(self, i): + return IZHIKEVICH_VTHRESH_MV + + class PointCellDescriptionBuilder(BaseCellDescriptionBuilder): """Builds an Arbor cable cell that behaves as a PyNN point neuron. diff --git a/pyNN/arbor/nmodl/hh_traub.mod b/pyNN/arbor/nmodl/hh_traub.mod new file mode 100644 index 00000000..d2c6d1d0 --- /dev/null +++ b/pyNN/arbor/nmodl/hh_traub.mod @@ -0,0 +1,64 @@ +: Traub-modified Hodgkin-Huxley channels, for building HH_cond_exp as an Arbor +: single-compartment cable cell. A near-verbatim port of the NEURON backend's +: hh_traub.mod: the TABLE statement is dropped (Arbor's modcc has no TABLE; the +: rates are computed exactly each step) and the rate "trap" is written with Arbor's +: exprelr (vtrap(x,y) = y*exprelr(x/y) = x/(exp(x/y)-1), singularity-safe). The +: NEST-compatible initial state m = h = n = 0 is kept. Reversal potentials for Na/K +: come from the ions (set per cell); the leak is a NONSPECIFIC_CURRENT. +NEURON { + SUFFIX hh_traub + USEION na READ ena WRITE ina + USEION k READ ek WRITE ik + NONSPECIFIC_CURRENT il + RANGE gnabar, gkbar, gl, el, vT +} + +UNITS { + (mV) = (millivolt) + (S) = (siemens) +} + +PARAMETER { + gnabar = 0.02 (S/cm2) + gkbar = 0.006 (S/cm2) + gl = 0.00001 (S/cm2) + el = -60.0 (mV) + vT = -63.0 (mV) +} + +STATE { m h n } + +BREAKPOINT { + SOLVE states METHOD cnexp + ina = gnabar * m * m * m * h * (v - ena) + ik = gkbar * n * n * n * n * (v - ek) + il = gl * (v - el) +} + +INITIAL { + : for compatibility with NEST, start from m = h = n = 0 + m = 0 + h = 0 + n = 0 +} + +DERIVATIVE states { + LOCAL u, alpha, beta + u = v - vT + : sodium activation + alpha = 0.32 * vtrap(13 - u, 4) + beta = 0.28 * vtrap(u - 40, 5) + m' = alpha - m * (alpha + beta) + : sodium inactivation + alpha = 0.128 * exp((17 - u) / 18) + beta = 4 / (exp((40 - u) / 5) + 1) + h' = alpha - h * (alpha + beta) + : potassium activation + alpha = 0.032 * vtrap(15 - u, 5) + beta = 0.5 * exp((10 - u) / 40) + n' = alpha - n * (alpha + beta) +} + +FUNCTION vtrap(x, y) { + vtrap = y * exprelr(x / y) +} diff --git a/pyNN/arbor/nmodl/izhikevich.mod b/pyNN/arbor/nmodl/izhikevich.mod new file mode 100644 index 00000000..d3b8674e --- /dev/null +++ b/pyNN/arbor/nmodl/izhikevich.mod @@ -0,0 +1,75 @@ +: Izhikevich (2003) quadratic integrate-and-fire neuron, for building the PyNN +: Izhikevich cell type as an Arbor cable cell: +: dv/dt = 0.04 v^2 + 5 v + 140 - u + I +: du/dt = a (b v - u) +: with reset, when v reaches vthresh: v -> c, u -> u + d. +: +: v is integrated by Arbor's cable solver: with the cell's absolute capacitance set +: equal to Cm, this POINT_PROCESS supplies i = -Cm(0.04 v^2 + 5 v + 140 - u) so that +: Cm dv/dt = -i reproduces the Izhikevich dv/dt (the injected current I enters via a +: separate iclamp for i_offset, contributing I/Cm as in the NEURON backend). There +: is no leak mechanism. u is a STATE solved by cnexp. +: +: Arbor mechanisms cannot write v, so the reset (v -> c) is done as in lif.mod: the +: threshold_detector (set to vthresh) delivers a POST_EVENT that starts a very short +: countdown during which a strong clamp holds v at c; the countdown t_reset is tiny +: (<= one timestep) so there is no artificial refractory period, only a one-step +: reset. u += d is applied in the POST_EVENT. The clamp is gated by the +: voltage-independent 0/1 multiplier `clamp` so modcc derives the conductance +: correctly (see lif.mod). +NEURON { + POINT_PROCESS izhikevich + RANGE a, b, c, d, vthresh, Cm, uinit, t_reset, g_reset + NONSPECIFIC_CURRENT i +} + +UNITS { + (mV) = (millivolt) + (nA) = (nanoamp) + (nF) = (nanofarad) + (uS) = (microsiemens) + (ms) = (millisecond) +} + +PARAMETER { + a = 0.02 (/ms) + b = 0.2 (/ms) + c = -65 (mV) : reset potential + d = 2 (mV/ms) : reset increment of u + vthresh = 30 (mV) : spike / reset threshold + Cm = 0.001 (nF) : capacitance (matches the cell's absolute cm) + uinit = -14 (mV/ms) + t_reset = 0.001 (ms) : clamp duration (<= dt: a one-step reset) + g_reset = 1000 (uS) +} + +STATE { + u (mV/ms) + refrac (ms) +} + +ASSIGNED { clamp } + +INITIAL { + u = uinit + refrac = 0 + clamp = 0 +} + +BREAKPOINT { + SOLVE states METHOD cnexp + clamp = 0 + if (refrac > 0) { clamp = 1 } + i = clamp * g_reset * (v - c) + + (1 - clamp) * (-Cm * (0.04 * v * v + 5 * v + 140 - u)) +} + +DERIVATIVE states { + u' = a * (b * v - u) + refrac' = -1 +} + +POST_EVENT(time) { + refrac = t_reset + u = u + d +} diff --git a/pyNN/arbor/nmodl/lif_gsfa_grr.mod b/pyNN/arbor/nmodl/lif_gsfa_grr.mod new file mode 100644 index 00000000..5f3ffab9 --- /dev/null +++ b/pyNN/arbor/nmodl/lif_gsfa_grr.mod @@ -0,0 +1,73 @@ +: Integrate-and-fire reset with conductance-based spike-frequency adaptation (g_s) +: and a conductance-based relative-refractory mechanism (g_r), for building +: IF_cond_exp_gsfa_grr as an Arbor cable cell (Muller et al. 2007). Combines the +: lif.mod reset (refractory countdown + clamp, driven by the threshold_detector via +: POST_EVENT) with two spike-triggered conductances that are incremented by q_s/q_r +: on each spike and decay with time constants tau_s/tau_r, pulling the membrane +: towards E_s/E_r. The sub-threshold leak is a separate `pas` mechanism, as for lif. +: +: The adaptation current is 0.001 * (g_s*(v-E_s) + g_r*(v-E_r)): with g in nS and v +: in mV, g*(v-E) is in pA, and the 0.001 factor converts it to nA (matching the +: NEURON backend's gsfa_grr.mod). The refractory clamp is gated by the +: voltage-independent 0/1 multiplier `clamp` so modcc derives the conductance +: correctly (see lif.mod). +NEURON { + POINT_PROCESS lif_gsfa_grr + RANGE v_reset, t_ref, g_reset, E_s, E_r, tau_s, tau_r, q_s, q_r + NONSPECIFIC_CURRENT i +} + +UNITS { + (mV) = (millivolt) + (nA) = (nanoamp) + (uS) = (microsiemens) + (nS) = (nanosiemens) + (ms) = (millisecond) +} + +PARAMETER { + v_reset = -65 (mV) + t_ref = 0.1 (ms) + g_reset = 1000 (uS) + E_s = -75 (mV) + E_r = -75 (mV) + tau_s = 100 (ms) + tau_r = 2 (ms) + q_s = 15 (nS) + q_r = 3000 (nS) +} + +STATE { + refrac (ms) + g_s (nS) + g_r (nS) +} + +ASSIGNED { clamp } + +INITIAL { + refrac = 0 + g_s = 0 + g_r = 0 + clamp = 0 +} + +BREAKPOINT { + SOLVE states METHOD cnexp + clamp = 0 + if (refrac > 0) { clamp = 1 } + i = clamp * g_reset * (v - v_reset) + + (0.001) * (g_s * (v - E_s) + g_r * (v - E_r)) +} + +DERIVATIVE states { + refrac' = -1 + g_s' = -g_s/tau_s + g_r' = -g_r/tau_r +} + +POST_EVENT(time) { + refrac = t_ref + g_s = g_s + q_s + g_r = g_r + q_r +} diff --git a/pyNN/arbor/standardmodels.py b/pyNN/arbor/standardmodels.py index 31d02f5a..5bb8b761 100644 --- a/pyNN/arbor/standardmodels.py +++ b/pyNN/arbor/standardmodels.py @@ -16,7 +16,8 @@ from ..parameters import ParameterSpace, IonicSpecies, Sequence from ..morphology import Morphology, NeuriteDistribution, LocationGenerator from .cells import (CellDescriptionBuilder, PointCellDescriptionBuilder, - LIFDynamics, AdExpDynamics) + LIFDynamics, AdExpDynamics, GsfaGrrDynamics, IzhikevichDynamics, + HHDynamics) from .simulator import state from .morphology import LabelledLocations @@ -713,3 +714,88 @@ def translate(self, parameters, copy=True): "excitatory": ("alphasyn", {"tau": "tau_syn_E", "e": "e_rev_E"}), "inhibitory": ("alphasyn", {"tau": "tau_syn_I", "e": "e_rev_I"}), }, parameters.shape, dynamics_class=AdExpDynamics) + + +class IF_cond_exp_gsfa_grr(cells.IF_cond_exp_gsfa_grr): + __doc__ = cells.IF_cond_exp_gsfa_grr.__doc__ + + translations = build_translations( + ('v_rest', 'E_L'), + ('v_reset', 'E_R'), + ('v_thresh', 'V_th'), + ('tau_refrac', 't_ref'), + ('tau_m', 'tau_m'), + ('cm', 'C_m'), + ('i_offset', 'i_offset'), + ('tau_syn_E', 'tau_syn_E'), + ('tau_syn_I', 'tau_syn_I'), + ('e_rev_E', 'e_rev_E'), + ('e_rev_I', 'e_rev_I'), + ('tau_sfa', 'tau_s'), + ('e_rev_sfa', 'E_s'), + ('q_sfa', 'q_s'), + ('tau_rr', 'tau_r'), + ('e_rev_rr', 'E_r'), + ('q_rr', 'q_r'), + ) + arbor_cell_kind = arbor.cell_kind.cable + + def translate(self, parameters, copy=True): + native = super().translate(parameters, copy) + return _point_cell_description(native, { + "excitatory": ("expsyn", {"tau": "tau_syn_E", "e": "e_rev_E"}), + "inhibitory": ("expsyn", {"tau": "tau_syn_I", "e": "e_rev_I"}), + }, parameters.shape, dynamics_class=GsfaGrrDynamics) + + +class Izhikevich(cells.Izhikevich): + __doc__ = cells.Izhikevich.__doc__ + + translations = build_translations( + ('a', 'a'), + ('b', 'b'), + ('c', 'c'), + ('d', 'd'), + ('i_offset', 'i_offset'), + ) + arbor_cell_kind = arbor.cell_kind.cable + # Izhikevich uses voltage-step (delta) synapses, which an Arbor cable cell + # cannot express (a mechanism cannot write v); synaptic input is not yet + # supported, so no receptors are placed. Current injection (i_offset, DCSource) + # drives the cell. + receptor_types = () + + def translate(self, parameters, copy=True): + native = super().translate(parameters, copy) + return _point_cell_description( + native, {}, parameters.shape, dynamics_class=IzhikevichDynamics) + + +class HH_cond_exp(cells.HH_cond_exp): + __doc__ = cells.HH_cond_exp.__doc__ + + translations = build_translations( + ('gbar_Na', 'gnabar'), + ('gbar_K', 'gkbar'), + ('g_leak', 'gl'), + ('e_rev_Na', 'ena'), + ('e_rev_K', 'ek'), + ('e_rev_leak', 'el'), + ('v_offset', 'vT'), + ('cm', 'C_m'), + ('i_offset', 'i_offset'), + ('tau_syn_E', 'tau_syn_E'), + ('tau_syn_I', 'tau_syn_I'), + ('e_rev_E', 'e_rev_E'), + ('e_rev_I', 'e_rev_I'), + ) + arbor_cell_kind = arbor.cell_kind.cable + # the standard model also lists a gap-junction receptor, not supported here + receptor_types = ('excitatory', 'inhibitory') + + def translate(self, parameters, copy=True): + native = super().translate(parameters, copy) + return _point_cell_description(native, { + "excitatory": ("expsyn", {"tau": "tau_syn_E", "e": "e_rev_E"}), + "inhibitory": ("expsyn", {"tau": "tau_syn_I", "e": "e_rev_I"}), + }, parameters.shape, dynamics_class=HHDynamics) diff --git a/test/system/scenarios/test_cell_types.py b/test/system/scenarios/test_cell_types.py index fab2585c..0ad6b70c 100644 --- a/test/system/scenarios/test_cell_types.py +++ b/test/system/scenarios/test_cell_types.py @@ -38,7 +38,7 @@ def test_EIF_cond_alpha_isfa_ista(sim, plot_figure=False): return data -@run_with_simulators("nest", "neuron", "brian2") +@run_with_simulators("arbor", "nest", "neuron", "brian2") def test_HH_cond_exp(sim, plot_figure=False): sim.setup(timestep=0.001, min_delay=0.1) cellparams = { @@ -458,6 +458,51 @@ def test_EIF_cond_exp_isfa_ista_adaptation(sim): assert isis[-1] > 1.5 * isis[0], isis +@run_with_simulators("arbor", "neuron") +def test_Izhikevich_regular_spiking(sim): + """A constant current into a regular-spiking Izhikevich neuron produces a + steady spike train. Arbor's quadratic dynamics and one-step reset are + cross-checked against NEURON to <0.03 ms on spike times / <0.1% on ISI.""" + sim.setup(timestep=0.01, min_delay=0.1) + cell = sim.Population(1, sim.Izhikevich(a=0.02, b=0.2, c=-65.0, d=8.0, i_offset=0.01)) + cell.initialize(v=-70.0, u=-14.0) + cell.record('spikes') + sim.run(300.0) + st = cell.get_data().segments[0].spiketrains[0].magnitude + sim.end() + assert 6 <= len(st) <= 10, len(st) + isis = np.diff(st) + # after the first ISI the train is regular (steady inter-spike interval) + assert abs(isis[-1] - isis[-2]) < 0.1, isis + assert 40.0 < isis[-1] < 50.0, isis + + +@run_with_simulators("arbor", "neuron") +def test_IF_cond_exp_gsfa_grr_adaptation(sim): + """A constant supra-threshold current into IF_cond_exp_gsfa_grr produces spike- + frequency adaptation from the g_s conductance: the ISIs lengthen and then settle + to a steady value. Arbor is cross-checked against NEURON to <0.01 ms on spike + times; nest/brian2 are excluded here because their gsfa_grr implementations + adapt differently for this strongly-driven regime.""" + sim.setup(timestep=0.01, min_delay=0.1) + cell = sim.Population(1, sim.IF_cond_exp_gsfa_grr( + v_rest=-65.0, v_reset=-65.0, v_thresh=-57.0, tau_m=10.0, cm=0.25, + tau_refrac=2.0, i_offset=1.0, + tau_sfa=100.0, q_sfa=15.0, e_rev_sfa=-75.0, + tau_rr=2.0, q_rr=3000.0, e_rev_rr=-75.0)) + cell.initialize(v=-65.0) + cell.record('spikes') + sim.run(400.0) + st = cell.get_data().segments[0].spiketrains[0].magnitude + sim.end() + assert len(st) >= 8, len(st) + isis = np.diff(st) + # early ISIs lengthen (adaptation) and later ones settle to a steady rate + assert isis[3] > isis[0], isis + assert isis[-1] > 2.0 * isis[0], isis # substantial adaptation + assert abs(isis[-1] - isis[-2]) < 0.1, isis # steady state reached + + @run_with_simulators("arbor", "neuron", "nest", "brian2") def test_IF_curr_alpha_EPSP(sim): """A single presynaptic spike into IF_curr_alpha produces an alpha-shaped EPSP diff --git a/test/unittests/test_arbor.py b/test/unittests/test_arbor.py index 2be63188..e2c21c9f 100644 --- a/test/unittests/test_arbor.py +++ b/test/unittests/test_arbor.py @@ -273,6 +273,49 @@ def test_eif_cond_exp_translation(self): def test_adexp_component_uses_adexp_dynamics(self): self.assertIs(arbor_standardmodels.AdExp.dynamics_class, arbor_cells.AdExpDynamics) + def test_hh_cond_exp_translation(self): + ct = arbor_standardmodels.HH_cond_exp(gbar_Na=20.0, e_rev_Na=50.0, tau_syn_E=0.2) + self.assertIs(ct.arbor_cell_kind, arbor.cell_kind.cable) + self.assertTrue(ct.conductance_based) + self.assertEqual(tuple(ct.receptor_types), ("excitatory", "inhibitory")) + builder = ct.native_parameters["cell_description"].base_value + builder.set_shape((1,)) + self.assertIsInstance(builder.dynamics, arbor_cells.HHDynamics) + p = builder.dynamics.neuron_parameters + self.assertAlmostEqual(p["gnabar"][0], 20.0) # whole-cell uS; -> S/cm2 at build + self.assertAlmostEqual(p["ena"][0], 50.0) + exc_model, exc_params = builder.post_synaptic_receptors["excitatory"] + self.assertEqual(exc_model, "expsyn") + + def test_hh_specific_conductance_conversion(self): + # 20 uS whole-cell over the 1e-3 cm2 point-cell area -> 0.02 S/cm2 + dyn = arbor_cells.HHDynamics(None) + self.assertAlmostEqual(dyn._specific_g(20.0), 0.02) + + def test_izhikevich_translation(self): + ct = arbor_standardmodels.Izhikevich(a=0.02, b=0.2, c=-65.0, d=8.0) + self.assertIs(ct.arbor_cell_kind, arbor.cell_kind.cable) + self.assertEqual(tuple(ct.receptor_types), ()) # voltage-step synapses unsupported + builder = ct.native_parameters["cell_description"].base_value + builder.set_shape((1,)) + self.assertIsInstance(builder.dynamics, arbor_cells.IzhikevichDynamics) + p = builder.dynamics.neuron_parameters + self.assertAlmostEqual(p["a"][0], 0.02) + self.assertAlmostEqual(p["d"][0], 8.0) + self.assertEqual(builder.post_synaptic_receptors, {}) + + def test_if_cond_exp_gsfa_grr_translation(self): + ct = arbor_standardmodels.IF_cond_exp_gsfa_grr( + tau_sfa=120.0, q_sfa=12.0, e_rev_sfa=-70.0, tau_rr=3.0) + builder = ct.native_parameters["cell_description"].base_value + builder.set_shape((1,)) + self.assertIsInstance(builder.dynamics, arbor_cells.GsfaGrrDynamics) + p = builder.dynamics.neuron_parameters + self.assertAlmostEqual(p["tau_s"][0], 120.0) + self.assertAlmostEqual(p["q_s"][0], 12.0) + self.assertAlmostEqual(p["E_s"][0], -70.0) + self.assertAlmostEqual(p["tau_r"][0], 3.0) + def test_reset_and_current_synapse_mechanisms_in_catalogue(self): cat = arbor.load_catalogue(arbor_simulator.catalogue_path()) mechs = list(cat) @@ -281,6 +324,9 @@ def test_reset_and_current_synapse_mechanisms_in_catalogue(self): self.assertIn("alphasyn", mechs) self.assertIn("alphasyn_curr", mechs) self.assertIn("adexp", mechs) + self.assertIn("lif_gsfa_grr", mechs) + self.assertIn("izhikevich", mechs) + self.assertIn("hh_traub", mechs) @unittest.skipUnless(have_arbor, "Requires Arbor")