diff --git a/pyNN/arbor/cells.py b/pyNN/arbor/cells.py index 6cbbfb33..0e186783 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 @@ -10,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) @@ -33,9 +26,67 @@ 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, 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({ + "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): 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 +101,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 @@ -134,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") @@ -152,18 +191,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 +210,330 @@ 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 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 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 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. + + 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. + + ``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_* standard models (from their translations). + """ + + def __init__(self, dynamics, post_synaptic_receptors): + super().__init__() + self.dynamics = dynamics + self.post_synaptic_receptors = post_synaptic_receptors + self.shape = None + + def set_shape(self, value): + self.shape = value + self.dynamics.set_shape(value) + for (_model, synapse_parameters) in self.post_synaptic_receptors.values(): + synapse_parameters.shape = value + + 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): + labels = { + "all": "(all)", "soma": "(tag 1)", "root": "(root)", + } + decor = arbor.decor() + # 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(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 + 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 = self.dynamics.i_offset(i) + if i_offset != 0.0: + self._place_iclamp(decor, "(root)", + tstart=0.0, duration=1e12, current=i_offset, + label="i_offset") + # 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) + + 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/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/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/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.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/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/populations.py b/pyNN/arbor/populations.py index 02237271..0936a469 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,12 +168,23 @@ 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 - # 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/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..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 @@ -90,7 +92,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..5bb8b761 100644 --- a/pyNN/arbor/standardmodels.py +++ b/pyNN/arbor/standardmodels.py @@ -13,9 +13,11 @@ 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 +from .cells import (CellDescriptionBuilder, PointCellDescriptionBuilder, + LIFDynamics, AdExpDynamics, GsfaGrrDynamics, IzhikevichDynamics, + HHDynamics) from .simulator import state from .morphology import LabelledLocations @@ -48,52 +50,152 @@ class SpikeSourceArray(cells.SpikeSourceArray): arbor_schedule_units = {"times": U.ms} -class BaseCurrentSource(object): - pass - - -class DCSource(BaseCurrentSource, electrodes.DCSource): - __doc__ = electrodes.DCSource.__doc__ +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( - ('amplitude', 'current'), - ('start', 'tstart'), - ('stop', 'duration', "stop - start", "tstart + duration") + ('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, + } + + +# --- 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 + + +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` ? if hasattr(cells, "parent"): - cell_descr = cells.parent._arbor_cell_description.base_value - index = cells.parent.id_to_index(cells.all_cells.astype(int)) + 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 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: - 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): - # 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__ @@ -103,6 +205,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__ @@ -116,8 +225,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'), @@ -127,6 +249,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__ @@ -320,3 +456,346 @@ 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 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__ + + 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"} + # the PointNeuronDynamics that realises this neuron component as a cable cell + 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. + + 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. + """ + 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(dynamics, 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 + + +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. + + ``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 dynamics_class.native_names}, shape=shape) + dynamics = dynamics_class(neuron_parameters) + 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(dynamics, 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) + + +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) + + +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__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 96e3cb8e..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 = { @@ -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,251 @@ 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() + + +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", "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") +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 + 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 + 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..d290d408 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) @@ -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/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 5bc81e05..e2c21c9f 100644 --- a/test/unittests/test_arbor.py +++ b/test/unittests/test_arbor.py @@ -87,11 +87,14 @@ 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) + 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") @@ -125,6 +128,283 @@ 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 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.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_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) + self.assertIn("lif", mechs) + self.assertIn("expsyn_curr", mechs) + 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") +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):