diff --git a/.github/workflows/full-test.yml b/.github/workflows/full-test.yml
index ca2c2507..be0a6196 100644
--- a/.github/workflows/full-test.yml
+++ b/.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
diff --git a/.gitignore b/.gitignore
index 4702d448..92cfe240 100644
--- a/.gitignore
+++ b/.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_*/
\ No newline at end of file
diff --git a/Makefile b/Makefile
index ccc281db..a98014c8 100644
--- a/Makefile
+++ b/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,
diff --git a/doc/backends/NEST.txt b/doc/backends/NEST.txt
index 13ef1624..3995f7df 100644
--- a/doc/backends/NEST.txt
+++ b/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
\ No newline at end of file
diff --git a/examples/nestml/iaf_psc_exp_neuron.nestml b/examples/nestml/iaf_psc_exp_neuron.nestml
new file mode 100644
index 00000000..1134f20b
--- /dev/null
+++ b/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()
diff --git a/examples/nestml/stdp.py b/examples/nestml/stdp.py
new file mode 100644
index 00000000..99f918a7
--- /dev/null
+++ b/examples/nestml/stdp.py
@@ -0,0 +1,117 @@
+"""
+Example of using a cell type defined in NESTML.
+
+This example requires that PyNN be installed using the "NESTML" option, i.e.
+
+  $ pip install PyNN[NESTML]
+
+or you can install NESTML directly:
+
+  $ pip install nestml
+
+
+Usage: python wang_buzsaki_synaptic_input.py [-h] [--plot-figure] simulator
+
+positional arguments:
+  simulator      nest or spinnaker
+
+optional arguments:
+  -h, --help     show this help message and exit
+  --plot-figure  Plot the simulation results to a file
+  --debug        Print debugging information
+
+"""
+
+import os
+from pyNN.utility import get_simulator, init_logging, normalized_filename, SimulationProgressBar
+from pyNN.random import NumpyRNG, RandomDistribution
+
+
+# === Configure the simulator ================================================
+
+sim, options = get_simulator(
+    (
+        "--plot-figure",
+        "Plot the simulation results to a file.",
+        {"action": "store_true"},
+    ),
+    ("--debug", "Print debugging information"),
+)
+
+if options.debug:
+    init_logging(None, debug=True)
+
+# === Register the NESTML synapse (and co-generated neuron) before sim.setup()
+
+current_dir = os.path.dirname(os.path.abspath(__file__))
+synapse_type = sim.nestml.nestml_synapse_type(
+    "stdp_synapse",
+    os.path.join(current_dir, "stdp_synapse.nestml"),
+    postsynaptic_neuron_nestml_description=os.path.join(current_dir, "iaf_psc_exp_neuron.nestml"),
+)
+post_cell_type = synapse_type.postsynaptic_cell_type
+
+sim.setup(timestep=0.01, min_delay=1.0)
+
+
+# === Build and instrument the network =======================================
+
+p2 = sim.Population(10, sim.SpikeSourcePoisson())
+p1 = sim.Population(10, post_cell_type())
+
+connector = sim.AllToAllConnector()
+
+prj = sim.Projection(p2, p1, connector, receptor_type='excitatory', synapse_type=synapse_type())
+
+p1.record(["V_m", "I_syn_exc", "I_syn_inh"])
+
+# === Run the simulation =====================================================
+
+print("Running simulation")
+t_stop = 100.0
+pb = SimulationProgressBar(t_stop / 10, t_stop)
+
+sim.run(t_stop, callbacks=[pb])
+
+
+# === Save the results, optionally plot a figure =============================
+
+print("Saving results")
+filename = normalized_filename("Results", "nestml_example", "pkl", options.simulator)
+p1.write_data(filename, annotations={"script_name": __file__})
+
+if options.plot_figure:
+    from pyNN.utility.plotting import Figure, Panel
+
+    figure_filename = filename.replace("pkl", "png")
+    Figure(
+        Panel(
+            p1.get_data().segments[0].filter(name="V_m")[0],
+            ylabel="Membrane potential (mV)",
+            data_labels=[p1.label],
+            yticks=True,
+        ),
+        Panel(
+            p1.get_data().segments[0].filter(name="I_syn_exc")[0],
+            ylabel="Excitatory synaptic current (pA)",
+            data_labels=[p1.label],
+            yticks=True,
+            #ylim=(0, 1),
+        ),
+        Panel(
+            p1.get_data().segments[0].filter(name="I_syn_inh")[0],
+            xticks=True,
+            xlabel="Time (ms)",
+            ylabel="Inhibitory synaptic current (pA)",
+            data_labels=[p1.label],
+            yticks=True,
+            #ylim=(0, 1),
+        ),
+        title="Responses of Wang-Buzsaki neuron model, defined in NESTML, to synaptic input",
+        annotations="Simulated with %s" % options.simulator.upper(),
+    ).save(figure_filename)
+    print(figure_filename)
+
+# === Clean up and quit ========================================================
+
+sim.end()
diff --git a/examples/nestml/stdp_synapse.nestml b/examples/nestml/stdp_synapse.nestml
new file mode 100644
index 00000000..608ceb13
--- /dev/null
+++ b/examples/nestml/stdp_synapse.nestml
@@ -0,0 +1,108 @@
+# stdp_synapse - Synapse model for spike-timing dependent plasticity
+# ##################################################################
+#
+#
+# Description
+# +++++++++++
+#
+# stdp_synapse is a synapse with spike-timing dependent plasticity (as defined in [1]_). Here the weight dependence exponent can be set separately for potentiation and depression. Examples:
+#
+# =================== ==== =============================
+# Multiplicative STDP [2]_ mu_plus = mu_minus = 1
+# Additive STDP       [3]_ mu_plus = mu_minus = 0
+# Guetig STDP         [1]_ mu_plus, mu_minus in [0, 1]
+# Van Rossum STDP     [4]_ mu_plus = 0 mu_minus = 1
+# =================== ==== =============================
+#
+# Note that in this particular STDP synapse model, the weight is not allowed to be negative. In case an inhibitory STDP synapse needs to be modeled, this model (with weight >= 0 at all times) can be connected to a postsynaptic neuron at its appropriate (inhibitory) input port. The sign of the postsynaptic response is thus handled in the postsynaptic neuron. In principle, an STDP synapse model can be defined that allows for negative weights, but in this case, care should be taken to prevent the sign of the weight from changing during learning, as a biological synapse cannot simply switch from one type to another, say, from glutamatergic to GABAergic.
+#
+#
+# References
+# ++++++++++
+#
+# .. [1] Guetig et al. (2003) Learning Input Correlations through Nonlinear
+#        Temporally Asymmetric Hebbian Plasticity. Journal of Neuroscience
+#
+# .. [2] Rubin, J., Lee, D. and Sompolinsky, H. (2001). Equilibrium
+#        properties of temporally asymmetric Hebbian plasticity, PRL
+#        86,364-367
+#
+# .. [3] Song, S., Miller, K. D. and Abbott, L. F. (2000). Competitive
+#        Hebbian learning through spike-timing-dependent synaptic
+#        plasticity,Nature Neuroscience 3:9,919--926
+#
+# .. [4] van Rossum, M. C. W., Bi, G-Q and Turrigiano, G. G. (2000).
+#        Stable Hebbian learning from spike timing-dependent
+#        plasticity, Journal of Neuroscience, 20:23,8812--8821
+#
+#
+# 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 stdp_synapse:
+    state:
+        w real = 1 [[w >= 0]]   # Synaptic weight
+        pre_trace real = 0.
+        post_trace real = 0.
+
+    parameters:
+        d ms = 1 ms    # Synaptic transmission delay
+        lambda real =  0.01         # (dimensionless) learning rate for causal updates
+        alpha real = 1              # relative learning rate for acausal firing
+        tau_tr_pre ms = 20 ms       # time constant of presynaptic trace
+        tau_tr_post ms = 20 ms      # time constant of postsynaptic trace
+        mu_plus real = 1            # weight dependence exponent for causal updates
+        mu_minus real = 1           # weight dependence exponent for acausal updates
+
+        Wmin real = 0. [[Wmin >= 0]]             # minimum absolute value of synaptic weight
+        Wmax real = 100. [[Wmax >= 0]]           # maximum absolute value of synaptic weight
+
+    equations:
+        pre_trace' = -pre_trace / tau_tr_pre
+        post_trace' = -post_trace / tau_tr_post
+
+    input:
+        pre_spikes <- spike
+        post_spikes <- spike
+
+    output:
+        spike(weight real, delay ms)
+
+    onReceive(post_spikes):
+        post_trace += 1
+
+        # potentiate synapse
+        # note that w is at all times >= 0
+        w_ real = Wmax * (w / Wmax + (lambda * (1. - (w / Wmax))**mu_plus * pre_trace))
+        w = min(Wmax, w_)
+
+    onReceive(pre_spikes):
+        pre_trace += 1
+
+        # depress synapse
+        # note that w is at all times >= 0
+        w_ real = Wmax * (w / Wmax - (alpha * lambda * (w / Wmax)**mu_minus * post_trace))
+        w = max(Wmin, w_)
+
+        # deliver spike to postsynaptic partner
+        emit_spike(w, d)
+
+    update:
+        integrate_odes()
diff --git a/examples/nestml/tsodyks_synapse.nestml b/examples/nestml/tsodyks_synapse.nestml
new file mode 100644
index 00000000..e155be3f
--- /dev/null
+++ b/examples/nestml/tsodyks_synapse.nestml
@@ -0,0 +1,42 @@
+model tsodyks_synapse_nestml:
+    parameters:
+        w real = 1          # Synaptic weight
+        d ms = 1 ms         # Dendritic delay (required by PyNESTML; NEST applies transmission delay at connection level)
+        tau_psc ms = 3 ms   # Time constant of postsynaptic current
+        tau_fac ms = 0 ms   # Time constant for facilitation (0 = no facilitation)
+        tau_rec ms = 800 ms # Time constant for recovery from depression
+        U real = 0.5        # Utilisation of synaptic efficacy
+
+    state:
+        x real = 1          # Fraction of synaptic resources available
+        y real = 0          # Fraction of resources in use
+        u real = 0.5        # Running value of utilisation
+        t_last_update ms = 0 ms
+
+    input:
+        pre_spikes <- spike
+
+    output:
+        spike(weight real, delay ms)
+
+    onReceive(pre_spikes):
+        dt ms = t - t_last_update
+        t_last_update = t
+
+        Puu real = tau_fac == 0 ms ? 0 : exp(-dt / tau_fac)
+        Pyy real = exp(-dt / tau_psc)
+        Pzz real = exp(-dt / tau_rec)
+        Pxy real = ((Pzz - 1) * tau_rec - (Pyy - 1) * tau_psc) / (tau_psc - tau_rec)
+        Pxz real = 1 - Pzz
+        z real = 1 - x - y
+
+        u *= Puu
+        x += Pxy * y + Pxz * z
+        y *= Pyy
+        u += U * (1 - u)
+
+        delta_y_tsp real = u * x
+        x -= delta_y_tsp
+        y += delta_y_tsp
+
+        emit_spike(delta_y_tsp * w, d)
diff --git a/examples/nestml/wang_buzsaki_current_injection.py b/examples/nestml/wang_buzsaki_current_injection.py
new file mode 100644
index 00000000..cab3384b
--- /dev/null
+++ b/examples/nestml/wang_buzsaki_current_injection.py
@@ -0,0 +1,208 @@
+"""
+Example of using a cell type defined in NESTML.
+
+This example requires that PyNN be installed using the "NESTML" option, i.e.
+
+  $ pip install PyNN[NESTML]
+
+or you can install NESTML directly:
+
+  $ pip install nestml
+
+
+Usage: python wang_buzsaki_current_injection.py [-h] [--plot-figure] simulator
+
+positional arguments:
+  simulator      nest or spinnaker
+
+optional arguments:
+  -h, --help     show this help message and exit
+  --plot-figure  Plot the simulation results to a file
+  --debug        Print debugging information
+
+"""
+
+from pyNN.utility import get_simulator, init_logging, normalized_filename
+from pyNN.random import NumpyRNG, RandomDistribution
+
+
+# === Configure the simulator ================================================
+
+sim, options = get_simulator(
+    (
+        "--plot-figure",
+        "Plot the simulation results to a file.",
+        {"action": "store_true"},
+    ),
+    ("--debug", "Print debugging information"),
+)
+
+if options.debug:
+    init_logging(None, debug=True)
+
+# === Register the NESTML cell type (must happen before sim.setup())
+
+# note that the NESTML definition can also be stored in a separate file
+
+nestml_description = """
+model wb_cond_exp:
+    state:
+        r integer = 0 # number of steps in the current refractory phase
+
+        V_m mV = E_L    # Membrane potential
+
+        Inact_h real = alpha_h_init / ( alpha_h_init + beta_h_init )
+        Act_n real = alpha_n_init / ( alpha_n_init + beta_n_init )
+
+    equations:
+        # synapses: exponential conductance
+        kernel g_inh = exp(-t / tau_syn_inh)
+        kernel g_exc = exp(-t / tau_syn_exc)
+
+        recordable inline I_syn_exc pA = convolve(g_exc, exc_spikes) * nS * ( V_m - E_exc )
+        recordable inline I_syn_inh pA = convolve(g_inh, inh_spikes) * nS * ( V_m - E_inh )
+
+        inline I_Na  pA = g_Na * _subexpr(V_m) * Inact_h * ( V_m - E_Na )
+        inline I_K   pA = g_K * Act_n**4 * ( V_m - E_K )
+        inline I_L   pA = g_L * ( V_m - E_L )
+
+        V_m' =( -( I_Na + I_K + I_L ) + I_e + I_stim + I_syn_exc - I_syn_inh ) / C_m
+        Act_n' = ( alpha_n(V_m) * ( 1 - Act_n ) - beta_n(V_m) * Act_n )  # n-variable
+        Inact_h' = ( alpha_h(V_m) * ( 1 - Inact_h ) - beta_h(V_m) * Inact_h ) # h-variable
+
+    parameters:
+        t_ref ms = 2 ms           # Refractory period
+        g_Na nS = 3500 nS         # Sodium peak conductance
+        g_K nS = 900 nS           # Potassium peak conductance
+        g_L nS = 10 nS            # Leak conductance
+        C_m pF = 100 pF           # Membrane capacitance
+        E_Na mV = 55 mV           # Sodium reversal potential
+        E_K mV = -90 mV           # Potassium reversal potential
+        E_L mV = -65 mV           # Leak reversal potential (aka resting potential)
+        V_Tr mV = -55 mV          # Spike threshold
+        tau_syn_exc ms = 0.2 ms   # Rise time of the excitatory synaptic alpha function
+        tau_syn_inh ms = 10 ms    # Rise time of the inhibitory synaptic alpha function
+        E_exc mV = 0 mV           # Excitatory synaptic reversal potential
+        E_inh mV = -75 mV         # Inhibitory synaptic reversal potential
+
+        # constant external input current
+        I_e pA = 0 pA
+
+    internals:
+        RefractoryCounts integer = steps(t_ref) # refractory time in steps
+
+        alpha_n_init 1/ms = -0.05/(ms*mV) * (E_L + 34.0 mV) / (exp(-0.1 * (E_L + 34.0 mV)) - 1.0)
+        beta_n_init  1/ms = 0.625/ms * exp(-(E_L + 44.0 mV) / 80.0 mV)
+        alpha_h_init 1/ms = 0.35/ms * exp(-(E_L + 58.0 mV) / 20.0 mV)
+        beta_h_init  1/ms = 5.0 / (exp(-0.1 / mV * (E_L + 28.0 mV)) + 1.0) /ms
+
+    input:
+        inh_spikes <- inhibitory spike
+        exc_spikes <- excitatory spike
+        I_stim pA <- continuous
+
+    output:
+        spike
+
+    update:
+        U_old mV = V_m
+        integrate_odes()
+        # sending spikes: crossing 0 mV, pseudo-refractoriness and local maximum...
+        if r > 0: # is refractory?
+            r -= 1
+        elif V_m > V_Tr and U_old > V_m: # threshold && maximum
+            r = RefractoryCounts
+            emit_spike()
+
+    function _subexpr(V_m mV) real:
+        return alpha_m(V_m)**3 / ( alpha_m(V_m) + beta_m(V_m) )**3
+
+    function alpha_m(V_m mV) 1/ms:
+        return 0.1/(ms*mV) * (V_m + 35.0 mV) / (1.0 - exp(-0.1 mV * (V_m + 35.0 mV)))
+
+    function beta_m(V_m mV) 1/ms:
+        return 4.0/(ms) * exp(-(V_m + 60.0 mV) / 18.0 mV)
+
+    function alpha_n(V_m mV) 1/ms:
+        return -0.05/(ms*mV) * (V_m + 34.0 mV) / (exp(-0.1 * (V_m + 34.0 mV)) - 1.0)
+
+    function beta_n(V_m mV) 1/ms:
+        return 0.625/ms * exp(-(V_m + 44.0 mV) / 80.0 mV)
+
+    function alpha_h(V_m mV) 1/ms:
+        return 0.35/ms * exp(-(V_m + 58.0 mV) / 20.0 mV)
+
+    function beta_h(V_m mV) 1/ms:
+        return 5.0 / (exp(-0.1 / mV * (V_m + 28.0 mV)) + 1.0) /ms
+"""
+
+celltype_cls = sim.nestml.nestml_cell_type("wb_cond_exp", nestml_description)
+
+sim.setup(timestep=0.01, min_delay=1.0)
+
+# add some variability between neurons
+rng = NumpyRNG(seed=1309463846)
+rnd = lambda min, max: RandomDistribution("uniform", (min, max), rng=rng)
+
+celltype = celltype_cls(
+    g_Na=rnd(3400, 3600),  # nS
+    g_K=rnd(850, 950),  # nS
+    g_L=rnd(8, 12),  # nS
+    C_m=rnd(80, 120),  # pF
+    V_Tr=rnd(-56, -54),  # mV
+    I_e=100.0,  # pA
+)
+
+
+# === Build and instrument the network =======================================
+
+cells = sim.Population(5, celltype, label=celltype_cls.__name__)
+
+cells.record(["V_m", "Act_n", "Inact_h"])
+
+
+# === Run the simulation =====================================================
+
+sim.run(100.0)
+
+
+# === Save the results, optionally plot a figure =============================
+
+filename = normalized_filename("Results", "nestml_example", "pkl", options.simulator)
+cells.write_data(filename, annotations={"script_name": __file__})
+
+if options.plot_figure:
+    from pyNN.utility.plotting import Figure, Panel
+
+    figure_filename = filename.replace("pkl", "png")
+    Figure(
+        Panel(
+            cells.get_data().segments[0].filter(name="V_m")[0],
+            ylabel="Membrane potential (mV)",
+            data_labels=[cells.label],
+            yticks=True,
+        ),
+        Panel(
+            cells.get_data().segments[0].filter(name="Act_n")[0],
+            ylabel="Activation variable (n)",
+            data_labels=[cells.label],
+            yticks=True,
+            ylim=(0, 1),
+        ),
+        Panel(
+            cells.get_data().segments[0].filter(name="Inact_h")[0],
+            xticks=True,
+            xlabel="Time (ms)",
+            ylabel="Inactivation variable (h)",
+            data_labels=[cells.label],
+            yticks=True,
+            ylim=(0, 1),
+        ),
+        title="Responses of Wang-Buzsaki neuron model, defined in NESTML, to current injection",
+        annotations="Simulated with %s" % options.simulator.upper(),
+    ).save(figure_filename)
+    print(figure_filename)
+
+# === Clean up and quit ========================================================
+
+sim.end()
diff --git a/examples/nestml/wang_buzsaki_synaptic_input.py b/examples/nestml/wang_buzsaki_synaptic_input.py
new file mode 100644
index 00000000..72f155b7
--- /dev/null
+++ b/examples/nestml/wang_buzsaki_synaptic_input.py
@@ -0,0 +1,125 @@
+"""
+Example of using a cell type defined in NESTML.
+
+This example requires that PyNN be installed using the "NESTML" option, i.e.
+
+  $ pip install PyNN[NESTML]
+
+or you can install NESTML directly:
+
+  $ pip install nestml
+
+
+Usage: python wang_buzsaki_synaptic_input.py [-h] [--plot-figure] simulator
+
+positional arguments:
+  simulator      nest or spinnaker
+
+optional arguments:
+  -h, --help     show this help message and exit
+  --plot-figure  Plot the simulation results to a file
+  --debug        Print debugging information
+
+"""
+
+import os
+from pyNN.utility import get_simulator, init_logging, normalized_filename, SimulationProgressBar
+from pyNN.random import NumpyRNG, RandomDistribution
+
+
+# === Configure the simulator ================================================
+
+sim, options = get_simulator(
+    (
+        "--plot-figure",
+        "Plot the simulation results to a file.",
+        {"action": "store_true"},
+    ),
+    ("--debug", "Print debugging information"),
+)
+
+if options.debug:
+    init_logging(None, debug=True)
+
+# === Register the NESTML cell type (must happen before sim.setup())
+
+current_dir = os.path.dirname(os.path.abspath(__file__))
+input_path = os.path.join(current_dir, "wb_cond_exp_neuron.nestml")
+celltype_cls = sim.nestml.nestml_cell_type("wb_cond_exp_neuron", input_path)
+
+sim.setup(timestep=0.01, min_delay=1.0)
+
+# add some variability between neurons
+rng = NumpyRNG(seed=1309463846)
+rnd = lambda min, max: RandomDistribution("uniform", (min, max), rng=rng)
+
+celltype = celltype_cls(
+    g_Na=rnd(3400, 3600),  # nS
+    g_K=rnd(850, 950),  # nS
+    g_L=rnd(8, 12),  # nS
+    C_m=rnd(80, 120),  # pF
+    V_Tr=rnd(-56, -54),  # mV
+    I_e=0,  # pA
+)
+
+
+# === Build and instrument the network =======================================
+
+inputs = sim.Population(5, sim.SpikeSourcePoisson(rate=100.0), label="input spikes")
+cells = sim.Population(5, celltype, label=celltype_cls.__name__)
+cells.record(["V_m", "I_syn_exc", "I_syn_inh"])
+
+print("Connecting cells")
+syn = sim.StaticSynapse(weight=0.01, delay=1.0)
+prj = sim.Projection(inputs, cells, sim.FixedNumberPostConnector(3), syn, receptor_type="excitatory")
+
+# === Run the simulation =====================================================
+
+print("Running simulation")
+t_stop = 100.0
+pb = SimulationProgressBar(t_stop / 10, t_stop)
+
+sim.run(t_stop, callbacks=[pb])
+
+
+# === Save the results, optionally plot a figure =============================
+
+print("Saving results")
+filename = normalized_filename("Results", "nestml_example", "pkl", options.simulator)
+cells.write_data(filename, annotations={"script_name": __file__})
+
+if options.plot_figure:
+    from pyNN.utility.plotting import Figure, Panel
+
+    figure_filename = filename.replace("pkl", "png")
+    Figure(
+        Panel(
+            cells.get_data().segments[0].filter(name="V_m")[0],
+            ylabel="Membrane potential (mV)",
+            data_labels=[cells.label],
+            yticks=True,
+        ),
+        Panel(
+            cells.get_data().segments[0].filter(name="I_syn_exc")[0],
+            ylabel="Excitatory synaptic current (pA)",
+            data_labels=[cells.label],
+            yticks=True,
+            #ylim=(0, 1),
+        ),
+        Panel(
+            cells.get_data().segments[0].filter(name="I_syn_inh")[0],
+            xticks=True,
+            xlabel="Time (ms)",
+            ylabel="Inhibitory synaptic current (pA)",
+            data_labels=[cells.label],
+            yticks=True,
+            #ylim=(0, 1),
+        ),
+        title="Responses of Wang-Buzsaki neuron model, defined in NESTML, to synaptic input",
+        annotations="Simulated with %s" % options.simulator.upper(),
+    ).save(figure_filename)
+    print(figure_filename)
+
+# === Clean up and quit ========================================================
+
+sim.end()
diff --git a/examples/nestml/wb_cond_exp_neuron.nestml b/examples/nestml/wb_cond_exp_neuron.nestml
new file mode 100644
index 00000000..e5c3be61
--- /dev/null
+++ b/examples/nestml/wb_cond_exp_neuron.nestml
@@ -0,0 +1,142 @@
+# wb_cond_exp - Wang-Buzsaki model
+# ################################
+#
+# Description
+# +++++++++++
+#
+# wb_cond_exp is an implementation of a modified Hodkin-Huxley model.
+#
+# (1) Post-synaptic currents: Incoming spike events induce a post-synaptic change
+# of conductance modeled by an exponential function.
+#
+# (2) Spike Detection: Spike detection is done by a combined threshold-and-local-
+# maximum search: if there is a local maximum above a certain threshold of
+# the membrane potential, it is considered a spike.
+#
+# References
+# ++++++++++
+#
+# .. [1] Wang, X.J. and Buzsaki, G., (1996) Gamma oscillation by synaptic
+#        inhibition in a hippocampal interneuronal network model. Journal of
+#        neuroscience, 16(20), pp.6402-6413.
+#
+# See Also
+# ++++++++
+#
+# hh_cond_exp_traub, wb_cond_multisyn
+#
+#
+# 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 wb_cond_exp_neuron:
+    state:
+        V_m mV = E_L    # Membrane potential
+        V_m_old mV = E_L    # Membrane potential at previous timestep for threshold check
+        refr_t ms = 0 ms    # Refractory period timer
+
+        Inact_h real = alpha_h_init / ( alpha_h_init + beta_h_init )
+        Act_n real = alpha_n_init / ( alpha_n_init + beta_n_init )
+
+    equations:
+        # synapses: exponential conductance
+        kernel g_inh = exp(-t / tau_syn_inh)
+        kernel g_exc = exp(-t / tau_syn_exc)
+
+        recordable inline I_syn_exc pA = convolve(g_exc, exc_spikes) * nS * ( V_m - E_exc )
+        recordable inline I_syn_inh pA = convolve(g_inh, inh_spikes) * nS * ( V_m - E_inh )
+
+        inline I_Na  pA = g_Na * _subexpr(V_m) * Inact_h * ( V_m - E_Na )
+        inline I_K   pA = g_K * Act_n**4 * ( V_m - E_K )
+        inline I_L   pA = g_L * ( V_m - E_L )
+
+        V_m' =( -( I_Na + I_K + I_L ) + I_e + I_stim - I_syn_exc - I_syn_inh ) / C_m
+        Act_n' = ( alpha_n(V_m) * ( 1 - Act_n ) - beta_n(V_m) * Act_n )  # n-variable
+        Inact_h' = ( alpha_h(V_m) * ( 1 - Inact_h ) - beta_h(V_m) * Inact_h ) # h-variable
+        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 = 100 pF           # Membrane capacitance
+        g_Na nS = 3500 nS         # Sodium peak conductance
+        g_K nS = 900 nS           # Potassium peak conductance
+        g_L nS = 10 nS            # Leak conductance
+        E_Na mV = 55 mV           # Sodium reversal potential
+        E_K mV = -90 mV           # Potassium reversal potential
+        E_L mV = -65 mV           # Leak reversal potential (aka resting potential)
+        V_Tr mV = -55 mV          # Spike threshold
+        refr_T ms = 2 ms         # Duration of refractory period
+
+        tau_syn_exc ms = 0.2 ms   # Rise time of the excitatory synaptic alpha function
+        tau_syn_inh ms = 10 ms    # Rise time of the inhibitory synaptic alpha function
+        E_exc mV = 0 mV           # Excitatory synaptic reversal potential
+        E_inh mV = -75 mV         # Inhibitory synaptic reversal potential
+
+        # constant external input current
+        I_e pA = 0 pA
+
+    internals:
+        alpha_n_init 1/ms = -0.05/(ms*mV) * (E_L + 34.0 mV) / (exp(-0.1 * (E_L + 34.0 mV)) - 1.0)
+        beta_n_init  1/ms = 0.625/ms * exp(-(E_L + 44.0 mV) / 80.0 mV)
+        alpha_h_init 1/ms = 0.35/ms * exp(-(E_L + 58.0 mV) / 20.0 mV)
+        beta_h_init  1/ms = 5.0 / (exp(-0.1 / mV * (E_L + 28.0 mV)) + 1.0) /ms
+
+    input:
+        exc_spikes <- excitatory spike
+        inh_spikes <- inhibitory spike
+        I_stim pA <- continuous
+
+    output:
+        spike
+
+    update:
+        # Hodgkin-Huxley type model: ODEs are always integrated, regardless of refractory state
+        V_m_old = V_m
+        if refr_t > 0 ms:
+            # neuron is absolute refractory
+            integrate_odes(V_m, Act_n, Inact_h, refr_t)
+        else:
+            # neuron not refractory
+            integrate_odes(V_m, Act_n, Inact_h)
+
+    onCondition(refr_t <= 0 ms and V_m > 0 mV and V_m_old > V_m):
+        # threshold crossing and maximum
+        refr_t = refr_T    # start of the refractory period
+        emit_spike()
+
+    function _subexpr(V_m mV) real:
+        return alpha_m(V_m)**3 / ( alpha_m(V_m) + beta_m(V_m) )**3
+
+    function alpha_m(V_m mV) 1/ms:
+        return 0.1/(ms*mV) * (V_m + 35.0 mV) / (1.0 - exp(-0.1 mV * (V_m + 35.0 mV)))
+
+    function beta_m(V_m mV) 1/ms:
+        return 4.0/(ms) * exp(-(V_m + 60.0 mV) / 18.0 mV)
+
+    function alpha_n(V_m mV) 1/ms:
+        return -0.05/(ms*mV) * (V_m + 34.0 mV) / (exp(-0.1 * (V_m + 34.0 mV)) - 1.0)
+
+    function beta_n(V_m mV) 1/ms:
+        return 0.625/ms * exp(-(V_m + 44.0 mV) / 80.0 mV)
+
+    function alpha_h(V_m mV) 1/ms:
+        return 0.35/ms * exp(-(V_m + 58.0 mV) / 20.0 mV)
+
+    function beta_h(V_m mV) 1/ms:
+        return 5.0 / (exp(-0.1 / mV * (V_m + 28.0 mV)) + 1.0) /ms
diff --git a/pyNN/models.py b/pyNN/models.py
index 477909cc..9b5f06e8 100644
--- a/pyNN/models.py
+++ b/pyNN/models.py
@@ -116,6 +116,9 @@ class BaseSynapseType(BaseModelType):
     # override for synapses that include an active presynaptic components
     has_presynaptic_components = False
 
+    delay_variable = 'delay'
+    weight_variable = 'weight'
+
     def __init__(self, **parameters):
         """
         `parameters` should be a mapping object, e.g. a dict
@@ -124,12 +127,12 @@ def __init__(self, **parameters):
         if parameters:
             all_parameters.update(**parameters)
         try:
-            if all_parameters['delay'] is None:
-                all_parameters['delay'] = self._get_minimum_delay()
-            if all_parameters['weight'] is None:
-                all_parameters['weight'] = 0.
+            if all_parameters[self.delay_variable] is None:
+                all_parameters[self.delay_variable] = self._get_minimum_delay()
+            if all_parameters[self.weight_variable] is None:
+                all_parameters[self.weight_variable] = 0.
         except KeyError as e:
-            if e.args[0] != 'delay':  # ElectricalSynapses don't have delays
+            if e.args[0] != self.delay_variable:  # ElectricalSynapses don't have delays
                 raise e
         self.parameter_space = ParameterSpace(all_parameters,
                                               self.get_schema(),
diff --git a/pyNN/nest/__init__.py b/pyNN/nest/__init__.py
index 4a5822e3..598a16f3 100644
--- a/pyNN/nest/__init__.py
+++ b/pyNN/nest/__init__.py
@@ -44,6 +44,7 @@
     rank,
 )
 from .procedural_api import create, connect, record, record_v, record_gsyn, set  # noqa: F401
+from . import nestml
 
 
 # ==============================================================================
diff --git a/pyNN/nest/control.py b/pyNN/nest/control.py
index 4600eda9..ca64437a 100644
--- a/pyNN/nest/control.py
+++ b/pyNN/nest/control.py
@@ -79,6 +79,10 @@ def setup(timestep=DEFAULT_TIMESTEP, min_delay=DEFAULT_MIN_DELAY,
     # Set min_delay and max_delay
     simulator.state.set_delays(min_delay, max_delay)
     nest.SetDefaults('spike_generator', {'precise_times': True})
+
+    from . import nestml
+    nestml._compiled = False
+    nestml._compile_and_resolve()
     return rank()
 
 
diff --git a/pyNN/nest/nestml.py b/pyNN/nest/nestml.py
new file mode 100644
index 00000000..98421b9f
--- /dev/null
+++ b/pyNN/nest/nestml.py
@@ -0,0 +1,291 @@
+# -*- coding: utf-8 -*-
+"""
+Support for neuron and synapse models defined in NESTML.
+
+NESTML models must be registered before calling sim.setup(). setup() triggers a
+single compilation pass (via PyNESTML) covering all registered models, then installs
+the resulting NEST module. After setup() the returned classes behave identically to
+those returned by native_cell_type() / native_synapse_type().
+
+Functions:
+    nestml_cell_type    - register a NESTML neuron description; return a cell type class
+    nestml_synapse_type - register a NESTML synapse description; return a synapse type class
+
+:copyright: Copyright 2006-2024 by the PyNN team, see AUTHORS.
+:license: CeCILL, see LICENSE for details.
+"""
+
+import logging
+import os.path
+import re
+import shutil
+import tempfile
+import uuid
+
+from pyNN.nest.cells import native_cell_type
+from pyNN.models import BaseCellType, BaseSynapseType
+from pyNN.nest.synapses import NESTSynapseMixin
+import pyNN
+
+logger = logging.getLogger("PyNN")
+
+_MODULE_NAME = "pynnnestmlmodule"
+
+# Module-level registry — survives state.clear() so definitions persist across setup() calls
+_pending = []    # list of entry dicts, one per registered model
+_compiled = False  # True after the first successful compile + install
+
+
+def _check_not_compiled():
+    if _compiled:
+        raise RuntimeError(
+            "NESTML models must be registered before sim.setup() is called. "
+            "Cannot add new NESTML models after compilation."
+        )
+
+
+def _make_pending_cell_type_class(name):
+    """Return a stub cell type class that will be resolved at setup() time."""
+
+    class NESTMLCellType(BaseCellType):
+        _pending = True
+        _nestml_name = name
+        default_parameters = {}
+        default_initial_values = {}
+        nest_name = {"on_grid": name, "off_grid": name}
+        recordable = []
+        injectable = True
+        uses_parrot = False
+
+        def __init__(self, **parameters):
+            self._pending_parameters = parameters
+
+        def get_schema(self):
+            return {n: type(v) for n, v in self.default_parameters.items()}
+
+        @classmethod
+        def _resolve(cls, nest_name=None):
+            real = native_cell_type(nest_name or cls._nestml_name)
+            for attr in ("nest_name", "nest_model", "default_parameters", "default_initial_values",
+                         "recordable", "injectable", "uses_parrot", "receptor_types",
+                         "standard_receptor_type", "conductance_based", "always_local", "units"):
+                if hasattr(real, attr):
+                    setattr(cls, attr, getattr(real, attr))
+            cls._pending = False
+            cls.__init__ = BaseCellType.__init__
+
+    NESTMLCellType.__name__ = name
+    NESTMLCellType.__qualname__ = name
+    return NESTMLCellType
+
+
+def nestml_cell_type(name, nestml_description):
+    """
+    Register a NESTML neuron description and return a cell type class.
+
+    ``nestml_description`` may be a path to a .nestml file or a string containing
+    NESTML source code.
+
+    The returned class is a stub until sim.setup() is called. setup() compiles all
+    registered NESTML models together in a single PyNESTML invocation and installs
+    the resulting NEST module. After setup() the class behaves identically to one
+    returned by native_cell_type().
+    """
+    _check_not_compiled()
+    stub = _make_pending_cell_type_class(name)
+    _pending.append({"type": "neuron", "name": name, "description": nestml_description, "stub": stub})
+    return stub
+
+
+def _get_model_name(filename: str) -> str:
+    """Extract the model name from a NESTML source file."""
+    with open(filename, "r") as fp:
+        nestml_model = fp.read()
+    return re.findall(r"model [^:\s]*:", nestml_model)[0][6:-1]
+
+
+def nestml_synapse_type(
+        synapse_name,
+        nestml_description,
+        postsynaptic_neuron_nestml_description=None,
+        weight_variable="w",
+        delay_variable="d"
+):
+    """
+    Register a NESTML synapse description and return a synapse type class.
+
+    ``nestml_description`` and (optionally) ``postsynaptic_neuron_nestml_description``
+    may each be a path to a .nestml file or a string containing NESTML source code.
+
+    For plastic synapses that require co-generation with a postsynaptic neuron model,
+    provide ``postsynaptic_neuron_nestml_description``. The co-generated neuron class
+    is then accessible as ``synapse_class.postsynaptic_cell_type`` both before and
+    after sim.setup().
+
+    The returned class is a stub until sim.setup() triggers compilation (see
+    nestml_cell_type() for details).
+    """
+    _check_not_compiled()
+
+    # For co-generation, create the neuron stub now so callers can hold a reference to it
+    # before setup() is called. It will be resolved during _compile_and_resolve().
+    neuron_stub = None
+    if postsynaptic_neuron_nestml_description:
+        neuron_stub = _make_pending_cell_type_class(synapse_name + "_postsynaptic_neuron")
+
+    # Capture in local names to avoid same-name shadowing in the class body below.
+    _delay_var = delay_variable
+    _weight_var = weight_variable
+
+    class NESTMLSynapseType(NESTSynapseMixin, BaseSynapseType):
+        _pending = True
+        _nestml_synapse_name = synapse_name
+        default_parameters = {}
+        nest_name = synapse_name
+        delay_variable = _delay_var
+        weight_variable = _weight_var
+
+        def __init__(self, **parameters):
+            self._pending_parameters = parameters
+
+        def get_schema(self):
+            return {n: type(v) for n, v in self.default_parameters.items()}
+
+        @property
+        def native_parameters(self):
+            return self.parameter_space
+
+        def get_native_names(self, *names):
+            return names
+
+        @classmethod
+        def _resolve(cls, resolved_synapse_name, resolved_neuron_name=None):
+            real = pyNN.nest.native_synapse_type(resolved_synapse_name)
+            for attr in ("nest_name", "default_parameters"):
+                setattr(cls, attr, getattr(real, attr))
+            if resolved_neuron_name and hasattr(cls, 'postsynaptic_cell_type'):
+                cls.postsynaptic_cell_type._resolve(resolved_neuron_name)
+            cls._pending = False
+            cls.__init__ = BaseSynapseType.__init__
+
+    NESTMLSynapseType.__name__ = synapse_name
+    NESTMLSynapseType.__qualname__ = synapse_name
+
+    if neuron_stub is not None:
+        NESTMLSynapseType.postsynaptic_cell_type = neuron_stub
+
+    _pending.append({
+        "type": "synapse",
+        "name": synapse_name,
+        "description": nestml_description,
+        "neuron_description": postsynaptic_neuron_nestml_description,
+        "weight_variable": weight_variable,
+        "delay_variable": delay_variable,
+        "stub": NESTMLSynapseType,
+        "neuron_stub": neuron_stub,
+    })
+    return NESTMLSynapseType
+
+
+def _ensure_file(description, tmpdirs):
+    """
+    Return a path to a .nestml file for the given description.
+
+    If ``description`` is already a file path, return it unchanged. If it is inline
+    NESTML source code, write it to a temporary file and return that path (the temp
+    directory is appended to ``tmpdirs`` for later cleanup).
+    """
+    if os.path.isfile(description):
+        return description
+    tmpdir = tempfile.mkdtemp()
+    tmpdirs.append(tmpdir)
+    filepath = os.path.join(tmpdir, "model.nestml")
+    with open(filepath, "w") as fp:
+        fp.write(description)
+    return filepath
+
+
+def _compile_and_resolve():
+    """
+    Compile all pending NESTML models together and resolve their stub classes.
+
+    Called by sim.setup() after the NEST kernel is initialised. If no models have
+    been registered this is a no-op.
+    """
+    global _compiled, _pending
+
+    if not _pending:
+        _compiled = True
+        return
+
+    from pynestml.frontend.pynestml_frontend import generate_nest_target
+    import nest
+
+    tmpdirs = []
+    input_paths = []
+    codegen_opts = {}
+    synapse_pairs = []  # tracks co-generation relationships for name resolution
+
+    for entry in _pending:
+        filepath = _ensure_file(entry["description"], tmpdirs)
+        input_paths.append(filepath)
+
+        if entry["type"] == "synapse":
+            syn_model = _get_model_name(filepath)
+            entry["_syn_model"] = syn_model  # actual NESTML model name, used for resolution
+            codegen_opts.setdefault("synapse_models", []).append(syn_model)
+            codegen_opts.setdefault("weight_variable", {})[syn_model] = entry["weight_variable"]
+            if entry["delay_variable"] is not None:
+                codegen_opts.setdefault("delay_variable", {})[syn_model] = entry["delay_variable"]
+
+            if entry["neuron_description"]:
+                neuron_filepath = _ensure_file(entry["neuron_description"], tmpdirs)
+                input_paths.append(neuron_filepath)
+                neu_model = _get_model_name(neuron_filepath)
+                synapse_pairs.append({
+                    "neuron": neu_model,
+                    "synapse": syn_model,
+                    "post_ports": ["post_spikes"],
+                    "entry": entry,
+                })
+
+    if synapse_pairs:
+        codegen_opts["neuron_synapse_pairs"] = [
+            {"neuron": p["neuron"], "synapse": p["synapse"], "post_ports": p["post_ports"]}
+            for p in synapse_pairs
+        ]
+
+    # Use a unique module name to prevent filename conflicts when multiple compilations
+    # run concurrently (e.g. pytest-xdist parallel workers).
+    module_name = "pynnnestml_" + uuid.uuid4().hex[:8] + "module"
+
+    build_dir = tempfile.mkdtemp()
+    tmpdirs.append(build_dir)
+    try:
+        generate_nest_target(
+            input_path=input_paths,
+            target_path=build_dir,
+            install_path=None,
+            module_name=module_name,
+            codegen_opts=codegen_opts,
+        )
+        nest.Install(module_name)
+    finally:
+        for tmpdir in tmpdirs:
+            shutil.rmtree(tmpdir, ignore_errors=True)
+
+    for entry in _pending:
+        stub = entry["stub"]
+        if entry["type"] == "neuron":
+            stub._resolve()
+        else:
+            pair = next((p for p in synapse_pairs if p["entry"] is entry), None)
+            if pair:
+                resolved_syn = pair["synapse"] + "__with_" + pair["neuron"]
+                resolved_neu = pair["neuron"] + "__with_" + pair["synapse"]
+                stub._resolve(resolved_syn, resolved_neu)
+            else:
+                stub._resolve(entry["_syn_model"])
+
+    _compiled = True
+    _pending.clear()
diff --git a/pyNN/nest/projections.py b/pyNN/nest/projections.py
index 2b69b7f4..8934d63c 100644
--- a/pyNN/nest/projections.py
+++ b/pyNN/nest/projections.py
@@ -250,17 +250,31 @@ def _convergent_connect(self, presynaptic_indices, postsynaptic_index,
             try:
                 weights = connection_parameter_group.pop('weight')
                 delays = connection_parameter_group.pop('delay')
+
+                weight_var = getattr(self.synapse_type, 'weight_variable', None) or 'weight'
+                delay_var = getattr(self.synapse_type, 'delay_variable', None) or 'delay'
+
+                if weight_var != 'weight':
+                    # Route user weight into the custom weight variable; leave NEST internal weight alone
+                    connection_parameter_group[weight_var] = weights
+                    weights = None
+                if delay_var != 'delay':
+                    # Route user delay into the custom delay variable; leave NEST internal delay alone
+                    connection_parameter_group[delay_var] = delays
+                    delays = None
+
                 # nest.Connect expects a 2D array
-                if not np.isscalar(weights):
+                if weights is not None and not np.isscalar(weights):
                     weights = np.array([weights])
-                if not np.isscalar(delays):
+                if delays is not None and not np.isscalar(delays):
                     delays = np.array([delays])
-                syn_dict.update({'weight': weights, 'delay': delays})
+                if weights is not None or delays is not None:
+                    syn_dict.update({'weight': weights, 'delay': delays})
 
                 if postsynaptic_cell.celltype.standard_receptor_type:
                     # For Tsodyks-Markram synapses models we set the "tau_psc" parameter to match
                     # the relevant "tau_syn" parameter from the post-synaptic neuron.
-                    if 'tsodyks' in self.nest_synapse_model:
+                    if 'tsodyks' in self.nest_synapse_model and hasattr(postsynaptic_cell.celltype, 'translations'):
                         translations = postsynaptic_cell.celltype.translations
                         if self.receptor_type == 'inhibitory':
                             param_name = translations['tau_syn_I']['translated_name']
diff --git a/pyNN/nest/synapses.py b/pyNN/nest/synapses.py
index c07f0c0d..9c5b5abb 100644
--- a/pyNN/nest/synapses.py
+++ b/pyNN/nest/synapses.py
@@ -40,6 +40,17 @@ def _get_nest_synapse_model(self):
             if value.is_homogeneous:
                 value.shape = (1,)
                 synapse_defaults[name] = value.evaluate(simplify=True)
+
+        if self.weight_variable and self.weight_variable != "weight":
+            w = synapse_defaults.pop("weight", None)
+            if w is not None:
+                synapse_defaults[self.weight_variable] = w
+
+        if self.delay_variable and self.delay_variable != "delay":
+            d = synapse_defaults.pop("delay", None)
+            if d is not None:
+                synapse_defaults[self.delay_variable] = d
+
         synapse_defaults = make_sli_compatible(synapse_defaults)
         synapse_defaults.pop("tau_minus", None)
         try:
diff --git a/pyNN/recording/__init__.py b/pyNN/recording/__init__.py
index 66630959..efd6c7a6 100644
--- a/pyNN/recording/__init__.py
+++ b/pyNN/recording/__init__.py
@@ -351,6 +351,8 @@ def _get_current_segment(self, filter_ids=None, variables='all', clear=False):
                 if signal_array.size > 0:
                     # may be empty if none of the recorded cells are on this MPI node
                     units = self.population.find_units(variable)
+                    if units == "unknown":
+                        units = "dimensionless"
                     channel_ids = np.fromiter(ids, dtype=int)
                     if len(ids) == signal_array.shape[1]:  # one channel per neuron
                         channel_index = np.array([self.population.id_to_index(id) for id in ids])
diff --git a/pyproject.toml b/pyproject.toml
index af831fca..b053c412 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -35,7 +35,7 @@ dependencies = [
 
 [project.optional-dependencies]
 test = ["pytest", "pytest-xdist", "pytest-cov", "flake8", "wheel", "mpi4py", "scipy", "matplotlib", "Cheetah3", "h5py", "Jinja2"]
-doc = ["sphinx"]
+doc = ["sphinx", "sphinxawesome_theme"]
 examples = ["matplotlib", "scipy"]
 plotting = ["matplotlib", "scipy"]
 MPI = ["mpi4py"]
@@ -45,6 +45,7 @@ brian2 = ["brian2"]
 arbor = ["arbor==0.9.0", "libNeuroML"]
 spiNNaker = ["spyNNaker"]
 neuroml = ["libNeuroML"]
+nestml = ["nestml"]
 
 [project.urls]
 homepage = "http://neuralensemble.org/PyNN/"
diff --git a/test/Dockerfile b/test/Dockerfile
index 3fb8f157..a7f83b77 100644
--- a/test/Dockerfile
+++ b/test/Dockerfile
@@ -51,7 +51,7 @@ COPY pyproject.toml /tmp/pynn-meta/
 RUN mkdir -p /tmp/pynn-meta/pyNN && touch /tmp/pynn-meta/pyNN/__init__.py && \
     pip install --no-cache-dir \
         "numpy<2" \
-        "/tmp/pynn-meta[test,neuron,brian2,arbor,MPI,sonata,neuroml]" && \
+        "/tmp/pynn-meta[test,neuron,brian2,arbor,MPI,sonata,neuroml,nestml]" && \
     rm -rf /tmp/pynn-meta
 
 COPY test/docker-entrypoint.sh /entrypoint.sh
diff --git a/test/environment-nest3.9.yml b/test/environment-nest3.9.yml
index d9d1b803..648d738f 100644
--- a/test/environment-nest3.9.yml
+++ b/test/environment-nest3.9.yml
@@ -21,6 +21,7 @@ channels:
 dependencies:
   - python=3.12
   - nest-simulator=3.9
+  - nestml
   - brian2
   - numpy
   - scipy
@@ -36,6 +37,7 @@ dependencies:
     - morphio
     - neuron
     - nrnutils
+    - pynestml
     - Cheetah3
     - Jinja2
     - pytest
diff --git a/test/system/test_nest_nestml.py b/test/system/test_nest_nestml.py
new file mode 100644
index 00000000..8f53c766
--- /dev/null
+++ b/test/system/test_nest_nestml.py
@@ -0,0 +1,212 @@
+import os
+import numpy as np
+import pytest
+
+try:
+    import pyNN.nest as sim
+    have_nest = True
+except ImportError:
+    have_nest = False
+
+try:
+    import pynestml  # noqa: F401  # pip install name is "nestml", import name is "pynestml"
+    have_pynestml = True
+except ImportError:
+    have_pynestml = False
+
+NESTML_MODEL_DIR = os.path.join(os.path.dirname(__file__), "..", "..", "examples", "nestml")
+
+
+@pytest.fixture(autouse=True)
+def reset_nestml_state():
+    """Reset NESTML module-level state before and after each test.
+
+    Resets both before (so state left by non-NESTML tests in the same worker doesn't bleed
+    in) and after (so state left by this test doesn't bleed into subsequent tests).
+    """
+    if have_nest:
+        from pyNN.nest import nestml
+        nestml._compiled = False
+        nestml._pending.clear()
+    yield
+    if have_nest:
+        from pyNN.nest import nestml
+        nestml._compiled = False
+        nestml._pending.clear()
+
+
+def test_nestml_cell_type_vm_trace():
+    """NESTML iaf_psc_exp_neuron V_m trace should be numerically identical to native NEST iaf_psc_exp."""
+    if not have_nest:
+        pytest.skip("nest not available")
+    if not have_pynestml:
+        pytest.skip("pynestml not available")
+
+    from pyNN.nest import nestml as pynn_nestml
+    iaf_path = os.path.join(NESTML_MODEL_DIR, "iaf_psc_exp_neuron.nestml")
+    NestmlIAF = pynn_nestml.nestml_cell_type("iaf_psc_exp_neuron", iaf_path)
+
+    sim.setup(timestep=0.1, min_delay=1.0)
+
+    nestml_pop = sim.Population(1, NestmlIAF(I_e=400.0), label="nestml")
+    native_pop = sim.Population(1, sim.native_cell_type("iaf_psc_exp")(I_e=400.0), label="native")
+    nestml_pop.record("V_m")
+    native_pop.record("V_m")
+
+    sim.run(100.0)
+
+    nestml_vm = nestml_pop.get_data().segments[0].filter(name="V_m")[0].magnitude
+    native_vm = native_pop.get_data().segments[0].filter(name="V_m")[0].magnitude
+
+    assert nestml_vm.shape == native_vm.shape
+    assert np.ptp(native_vm) > 5.0, "native iaf_psc_exp shows no dynamics — check I_e"
+    np.testing.assert_allclose(nestml_vm, native_vm, atol=1e-9,
+                               err_msg="V_m traces differ between NESTML and native iaf_psc_exp")
+
+    sim.end()
+
+
+def test_nestml_synapse_weight_changes():
+    """STDP synapse weights should change from their initial value after Poisson-driven activity."""
+    if not have_nest:
+        pytest.skip("nest not available")
+    if not have_pynestml:
+        pytest.skip("pynestml not available")
+
+    from pyNN.nest import nestml as pynn_nestml
+    iaf_path = os.path.join(NESTML_MODEL_DIR, "iaf_psc_exp_neuron.nestml")
+    stdp_path = os.path.join(NESTML_MODEL_DIR, "stdp_synapse.nestml")
+    stdp_cls = pynn_nestml.nestml_synapse_type(
+        "stdp_synapse", stdp_path,
+        postsynaptic_neuron_nestml_description=iaf_path,
+    )
+    PostCellType = stdp_cls.postsynaptic_cell_type
+
+    sim.setup(timestep=0.1, min_delay=1.0)
+
+    source = sim.Population(10, sim.SpikeSourcePoisson(rate=100.0), label="source")
+    target = sim.Population(10, PostCellType(), label="target")
+
+    initial_weight = 1.0
+    prj = sim.Projection(
+        source, target,
+        sim.AllToAllConnector(),
+        stdp_cls(weight=initial_weight, delay=1.0),
+        receptor_type="excitatory",
+    )
+
+    sim.run(1000.0)
+
+    weights = np.array(prj.get("weight", format="list"))[:, 2]
+    assert not np.allclose(weights, initial_weight), \
+        "STDP weights did not change from initial value — plasticity may not be active"
+
+    sim.end()
+
+
+def test_nestml_tsodyks_synapse_vm_trace():
+    """NESTML tsodyks_synapse postsynaptic V_m should be numerically identical to native NEST tsodyks_synapse."""
+    if not have_nest:
+        pytest.skip("nest not available")
+    if not have_pynestml:
+        pytest.skip("pynestml not available")
+
+    from pyNN.nest import nestml as pynn_nestml
+    tsodyks_path = os.path.join(NESTML_MODEL_DIR, "tsodyks_synapse.nestml")
+    TsodyksSyn = pynn_nestml.nestml_synapse_type(
+        "tsodyks_synapse_nestml", tsodyks_path,
+        weight_variable="w",
+        delay_variable="d",
+    )
+
+    sim.setup(timestep=0.1, min_delay=1.0)
+
+    spike_times = [50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0]
+    source = sim.Population(1, sim.SpikeSourceArray(spike_times=spike_times), label="source")
+
+    nestml_target = sim.Population(1, sim.native_cell_type("iaf_psc_exp")(), label="nestml_target")
+    native_target = sim.Population(1, sim.native_cell_type("iaf_psc_exp")(), label="native_target")
+
+    NativeTsodyks = sim.native_synapse_type("tsodyks_synapse")
+
+    sim.Projection(
+        source, nestml_target,
+        sim.AllToAllConnector(),
+        TsodyksSyn(weight=500.0, delay=1.0),
+        receptor_type="excitatory",
+    )
+    sim.Projection(
+        source, native_target,
+        sim.AllToAllConnector(),
+        NativeTsodyks(weight=500.0, delay=1.0, tau_psc=3.0, tau_fac=0.0, tau_rec=800.0, U=0.5),
+        receptor_type="excitatory",
+    )
+
+    nestml_target.record("V_m")
+    native_target.record("V_m")
+
+    sim.run(500.0)
+
+    nestml_vm = nestml_target.get_data().segments[0].filter(name="V_m")[0].magnitude
+    native_vm = native_target.get_data().segments[0].filter(name="V_m")[0].magnitude
+
+    assert nestml_vm.shape == native_vm.shape
+    assert np.ptp(native_vm) > 1.0, "native tsodyks_synapse target shows no response — check weight/spike_times"
+    np.testing.assert_allclose(nestml_vm, native_vm, atol=1e-6,
+                               err_msg="V_m traces differ between NESTML and native tsodyks_synapse")
+
+    sim.end()
+
+
+def test_nestml_cell_type_inline_string():
+    """nestml_cell_type() should accept inline NESTML source, not just file paths."""
+    if not have_nest:
+        pytest.skip("nest not available")
+    if not have_pynestml:
+        pytest.skip("pynestml not available")
+
+    from pyNN.nest import nestml as pynn_nestml
+    iaf_path = os.path.join(NESTML_MODEL_DIR, "iaf_psc_exp_neuron.nestml")
+    with open(iaf_path) as f:
+        iaf_source = f.read()
+
+    NestmlIAF = pynn_nestml.nestml_cell_type("iaf_psc_exp_neuron", iaf_source)
+    sim.setup(timestep=0.1, min_delay=1.0)
+
+    pop = sim.Population(1, NestmlIAF(I_e=400.0), label="nestml_inline")
+    pop.record("V_m")
+    sim.run(100.0)
+
+    vm = pop.get_data().segments[0].filter(name="V_m")[0].magnitude
+    assert np.ptp(vm) > 5.0, "No membrane potential dynamics — inline model may not have compiled"
+
+    sim.end()
+
+
+def test_nestml_register_after_setup_raises():
+    """Calling nestml_cell_type() after sim.setup() should raise RuntimeError."""
+    if not have_nest:
+        pytest.skip("nest not available")
+
+    from pyNN.nest import nestml as pynn_nestml
+    sim.setup(timestep=0.1, min_delay=1.0)
+    iaf_path = os.path.join(NESTML_MODEL_DIR, "iaf_psc_exp_neuron.nestml")
+
+    with pytest.raises(RuntimeError, match="before sim.setup"):
+        pynn_nestml.nestml_cell_type("iaf_psc_exp_neuron", iaf_path)
+
+    sim.end()
+
+
+def test_nestml_setup_without_models():
+    """sim.setup() should succeed and set _compiled even when no NESTML models are registered."""
+    if not have_nest:
+        pytest.skip("nest not available")
+
+    from pyNN.nest import nestml as pynn_nestml
+    sim.setup(timestep=0.1, min_delay=1.0)
+
+    assert pynn_nestml._compiled, "_compiled should be True after setup() even with no models"
+    assert pynn_nestml._pending == [], "_pending should be empty after setup()"
+
+    sim.end()