Add support for NESTML code generation for NEST Simulator

.github/workflows/full-test.yml

@@ -56,7 +56,7 @@ jobs:
         python -m pip install arbor==0.9.0 libNeuroML
     - name: Install PyNN itself
       run: |
-        python -m pip install -e ".[test]"
+        python -m pip install -e ".[test,nestml]"
     - name: Test installation has worked (Ubuntu)
       # this is needed because the PyNN tests are just skipped if the simulator
       # fails to install, so we need to catch import failures separately

.gitignore

@@ -93,3 +93,18 @@ tmp
 doc/logos
 docker_*.eggs/
 /test/system/tmp_serialization_test
+
+# Test coverage output
+coverage.lcov
+coverage.xml
+
+# Compiled shared libraries (NEURON, Arbor, NESTML)
+*.so
+pyNN/neuron/nmodl/*/
+
+# Simulation output and NESTML compilation targets
+**/Results/
+**/target/
+
+# Scratch directories
+tmp_*/

Makefile

@@ -99,7 +99,7 @@ $(NEST_STAMP): $(NEST_VENV_STAMP) $(NEST_SRC_UNPACKED)
 	$(NEST_PIP) install \
 	    "neuron>=9.0.0" nrnutils "arbor==0.9.0" \
 	    brian2 libNeuroML scipy matplotlib Cheetah3 h5py Jinja2 \
-	    pytest pytest-xdist pytest-cov flake8
+	    pytest pytest-xdist pytest-cov flake8 morphio nestml
 	$(NEST_PIP) install -e .
 	# Compile NEURON .mod mechanisms against the venv's NEURON.
 	# The compiled arm64/ dir lives in the source tree and is version-specific,

doc/backends/NEST.txt

@@ -165,4 +165,137 @@ implications for their usage in PyNN:
   However, they will only deviate from the default when changed manually.
 
 
+Using models defined in NESTML
+==============================
+
+`NESTML`_ is a domain-specific language for defining neuron and synapse models that can be
+compiled to efficient C++ code for use in NEST. PyNN wraps PyNESTML to compile and install
+NESTML models automatically at ``sim.setup()`` time.
+
+Installation
+------------
+
+NESTML support requires the ``nestml`` package (imported as ``pynestml``), which is not
+installed by default:
+
+.. code-block:: bash
+
+    pip install PyNN[nestml]   # installs PyNN together with nestml
+    # or
+    pip install nestml         # installs nestml into an existing PyNN environment
+
+Important: register models before ``sim.setup()``
+--------------------------------------------------
+
+All NESTML models registered for a simulation are compiled together in a single PyNESTML
+pass when ``sim.setup()`` is called. Both :func:`nestml_cell_type` and
+:func:`nestml_synapse_type` must therefore be called *before* ``sim.setup()``. Calling
+them afterwards raises a ``RuntimeError``.
+
+The required call order is:
+
+.. code-block:: python
+
+    import pyNN.nest as sim
+    from pyNN.nest import nestml
+
+    # 1. Register NESTML models — no compilation yet
+    MyNeuron = nestml.nestml_cell_type("my_neuron", "my_neuron.nestml")
+
+    # 2. setup() triggers the single compile + install pass
+    sim.setup(timestep=0.1, min_delay=1.0)
+
+    # 3. MyNeuron now behaves identically to native_cell_type() results
+    pop = sim.Population(100, MyNeuron(param=value))
+
+
+NESTML neuron models
+--------------------
+
+Use :func:`~pyNN.nest.nestml.nestml_cell_type` with the model name and either a path to a
+``.nestml`` file or a string containing NESTML source code inline:
+
+.. code-block:: python
+
+    # From a file
+    MyNeuron = nestml.nestml_cell_type("my_neuron", "/path/to/my_neuron.nestml")
+
+    # Inline NESTML source
+    nestml_source = """
+    model my_neuron:
+        ...
+    """
+    MyNeuron = nestml.nestml_cell_type("my_neuron", nestml_source)
+
+After ``sim.setup()`` the returned class behaves identically to one returned by
+:func:`native_cell_type`: parameters can be set, state variables initialised, and
+variables recorded in the usual way.
+
+See ``examples/nestml/wang_buzsaki_synaptic_input.py`` for a file-based example and
+``examples/nestml/wang_buzsaki_current_injection.py`` for the inline variant.
+
+
+NESTML synapse models
+---------------------
+
+Use :func:`~pyNN.nest.nestml.nestml_synapse_type` to register a synapse model. The
+``weight_variable`` and ``delay_variable`` arguments identify which variables in the NESTML
+model correspond to the connection weight and dendritic delay respectively:
+
+.. code-block:: python
+
+    TsodyksSyn = nestml.nestml_synapse_type(
+        "tsodyks_synapse_nestml",
+        "/path/to/tsodyks_synapse.nestml",
+        weight_variable="w",
+        delay_variable="d",
+    )
+    sim.setup(timestep=0.1, min_delay=1.0)
+
+    prj = sim.Projection(
+        source, target,
+        sim.AllToAllConnector(),
+        TsodyksSyn(weight=1.0, delay=1.0),
+        receptor_type="excitatory",
+    )
+
+
+Plastic synapses requiring co-generation
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Some NESTML plasticity models (such as STDP) must be co-generated with a specific
+postsynaptic neuron model so that PyNESTML can wire up the pre/post spike communication.
+Pass the neuron description via ``postsynaptic_neuron_nestml_description``; the
+co-generated neuron class is then accessible as an attribute of the synapse class before
+and after ``sim.setup()``:
+
+.. code-block:: python
+
+    stdp = nestml.nestml_synapse_type(
+        "stdp_synapse",
+        "stdp_synapse.nestml",
+        postsynaptic_neuron_nestml_description="iaf_psc_exp_neuron.nestml",
+    )
+    PostNeuron = stdp.postsynaptic_cell_type  # available immediately, before setup()
+
+    sim.setup(timestep=0.1, min_delay=1.0)
+
+    source = sim.Population(10, sim.SpikeSourcePoisson(rate=50.0))
+    target = sim.Population(10, PostNeuron())
+    prj = sim.Projection(source, target, sim.AllToAllConnector(),
+                         stdp(weight=1.0, delay=1.0), receptor_type="excitatory")
+
+See ``examples/nestml/stdp.py`` for a complete working example.
+
+
+Future backends
+---------------
+
+NESTML support is currently specific to the NEST backend. The SpiNNaker backend
+(developed separately as `sPyNNaker`_) also has partial NESTML support, and it is hoped
+that other backends may gain NESTML support in future.
+
+
+.. _`NESTML`: https://nestml.readthedocs.io/
+.. _`sPyNNaker`: https://github.com/SpiNNakerManchester/sPyNNaker
 .. _`NEST random number generator`: https://nest-simulator.readthedocs.io/en/stable/guides/random_numbers.html#select-the-type-of-random-number-generator

examples/nestml/iaf_psc_exp_neuron.nestml

@@ -0,0 +1,130 @@
+# iaf_psc_exp - Leaky integrate-and-fire neuron model
+# ###################################################
+#
+# Description
+# +++++++++++
+#
+# iaf_psc_exp is an implementation of a leaky integrate-and-fire model
+# with exponentially decaying synaptic currents according to [1]_.
+# Thus, postsynaptic currents have an infinitely short rise time.
+#
+# The threshold crossing is followed by an absolute refractory period
+# during which the membrane potential is clamped to the resting potential
+# and spiking is prohibited.
+#
+# The general framework for the consistent formulation of systems with
+# neuron like dynamics interacting by point events is described in
+# [1]_.  A flow chart can be found in [2]_.
+#
+# Critical tests for the formulation of the neuron model are the
+# comparisons of simulation results for different computation step
+# sizes.
+#
+# .. note::
+#
+#    If tau_m is very close to tau_syn_exc or tau_syn_inh, numerical problems
+#    may arise due to singularities in the propagator matrics. If this is
+#    the case, replace equal-valued parameters by a single parameter.
+#
+# For details, please see ``IAF_neurons_singularity.ipynb`` in
+# the NEST source code (``docs/model_details``).
+#
+#
+# References
+# ++++++++++
+#
+# .. [1] Rotter S,  Diesmann M (1999). Exact simulation of
+#        time-invariant linear systems with applications to neuronal
+#        modeling. Biologial Cybernetics 81:381-402.
+#        DOI: https://doi.org/10.1007/s004220050570
+# .. [2] Diesmann M, Gewaltig M-O, Rotter S, & Aertsen A (2001). State
+#        space analysis of synchronous spiking in cortical neural
+#        networks. Neurocomputing 38-40:565-571.
+#        DOI: https://doi.org/10.1016/S0925-2312(01)00409-X
+# .. [3] Morrison A, Straube S, Plesser H E, Diesmann M (2006). Exact
+#        subthreshold integration with continuous spike times in discrete time
+#        neural network simulations. Neural Computation, in press
+#        DOI: https://doi.org/10.1162/neco.2007.19.1.47
+#
+#
+# See also
+# ++++++++
+#
+# iaf_psc_delta, iaf_psc_alpha, iaf_cond_exp
+#
+#
+# Copyright statement
+# +++++++++++++++++++
+#
+# This file is part of NEST.
+#
+# Copyright (C) 2004 The NEST Initiative
+#
+# NEST is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 2 of the License, or
+# (at your option) any later version.
+#
+# NEST is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with NEST.  If not, see <http://www.gnu.org/licenses/>.
+#
+#
+model iaf_psc_exp_neuron:
+
+    state:
+        V_m mV = E_L    # Membrane potential
+        refr_t ms = 0 ms    # Refractory period timer
+        I_syn_exc pA = 0 pA
+        I_syn_inh pA = 0 pA
+
+    equations:
+        I_syn_exc' = -I_syn_exc / tau_syn_exc
+        I_syn_inh' = -I_syn_inh / tau_syn_inh
+        V_m' = -(V_m - E_L) / tau_m + (I_syn_exc - I_syn_inh + I_e + I_stim) / C_m
+        refr_t' = -1e3 * ms/s    # refractoriness is implemented as an ODE, representing a timer counting back down to zero. XXX: TODO: This should simply read ``refr_t' = -1 / s`` (see https://github.com/nest/nestml/issues/984)
+
+    parameters:
+        C_m pF = 250 pF           # Capacitance of the membrane
+        tau_m ms = 10 ms          # Membrane time constant
+        tau_syn_inh ms = 2 ms     # Time constant of inhibitory synaptic current
+        tau_syn_exc ms = 2 ms     # Time constant of excitatory synaptic current
+        refr_T ms = 2 ms          # Duration of refractory period
+        E_L mV = -70 mV           # Resting potential
+        V_reset mV = -70 mV       # Reset value of the membrane potential
+        V_th mV = -55 mV          # Spike threshold potential
+
+        # constant external input current
+        I_e pA = 0 pA
+
+    input:
+        exc_spikes <- excitatory spike
+        inh_spikes <- inhibitory spike
+        I_stim pA <- continuous
+
+    output:
+        spike
+
+    update:
+        if refr_t > 0 ms:
+            # neuron is absolute refractory, do not evolve V_m
+            integrate_odes(I_syn_exc, I_syn_inh, refr_t)
+        else:
+            # neuron not refractory
+            integrate_odes(I_syn_exc, I_syn_inh, V_m)
+
+    onReceive(exc_spikes):
+        I_syn_exc += exc_spikes * pA * s
+
+    onReceive(inh_spikes):
+        I_syn_inh += inh_spikes * pA * s
+
+    onCondition(refr_t <= 0 ms and V_m >= V_th):
+        # threshold crossing
+        refr_t = refr_T    # start of the refractory period
+        V_m = V_reset
+        emit_spike()