Fix accessing of third factor in modulatory spikes handler in synapse models (#1317)

clinssen · web-flow · commit c2669637bb5e · 2026-02-18T10:04:48.000+01:00
diff --git a/pynestml/codegeneration/resources_nest/point_neuron/common/SynapseHeader.h.jinja2 b/pynestml/codegeneration/resources_nest/point_neuron/common/SynapseHeader.h.jinja2
@@ -475,15 +475,15 @@ inline void set_{{ variable.get_name() }}(const {{ declarations.print_variable_t
    * Update internal state (``S_``) of the synapse according to the dynamical equations defined in the model and the statements in the ``update`` block.
   **/
   inline void
-  update_internal_state_(double t_start, double timestep, const {{synapseName}}CommonSynapseProperties& cp);
+  update_internal_state_(double t_start, double timestep, const {{ synapseName }}CommonSynapseProperties& cp);
 
   void recompute_internal_variables();
 
    std::string get_name() const;
 
 public:
   // this line determines which common properties to use
-  typedef {{synapseName}}CommonSynapseProperties CommonPropertiesType;
+  typedef {{ synapseName }}CommonSynapseProperties CommonPropertiesType;
 
   typedef Connection< targetidentifierT > ConnectionBase;
 
@@ -846,7 +846,7 @@ void get_entry_from_continuous_variable_history(double t,
 {#- post spike based: grab the entry from the post spiking history buffer #}
 {%-     if continuous_state_buffering_method == "post_spike_based" %}
 {%-         for var_name in paired_neuron.state_vars_that_need_continuous_buffering %}
-{%-         set var_name_post = utils.get_var_name_tuples_of_neuron_synapse_pair(continuous_post_ports, var_name) %}
+{%-             set var_name_post = utils.get_var_name_tuples_of_neuron_synapse_pair(continuous_post_ports, var_name) %}
         const double __{{ var_name }} = start->{{ var_name }}_;
 {%-         endfor %}
 {%-     endif %}
@@ -930,7 +930,7 @@ void get_entry_from_continuous_variable_history(double t,
      * update synapse internal state from `t_lastspike_` to `__t_spike`
     **/
 {%- if vt_ports is defined and vt_ports|length > 0  %}
-{%- set vt_port = vt_ports[0] %}
+{%-     set vt_port = vt_ports[0] %}
     process_{{vt_port}}_spikes_( vt_spikes, t_lastspike_, __t_spike, cp );
 {%- endif %}
 
@@ -1058,16 +1058,32 @@ constexpr ConnectionModelProperties {{synapseName}}< targetidentifierT >::proper
 {%- endif %}
 {%- if vt_ports is defined and vt_ports|length > 0  %}
 {%- set vt_port = vt_ports[0] %}
+
+/**
+ * Handler for volume transmitter spikes
+**/
 template < typename targetidentifierT >
 void
-{{synapseName}}< targetidentifierT >::process_{{vt_port}}_spikes_( const std::vector< spikecounter >& vt_spikes,
+{{synapseName}}< targetidentifierT >::process_{{ vt_port }}_spikes_( const std::vector< spikecounter >& vt_spikes,
     double t0,
     double t1,
     const {{synapseName}}CommonSynapseProperties& cp )
 {
 #ifdef DEBUG
   std::cout << "[synapse " << this << "] {{ synapseName }}::process_{{ vt_port }}_spikes_(): t0 = " << t0 << ", t1 = " << t1 << "\n";
 #endif
+
+  const size_t tid = kernel().vp_manager.get_thread_id();
+{%- if paired_neuron_name is not none and paired_neuron_name|length > 0 %}
+  typedef {{ paired_neuron_name }} post_neuron_t;
+  {{ paired_neuron_name }}* __target = static_cast< {{ paired_neuron_name }}* >(get_target(tid));
+  assert(__target);
+{%      for var_name in paired_neuron.state_vars_that_need_continuous_buffering %}
+{%-         set var_name_post = utils.get_var_name_tuples_of_neuron_synapse_pair(continuous_post_ports, var_name) %}
+  const double __{{ var_name }} = ((post_neuron_t*)(__target))->get_{{ var_name_post }}();
+{%-     endfor %}
+{%- endif %}
+
   // process dopa spikes in (t0, t1]
   // propagate weight from t0 to t1
   if ( ( vt_spikes.size() > vt_spikes_idx_ + 1 )
@@ -1529,8 +1545,9 @@ inline void
 }
 
 {%- if vt_ports is defined and vt_ports|length > 0  %}
+
 /**
- * Update to end of timestep ``t_trig``, while processing vt spikes and post spikes
+ * Update synapse internal state to end of timestep ``t_trig``, while processing vt spikes and post spikes
 **/
 template < typename targetidentifierT >
 {%- if not (nest_version.startswith("v2") or nest_version.startswith("v3.0") or nest_version.startswith("v3.1") or nest_version.startswith("v3.2") or nest_version.startswith("v3.3") or nest_version.startswith("v3.4")) %}
@@ -1554,8 +1571,15 @@ inline void
 
   const size_t tid = kernel().vp_manager.get_thread_id();
 {%- if paired_neuron_name is not none and paired_neuron_name|length > 0 %}
+  typedef {{ paired_neuron_name }} post_neuron_t;
   {{ paired_neuron_name }}* __target = static_cast< {{ paired_neuron_name }}* >(get_target(tid));
   assert(__target);
+
+{%      for var_name in paired_neuron.state_vars_that_need_continuous_buffering %}
+{%-         set var_name_post = utils.get_var_name_tuples_of_neuron_synapse_pair(continuous_post_ports, var_name) %}
+  const double __{{ var_name }} = ((post_neuron_t*)(__target))->get_{{ var_name_post }}();
+{%-     endfor %}
+
 {%- else %}
   Node* __target = get_target( tid );
 {%- endif %}
@@ -1602,7 +1626,7 @@ inline void
 {#- post spike based: grab the entry from the post spiking history buffer #}
 {%-     if continuous_state_buffering_method == "post_spike_based" %}
 {%-         for var_name in paired_neuron.state_vars_that_need_continuous_buffering %}
-{%-         set var_name_post = utils.get_var_name_tuples_of_neuron_synapse_pair(continuous_post_ports, var_name) %}
+{%-             set var_name_post = utils.get_var_name_tuples_of_neuron_synapse_pair(continuous_post_ports, var_name) %}
         const double __{{ var_name }} = start->{{ var_name }}_;
 {%-         endfor %}
 {%-     endif %}
@@ -1616,14 +1640,16 @@ inline void
 {%- filter indent(6, True) %}
 {%- if post_ports is defined %}
 {%-     for post_port in spiking_post_ports %}
+{%-         set dynamics = synapse.get_on_receive_block(post_port) %}
+{%-         if dynamics is not none %}
     /**
-     *  NESTML generated onReceive code block for postsynaptic port "{{post_port}}" begins here!
+     *  NESTML generated onReceive code block for postsynaptic port "{{ post_port }}" begins here!
     **/
 
-{%-         set dynamics = synapse.get_on_receive_block(post_port) %}
-{%-         with ast = dynamics.get_stmts_body() %}
-{%-             include "directives_cpp/StmtsBody.jinja2" %}
-{%-         endwith %}
+{%-             with ast = dynamics.get_stmts_body() %}
+{%-                 include "directives_cpp/StmtsBody.jinja2" %}
+{%-             endwith %}
+{%-         endif %}
 {%-     endfor %}
 {%- endif %}
 {%- endfilter %}
@@ -1659,7 +1685,6 @@ inline void
   t_lastspike_ = t_trig;
 }
 
-
 {%- endif %}
 
 } // namespace {{ synapseName }};
diff --git a/tests/nest_tests/resources/test_continuous_third_factor_in_modulatory_spikes_handler_neuron.nestml b/tests/nest_tests/resources/test_continuous_third_factor_in_modulatory_spikes_handler_neuron.nestml
@@ -0,0 +1,83 @@
+# test_continuous_third_factor_in_modulatory_spikes_handler_neuron
+# ################################################################
+#
+# Description
+# +++++++++++
+#
+# Like ``iaf_psc_exp``, but with a dummy postsynaptic third-factor variable.
+#
+# 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 test_continuous_third_factor_in_modulatory_spikes_handler_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
+        e_trace_dend real = 42
+
+    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/tests/nest_tests/resources/test_continuous_third_factor_in_modulatory_spikes_handler_synapse.nestml b/tests/nest_tests/resources/test_continuous_third_factor_in_modulatory_spikes_handler_synapse.nestml
@@ -0,0 +1,50 @@
+# test_continuous_third_factor_in_modulatory_spikes_handler_synapse
+# #################################################################
+#
+# 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 test_continuous_third_factor_in_modulatory_spikes_handler_synapse:
+    state:
+        w real = 1.0
+        e_trace_d_read real = 0.0
+        e_trace_d_mod real = 0.0
+
+    parameters:
+        d ms = 1 ms                 # Synaptic transmission delay
+
+    input:
+        pre_spikes <- spike
+        post_spikes <- spike
+        mod_spikes <- spike
+        e_trace_d real <- continuous
+
+    output:
+        spike(w real, d ms)
+
+    onReceive(pre_spikes):
+        e_trace_d_read = e_trace_d
+        println("e_trace_d_read: {e_trace_d_read}")
+        emit_spike(w, d)
+
+    onReceive(mod_spikes):
+        e_trace_d_mod = e_trace_d
+        println("e_trace_d_mod: {e_trace_d_mod}")
+
diff --git a/tests/nest_tests/test_continuous_third_factor_in_modulatory_spikes_handler.py b/tests/nest_tests/test_continuous_third_factor_in_modulatory_spikes_handler.py
@@ -0,0 +1,105 @@
+# -*- coding: utf-8 -*-
+#
+# test_continuous_third_factor_in_modulatory_spikes_handler.py
+#
+# 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/>.
+
+import numpy as np
+import pytest
+import re
+
+import nest
+
+from pynestml.codegeneration.nest_tools import NESTTools
+from pynestml.frontend.pynestml_frontend import generate_nest_target
+
+
+@pytest.mark.skipif(NESTTools.detect_nest_version().startswith("v2"),
+                    reason="This test does not support NEST 2")
+class TestContinuousThirdFactorInModulatorySpikesHandler:
+    r"""Test that continuous third factors can be accessed in the modulatory spikes handler in a synapse model."""
+
+    @pytest.fixture(scope="module", autouse=True)
+    def generate_all_models(self):
+        codegen_opts = {"continuous_state_buffering_method": "post_spike_based",
+                        "neuron_synapse_pairs": [{"neuron": "test_continuous_third_factor_in_modulatory_spikes_handler_neuron",
+                                                  "synapse": "test_continuous_third_factor_in_modulatory_spikes_handler_synapse",
+                                                  "vt_ports": ["mod_spikes"],
+                                                  "post_ports": ["post_spikes",
+                                                                 ["e_trace_d", "e_trace_dend"]]}],
+                        "delay_variable": {"test_continuous_third_factor_in_modulatory_spikes_handler_synapse": "d"},
+                        "weight_variable": {"test_continuous_third_factor_in_modulatory_spikes_handler_synapse": "w"}}
+
+        generate_nest_target(input_path=["tests/nest_tests/resources/test_continuous_third_factor_in_modulatory_spikes_handler_neuron.nestml",
+                                         "tests/nest_tests/resources/test_continuous_third_factor_in_modulatory_spikes_handler_synapse.nestml"],
+                             logging_level="DEBUG",
+                             suffix="_nestml",
+                             codegen_opts=codegen_opts)
+
+    def test_continuous_third_factor_in_modulatory_spikes_handler(self):
+        t_stop = 100.    # [ms]
+
+        nest.ResetKernel()
+        nest.Install("nestmlmodule")
+
+        # create spike_generator
+        vt_sg = nest.Create("poisson_generator",
+                            params={"rate": 20.})
+
+        # create volume transmitter
+        vt = nest.Create("volume_transmitter")
+        vt_parrot = nest.Create("parrot_neuron")
+        nest.Connect(vt_sg, vt_parrot)
+        sr = nest.Create("spike_recorder")
+        nest.Connect(vt_parrot, sr)
+        nest.Connect(vt_parrot, vt, syn_spec={"synapse_model": "static_synapse",
+                                              "weight": 1.,
+                                              "delay": 1.})   # delay is ignored!
+
+        nestml_model_name = "test_continuous_third_factor_in_modulatory_spikes_handler_neuron_nestml__with_test_continuous_third_factor_in_modulatory_spikes_handler_synapse_nestml"
+        synapse_model_name = "test_continuous_third_factor_in_modulatory_spikes_handler_synapse_nestml__with_test_continuous_third_factor_in_modulatory_spikes_handler_neuron_nestml"
+
+        neuron1 = nest.Create("parrot_neuron")
+        neuron2 = nest.Create(nestml_model_name)
+        nest.CopyModel(synapse_model_name, "my_nestml_synapse", {"volume_transmitter": vt})
+        nest.Connect(neuron1, neuron2, syn_spec={"synapse_model": "my_nestml_synapse", "weight": 1., "delay": 1.})
+
+        sg = nest.Create("spike_generator",
+                         params={"spike_times": [20., 50., 80.]})
+        nest.Connect(sg, neuron1, syn_spec={"synapse_model": "static_synapse", "weight": 1., "delay": 1.})
+
+        multimeter2 = nest.Create("multimeter")
+
+        V_m_specifier = "V_m"  # "delta_V_m"
+        nest.SetStatus(multimeter2, {"record_from": [V_m_specifier]})
+
+        nest.Connect(multimeter2, neuron2)
+
+        sd_pre_neuron = nest.Create("spike_recorder")
+        sd_post_neuron = nest.Create("spike_recorder")
+
+        nest.Connect(neuron1, sd_pre_neuron)
+        nest.Connect(neuron2, sd_post_neuron)
+
+        nest.Simulate(t_stop)
+
+        syn = nest.GetConnections(neuron1, neuron2)[0]
+        e_trace_d_read = syn.e_trace_d_read
+        e_trace_d_mod = syn.e_trace_d_mod
+        np.testing.assert_allclose(e_trace_d_read, 42)
+        np.testing.assert_allclose(e_trace_d_mod, 42)
