diff --git a/models/synapses/tsodyks_synapse.nestml b/models/synapses/tsodyks_synapse.nestml new file mode 100644 index 000000000..ab37a367e --- /dev/null +++ b/models/synapses/tsodyks_synapse.nestml @@ -0,0 +1,100 @@ +""" +tsodyks_synapse - Synapse type with short term plasticity +######################################################### + +Description ++++++++++++ + +This synapse model implements synaptic short-term depression and short-term +facilitation according to [1]_. In particular it solves Eqs (3) and (4) from +this paper in an exact manner. + +Synaptic depression is motivated by depletion of vesicles in the readily +releasable pool of synaptic vesicles (variable x in equation (3)). Synaptic +facilitation comes about by a presynaptic increase of release probability, +which is modeled by variable U in Eq (4). + +The original interpretation of variable y is the amount of glutamate +concentration in the synaptic cleft. In [1]_ this variable is taken to be +directly proportional to the synaptic current caused in the postsynaptic +neuron (with the synaptic weight w as a proportionality constant). In order +to reproduce the results of [1]_ and to use this model of synaptic plasticity +in its original sense, the user therefore has to ensure the following +conditions: + +1.) The postsynaptic neuron must be of type ``iaf_psc_exp`` or ``iaf_psc_exp_htum``, +because these neuron models have a postsynaptic current which decays +exponentially. + +2.) The time constant of each ``tsodyks_synapse`` targeting a particular neuron +must be chosen equal to that neuron's synaptic time constant. In particular +that means that all synapses targeting a particular neuron have the same +parameter ``tau_psc``. + +However, there are no technical restrictions using this model of synaptic +plasticity also in conjunction with neuron models that have a different +dynamics for their synaptic current or conductance. The effective synaptic +weight, which will be transmitted to the postsynaptic neuron upon occurrence +of a spike at time t is :math:`u(t) \cdot x(t) \cdot w`, where ``u(t)`` and ``x(t)`` +are defined in Eq (3) and (4), w is the synaptic weight specified upon connection. +The interpretation is as follows: The quantity :math:`u(t) \cdot x(t)` is the release +probability times the amount of releasable synaptic vesicles at time t of the +presynaptic neuron's spike, so this equals the amount of transmitter expelled +into the synaptic cleft. + +The amount of transmitter then relaxes back to 0 with time constant tau_psc +of the synapse's variable y. Since the dynamics of ``y(t)`` is linear, the +postsynaptic neuron can reconstruct from the amplitude of the synaptic +impulse :math:`u(t) \cdot x(t) \cdot w` the full shape of ``y(t)``. The postsynaptic +neuron, however, might choose to have a synaptic current that is not necessarily +identical to the concentration of transmitter ``y(t)`` in the synaptic cleft. It may +realize an arbitrary postsynaptic effect depending on ``y(t)``. + + +References +++++++++++ + +.. [1] Tsodyks M, Uziel A, Markram H (2000). Synchrony generation in recurrent + networks with frequency-dependent synapses. Journal of Neuroscience, + 20 RC50. URL: http://infoscience.epfl.ch/record/183402 +""" +model tsodyks_synapse: + state: + x real = 1 + y real = 0 + u real = 0 + + parameters: + w real = 1 @nest::weight # Synaptic weight + d ms = 1 ms @nest::delay # Synaptic transmission delay + tau_psc ms = 3 ms + tau_fac ms = 0 ms + tau_rec ms = 800 ms + U real = .5 + + input: + pre_spikes <- spike + + output: + spike + + onReceive(pre_spikes): + Puu real = tau_fac == 0 ? 0 : exp(-resolution() / tau_fac) + Pyy real = exp(-resolution() / tau_psc) + Pzz real = exp(-resolution() / 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 + + # depress synapse + 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 + + # deliver spike to postsynaptic partner + emit_spike(delta_y_tsp * w, d) diff --git a/pynestml/codegeneration/resources_nest/point_neuron/common/SynapseHeader.h.jinja2 b/pynestml/codegeneration/resources_nest/point_neuron/common/SynapseHeader.h.jinja2 index 3d727ea23..3554de941 100644 --- a/pynestml/codegeneration/resources_nest/point_neuron/common/SynapseHeader.h.jinja2 +++ b/pynestml/codegeneration/resources_nest/point_neuron/common/SynapseHeader.h.jinja2 @@ -671,14 +671,23 @@ public: send( Event& e, const thread tid, const {{synapseName}}CommonSynapseProperties& cp ) {%- endif %} { - const double __resolution = nest::Time::get_resolution().get_ms(); // do not remove, this is necessary for the resolution() function + const double __t_spike = e.get_stamp().get_ms(); + + double __resolution; // do not remove, this is necessary for the resolution() function + if (t_lastspike_ < 0.) + { + __resolution = __t_spike; + } + else + { + __resolution = __t_spike - t_lastspike_; + } auto get_thread = [tid]() { return tid; }; - const double __t_spike = e.get_stamp().get_ms(); #ifdef DEBUG std::cout << "{{synapseName}}::send(): handling pre spike at t = " << __t_spike << std::endl; #endif diff --git a/tests/nest_tests/tsodyks_synapse_test.py b/tests/nest_tests/tsodyks_synapse_test.py new file mode 100644 index 000000000..59d1850f3 --- /dev/null +++ b/tests/nest_tests/tsodyks_synapse_test.py @@ -0,0 +1,324 @@ +# -*- coding: utf-8 -*- +# +# tsodyks_synapse_test.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 . + +import numpy as np +import os +import pytest + +import nest + +from pynestml.codegeneration.nest_tools import NESTTools +from pynestml.frontend.pynestml_frontend import generate_nest_target + +try: + import matplotlib + matplotlib.use("Agg") + import matplotlib.ticker + import matplotlib.pyplot as plt + TEST_PLOTS = True +except Exception: + TEST_PLOTS = False + + +class TestNESTTsodyksSynapse: + + neuron_model_name = "iaf_psc_exp_neuron_nestml__with_tsodyks_synapse_nestml" + ref_neuron_model_name = "iaf_psc_exp_neuron_nestml_non_jit" + + synapse_model_name = "tsodyks_synapse_nestml__with_iaf_psc_exp_neuron_nestml" + ref_synapse_model_name = "tsodyks_synapse" + + @pytest.fixture(scope="module", autouse=True) + def setUp(self): + """Generate the model code""" + + jit_codegen_opts = {"neuron_synapse_pairs": [{"neuron": "iaf_psc_exp_neuron", + "synapse": "tsodyks_synapse"}]} + if not NESTTools.detect_nest_version().startswith("v2"): + jit_codegen_opts["neuron_parent_class"] = "StructuralPlasticityNode" + jit_codegen_opts["neuron_parent_class_include"] = "structural_plasticity_node.h" + + # generate the "jit" model (co-generated neuron and synapse), that does not rely on ArchivingNode + files = [os.path.join("models", "neurons", "iaf_psc_exp_neuron.nestml"), + os.path.join("models", "synapses", "tsodyks_synapse.nestml")] + input_path = [os.path.realpath(os.path.join(os.path.dirname(__file__), os.path.join( + os.pardir, os.pardir, s))) for s in files] + generate_nest_target(input_path=input_path, + target_path="/tmp/nestml-jit", + logging_level="INFO", + module_name="nestml_jit_module", + suffix="_nestml", + codegen_opts=jit_codegen_opts) + + if NESTTools.detect_nest_version().startswith("v2"): + non_jit_codegen_opts = {"neuron_parent_class": "Archiving_Node", + "neuron_parent_class_include": "archiving_node.h"} + else: + non_jit_codegen_opts = {"neuron_parent_class": "ArchivingNode", + "neuron_parent_class_include": "archiving_node.h"} + + # generate the "non-jit" model, that relies on ArchivingNode + generate_nest_target(input_path=os.path.realpath(os.path.join(os.path.dirname(__file__), + os.path.join(os.pardir, os.pardir, "models", "neurons", "iaf_psc_exp_neuron.nestml"))), + target_path="/tmp/nestml-non-jit", + logging_level="INFO", + module_name="nestml_non_jit_module", + suffix="_nestml_non_jit", + codegen_opts=non_jit_codegen_opts) + + def test_nest_tsodyks_synapse(self): + fname_snip = "" + + pre_spike_times = [1., 11., 21.] # [ms] + post_spike_times = [6., 16., 26.] # [ms] + + post_spike_times = np.sort(np.unique(1 + np.round(10 * np.sort(np.abs(np.random.randn(10)))))) # [ms] + pre_spike_times = np.sort(np.unique(1 + np.round(10 * np.sort(np.abs(np.random.randn(10)))))) # [ms] + + #post_spike_times = np.sort(np.unique(1 + np.round(100 * np.sort(np.abs(np.random.randn(100)))))) # [ms] + #pre_spike_times = np.sort(np.unique(1 + np.round(100 * np.sort(np.abs(np.random.randn(100)))))) # [ms] + + """pre_spike_times = np.array([2., 4., 7., 8., 12., 13., 19., 23., 24., 28., 29., 30., 33., 34., + 35., 36., 38., 40., 42., 46., 51., 53., 54., 55., 56., 59., 63., 64., + 65., 66., 68., 72., 73., 76., 79., 80., 83., 84., 86., 87., 90., 95., + 99., 100., 103., 104., 105., 111., 112., 126., 131., 133., 134., 139., 147., 150., + 152., 155., 172., 175., 176., 181., 196., 197., 199., 202., 213., 215., 217., 265.]) + post_spike_times = np.array([4., 5., 6., 7., 10., 11., 12., 16., 17., 18., 19., 20., 22., 23., + 25., 27., 29., 30., 31., 32., 34., 36., 37., 38., 39., 42., 44., 46., + 48., 49., 50., 54., 56., 57., 59., 60., 61., 62., 67., 74., 76., 79., + 80., 81., 83., 88., 93., 94., 97., 99., 100., 105., 111., 113., 114., 115., + 116., 119., 123., 130., 132., 134., 135., 145., 152., 155., 158., 166., 172., 174., + 188., 194., 202., 245., 249., 289., 454.])""" + + self.run_synapse_test(neuron_model_name=self.neuron_model_name, + ref_neuron_model_name=self.ref_neuron_model_name, + synapse_model_name=self.synapse_model_name, + ref_synapse_model_name=self.ref_synapse_model_name, + resolution=.5, # [ms] + delay=1.5, # [ms] + pre_spike_times=pre_spike_times, + post_spike_times=post_spike_times, + fname_snip=fname_snip) + + def run_synapse_test(self, neuron_model_name, + ref_neuron_model_name, + synapse_model_name, + ref_synapse_model_name, + resolution=1., # [ms] + delay=1., # [ms] + sim_time=None, # if None, computed from pre and post spike times + pre_spike_times=None, + post_spike_times=None, + fname_snip=""): + + if pre_spike_times is None: + pre_spike_times = [] + + if post_spike_times is None: + post_spike_times = [] + + if sim_time is None: + sim_time = max(np.amax(pre_spike_times), np.amax(post_spike_times)) + 5 * delay + + nest.ResetKernel() + nest.set_verbosity("M_ALL") + nest.Install("nestml_jit_module") + nest.Install("nestml_non_jit_module") + nest.SetKernelStatus({"resolution": resolution}) + + print("Pre spike times: " + str(pre_spike_times)) + print("Post spike times: " + str(post_spike_times)) + + wr = nest.Create("weight_recorder") + wr_ref = nest.Create("weight_recorder") + nest.CopyModel(synapse_model_name, "tsodyks_nestml_rec", + {"weight_recorder": wr[0], "w": 1., "d": 1., "receptor_type": 0}) + nest.CopyModel(ref_synapse_model_name, "tsodyks_ref_rec", + {"weight_recorder": wr_ref[0], "weight": 1., "delay": 1., "receptor_type": 0}) + + # create spike_generators with these times + pre_sg = nest.Create("spike_generator", + params={"spike_times": pre_spike_times}) + post_sg = nest.Create("spike_generator", + params={"spike_times": post_spike_times, + "allow_offgrid_times": True}) + + # create parrot neurons and connect spike_generators + pre_neuron = nest.Create("parrot_neuron") + post_neuron = nest.Create(neuron_model_name) + + pre_neuron_ref = nest.Create("parrot_neuron") + post_neuron_ref = nest.Create(ref_neuron_model_name) + + if NESTTools.detect_nest_version().startswith("v2"): + spikedet_pre = nest.Create("spike_detector") + spikedet_post = nest.Create("spike_detector") + else: + spikedet_pre = nest.Create("spike_recorder") + spikedet_post = nest.Create("spike_recorder") + mm = nest.Create("multimeter", params={"record_from": ["V_m"]}) + + if NESTTools.detect_nest_version().startswith("v2"): + spikedet_pre_ref = nest.Create("spike_detector") + spikedet_post_ref = nest.Create("spike_detector") + else: + spikedet_pre_ref = nest.Create("spike_recorder") + spikedet_post_ref = nest.Create("spike_recorder") + mm_ref = nest.Create("multimeter", params={"record_from": ["V_m"]}) + + if True: + nest.Connect(pre_sg, pre_neuron, "one_to_one", syn_spec={"delay": 1.}) + nest.Connect(post_sg, post_neuron, "one_to_one", syn_spec={"delay": 1., "weight": 9999.}) + if NESTTools.detect_nest_version().startswith("v2"): + nest.Connect(pre_neuron, post_neuron, "all_to_all", syn_spec={"model": "tsodyks_nestml_rec"}) + else: + nest.Connect(pre_neuron, post_neuron, "all_to_all", syn_spec={"synapse_model": "tsodyks_nestml_rec"}) + nest.Connect(mm, post_neuron) + nest.Connect(pre_neuron, spikedet_pre) + nest.Connect(post_neuron, spikedet_post) + + nest.Connect(pre_sg, pre_neuron_ref, "one_to_one", syn_spec={"delay": 1.}) + nest.Connect(post_sg, post_neuron_ref, "one_to_one", syn_spec={"delay": 1., "weight": 9999.}) + if NESTTools.detect_nest_version().startswith("v2"): + nest.Connect(pre_neuron_ref, post_neuron_ref, "all_to_all", + syn_spec={"model": ref_synapse_model_name}) + else: + nest.Connect(pre_neuron_ref, post_neuron_ref, "all_to_all", + syn_spec={"synapse_model": ref_synapse_model_name}) + nest.Connect(mm_ref, post_neuron_ref) + nest.Connect(pre_neuron_ref, spikedet_pre_ref) + nest.Connect(post_neuron_ref, spikedet_post_ref) + + # get tsodyks synapse and weight before protocol + syn = nest.GetConnections(source=pre_neuron, synapse_model="tsodyks_nestml_rec") + syn_ref = nest.GetConnections(source=pre_neuron_ref, synapse_model=ref_synapse_model_name) + + n_steps = int(np.ceil(sim_time / resolution)) + 1 + t = 0. + t_hist = [] + hist = {} + hist_ref = {} + for key in ["x", "y", "u"]: + hist[key] = [] + hist_ref[key] = [] + + while t <= sim_time: + nest.Simulate(resolution) + t += resolution + t_hist.append(t) + for key in ["x", "y", "u"]: + hist_ref[key].append(nest.GetStatus(syn_ref)[0][key]) + hist[key].append(nest.GetStatus(syn)[0][key]) + + # plot + if TEST_PLOTS: + fig, ax = plt.subplots(nrows=2) + ax1, ax2 = ax + + if True: + timevec = nest.GetStatus(mm, "events")[0]["times"] + V_m = nest.GetStatus(mm, "events")[0]["V_m"] + #ax2.plot(timevec, nest.GetStatus(mm, "events")[0]["post_trace__for_tsodyks_nestml"], label="post_tr nestml") + ax1.plot(timevec, V_m, label="nestml", alpha=.7, linestyle=":") + + pre_ref_spike_times_ = nest.GetStatus(spikedet_pre_ref, "events")[0]["times"] + timevec = nest.GetStatus(mm_ref, "events")[0]["times"] + V_m = nest.GetStatus(mm_ref, "events")[0]["V_m"] + ax1.plot(timevec, V_m, label="nest ref", alpha=.7) + ax1.set_ylabel("V_m") + + for _ax in ax: + _ax.grid(which="major", axis="both") + _ax.grid(which="minor", axis="x", linestyle=":", alpha=.4) + # _ax.minorticks_on() + _ax.set_xlim(0., sim_time) + _ax.legend() + fig.savefig("/tmp/tsodyks_synapse_test" + fname_snip + "_V_m.png", dpi=300) + + # plot + if TEST_PLOTS: + fig, ax = plt.subplots(nrows=5) + ax1, ax2, ax3, ax4, ax5 = ax + + pre_spike_times_ = nest.GetStatus(spikedet_pre, "events")[0]["times"] + print("Actual pre spike times: " + str(pre_spike_times_)) + pre_ref_spike_times_ = nest.GetStatus(spikedet_pre_ref, "events")[0]["times"] + print("Actual pre ref spike times: " + str(pre_ref_spike_times_)) + + n_spikes = len(pre_spike_times_) + for i in range(n_spikes): + if i == 0: + _lbl = "nestml" + else: + _lbl = None + ax1.plot(2 * [pre_spike_times_[i] + delay], [0, 1], linewidth=2, color="blue", alpha=.4, label=_lbl) + + post_spike_times_ = nest.GetStatus(spikedet_post, "events")[0]["times"] + print("Actual post spike times: " + str(post_spike_times_)) + post_ref_spike_times_ = nest.GetStatus(spikedet_post_ref, "events")[0]["times"] + print("Actual post ref spike times: " + str(post_ref_spike_times_)) + + if True: + n_spikes = len(pre_ref_spike_times_) + for i in range(n_spikes): + if i == 0: + _lbl = "nest ref" + else: + _lbl = None + ax1.plot(2 * [pre_ref_spike_times_[i] + delay], [0, 1], + linewidth=2, color="cyan", label=_lbl, alpha=.4) + ax1.set_ylabel("Pre spikes") + + if True: + n_spikes = len(post_spike_times_) + for i in range(n_spikes): + if i == 0: + _lbl = "nestml" + else: + _lbl = None + ax2.plot(2 * [post_spike_times_[i]], [0, 1], linewidth=2, color="black", alpha=.4, label=_lbl) + + n_spikes = len(post_ref_spike_times_) + for i in range(n_spikes): + if i == 0: + _lbl = "nest ref" + else: + _lbl = None + ax2.plot(2 * [post_ref_spike_times_[i]], [0, 1], linewidth=2, color="red", alpha=.4, label=_lbl) + #ax2.plot(timevec, nest.GetStatus(mm, "events")[0]["post_⎄trace__for_tsodyks_nestml"], label="nestml post tr") + ax2.set_ylabel("Post spikes") + + for i, key in enumerate(["x", "y", "u"]): + ax[2 + i].plot(t_hist, hist[key], marker="o", label="nestml") + ax[2 + i].plot(t_hist, hist_ref[key], linestyle="--", marker="x", label="ref") + ax[2 + i].set_ylabel(key) + + ax[-1].set_xlabel("Time [ms]") + + for _ax in ax: + _ax.grid(which="major", axis="both") + _ax.xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.arange(0, np.ceil(sim_time)))) + _ax.set_xlim(0., sim_time) + _ax.legend() + fig.savefig("/tmp/tsodyks_synapse_test" + fname_snip + ".png", dpi=300) + + # verify + # ...