diff --git a/models/iaf_psc_alpha_multisynapse.cpp b/models/iaf_psc_alpha_multisynapse.cpp index faca17691e..cee69d2889 100644 --- a/models/iaf_psc_alpha_multisynapse.cpp +++ b/models/iaf_psc_alpha_multisynapse.cpp @@ -110,7 +110,6 @@ iaf_psc_alpha_multisynapse::Parameters_::Parameters_() iaf_psc_alpha_multisynapse::State_::State_() : I_const_( 0.0 ) , V_m_( 0.0 ) - , current_( 0.0 ) , refractory_steps_( 0 ) { y1_syn_.clear(); @@ -338,11 +337,9 @@ iaf_psc_alpha_multisynapse::update( Time const& origin, const long from, const l // neuron not refractory S_.V_m_ = V_.P30_ * ( S_.I_const_ + P_.I_e_ ) + V_.P33_ * S_.V_m_; - S_.current_ = 0.0; for ( size_t i = 0; i < P_.n_receptors_(); i++ ) { S_.V_m_ += V_.P31_syn_[ i ] * S_.y1_syn_[ i ] + V_.P32_syn_[ i ] * S_.y2_syn_[ i ]; - S_.current_ += S_.y2_syn_[ i ]; } // lower bound of membrane potential diff --git a/models/iaf_psc_alpha_multisynapse.h b/models/iaf_psc_alpha_multisynapse.h index bb372b2df2..5191d20750 100644 --- a/models/iaf_psc_alpha_multisynapse.h +++ b/models/iaf_psc_alpha_multisynapse.h @@ -214,8 +214,6 @@ class iaf_psc_alpha_multisynapse : public ArchivingNode std::vector< double > y2_syn_; //! This is the membrane potential RELATIVE TO RESTING POTENTIAL. double V_m_; - double current_; //! This is the current in a time step. This is only here - //! to allow logging int refractory_steps_; //!< Number of refractory steps remaining @@ -300,7 +298,7 @@ class iaf_psc_alpha_multisynapse : public ArchivingNode } else if ( elem == State_::I ) { - return S_.current_; + return std::accumulate( S_.y2_syn_.begin(), S_.y2_syn_.end(), 0.0 ); } else { diff --git a/models/iaf_psc_exp_multisynapse.cpp b/models/iaf_psc_exp_multisynapse.cpp index 6b22778d4b..7c53921537 100644 --- a/models/iaf_psc_exp_multisynapse.cpp +++ b/models/iaf_psc_exp_multisynapse.cpp @@ -22,7 +22,6 @@ #include "iaf_psc_exp_multisynapse.h" - // Includes from libnestutil: #include "dict_util.h" #include "exceptions.h" @@ -107,7 +106,6 @@ iaf_psc_exp_multisynapse::Parameters_::Parameters_() iaf_psc_exp_multisynapse::State_::State_() : I_const_( 0.0 ) , V_m_( 0.0 ) - , current_( 0.0 ) , refractory_steps_( 0 ) { i_syn_.clear(); @@ -313,17 +311,16 @@ iaf_psc_exp_multisynapse::update( const Time& origin, const long from, const lon { S_.V_m_ = S_.V_m_ * V_.P22_ + ( P_.I_e_ + S_.I_const_ ) * V_.P20_; // not sure about this - S_.current_ = 0.0; for ( size_t i = 0; i < P_.n_receptors_(); i++ ) { S_.V_m_ += V_.P21_syn_[ i ] * S_.i_syn_[ i ]; - S_.current_ += S_.i_syn_[ i ]; // not sure about this } } else { --S_.refractory_steps_; // neuron is absolute refractory } + for ( size_t i = 0; i < P_.n_receptors_(); i++ ) { // exponential decaying PSCs diff --git a/models/iaf_psc_exp_multisynapse.h b/models/iaf_psc_exp_multisynapse.h index 785f187a14..bf03e3dbd3 100644 --- a/models/iaf_psc_exp_multisynapse.h +++ b/models/iaf_psc_exp_multisynapse.h @@ -211,9 +211,7 @@ class iaf_psc_exp_multisynapse : public ArchivingNode double I_const_; //!< synaptic dc input current, variable 0 std::vector< double > i_syn_; - double V_m_; //!< membrane potential, variable 2 - double current_; //!< This is the current in a time step. This is only - //!< here to allow logging + double V_m_; //!< membrane potential, variable 2 //! absolute refractory counter (no membrane potential propagation) int refractory_steps_; @@ -300,7 +298,7 @@ class iaf_psc_exp_multisynapse : public ArchivingNode } else if ( elem == State_::I ) { - return S_.current_; + return std::accumulate( S_.i_syn_.begin(), S_.i_syn_.end(), 0.0 ); } else { diff --git a/testsuite/pytests/test_iaf_psc_alpha_multisynapse.py b/testsuite/pytests/sli2py_neurons/test_iaf_psc_alpha_multisynapse.py similarity index 60% rename from testsuite/pytests/test_iaf_psc_alpha_multisynapse.py rename to testsuite/pytests/sli2py_neurons/test_iaf_psc_alpha_multisynapse.py index 91b2c8c497..c222ac8112 100644 --- a/testsuite/pytests/test_iaf_psc_alpha_multisynapse.py +++ b/testsuite/pytests/sli2py_neurons/test_iaf_psc_alpha_multisynapse.py @@ -28,6 +28,7 @@ import numpy as np import numpy.testing as nptest import pytest +from scipy.linalg import expm @pytest.fixture(autouse=True) @@ -37,8 +38,7 @@ def reset(): def alpha_fn(t, tau_syn): vals = np.zeros_like(t) - zero_inds = t <= 0.0 - nonzero_inds = ~zero_inds + nonzero_inds = t > 0.0 vals[nonzero_inds] = np.e / tau_syn * t[nonzero_inds] * np.exp(-t[nonzero_inds] / tau_syn) return vals @@ -50,6 +50,16 @@ def test_I_syn_1_in_recordables(): assert "I_syn_1" in nrn.get("recordables") +def alpha_psc_voltage_response(t, tau_syn, tau_m, C_m, w): + vals = np.zeros_like(t) + nonzero_inds = t > 0.0 + A = np.array([[-1.0 / tau_syn, 0.0, 0.0], [1.0, -1.0 / tau_syn, 0.0], [0.0, 1.0 / C_m, -1.0 / tau_m]]) + + expAt = expm(A[None, ...] * t[nonzero_inds, None, None]) # shape (t, 3, 3) + vals[nonzero_inds] = expAt[:, 2, 0] * w * np.e / tau_syn # first two state variables are 0 + return vals + + def test_resize_recordables(): """ Test resizing of recordables. @@ -80,40 +90,57 @@ def test_simulation_against_analytical_soln(): from multiple different synaptic ports are the same as the analytical solution. """ - tau_syn = [2.0, 20.0, 60.0, 100.0] - delays = [100.0, 200.0, 500.0, 1200.0] - weight = 1.0 - spike_time = 10.0 - simtime = 2500.0 + tau_syns = [2.0, 20.0, 60.0, 100.0] + delays = [7.0, 5.0, 2.0, 1.0] + weights = [30.0, 50.0, 20.0, 10.0] + C_m = 250.0 + tau_m = 15.0 + spike_time = 1.0 + simtime = 100.0 + dt = 1.0 + + nest.set(resolution=dt) nrn = nest.Create( "iaf_psc_alpha_multisynapse", params={ - "C_m": 250.0, + "C_m": C_m, "E_L": 0.0, "V_m": 0.0, "V_th": 1500.0, "I_e": 0.0, - "tau_m": 15.0, - "tau_syn": tau_syn, + "tau_m": tau_m, + "tau_syn": tau_syns, }, ) - sg = nest.Create("spike_generator", params={"spike_times": [spike_time]}) - for i, syn_id in enumerate(range(1, 5)): - syn_spec = {"synapse_model": "static_synapse", "delay": delays[i], "weight": weight, "receptor_type": syn_id} + sg = nest.Create("spike_generator", params={"spike_times": [spike_time]}) + for syn_idx, (delay, weight) in enumerate(zip(delays, weights)): + syn_spec = { + "synapse_model": "static_synapse", + "delay": delay, + "weight": weight, + "receptor_type": syn_idx + 1, + } nest.Connect(sg, nrn, conn_spec="one_to_one", syn_spec=syn_spec) - mm = nest.Create("multimeter", params={"record_from": ["I_syn_1", "I_syn_2", "I_syn_3", "I_syn_4"]}) + mm = nest.Create( + "multimeter", + params={"record_from": ["I_syn_1", "I_syn_2", "I_syn_3", "I_syn_4", "V_m", "I_syn"], "interval": dt}, + ) nest.Connect(mm, nrn) nest.Simulate(simtime) times = mm.get("events", "times") - I_syn = np.sum([mm.get("events", f"I_syn_{i}") for i in range(1, 5)], axis=0) - I_syn_analytical = np.zeros_like(times, dtype=np.float64) - for i in range(4): - I_syn_analytical += alpha_fn(times - delays[i] - spike_time, tau_syn[i]) + I_syns_analytical = [] + V_m_analytical = np.zeros_like(times) + for weight, delay, tau_s in zip(weights, delays, tau_syns): + I_syns_analytical.append(alpha_fn(times - delay - spike_time, tau_s) * weight) + V_m_analytical += alpha_psc_voltage_response(times - delay - spike_time, tau_s, tau_m, C_m, weight) + + for idx, I_syn_analytical in enumerate(I_syns_analytical): + nptest.assert_array_almost_equal(mm.get("events", f"I_syn_{idx+1}"), I_syn_analytical) - nptest.assert_array_almost_equal(I_syn, I_syn_analytical) + nptest.assert_array_almost_equal(mm.get("events", "V_m"), V_m_analytical) diff --git a/testsuite/pytests/sli2py_neurons/test_iaf_psc_exp_multisynapse.py b/testsuite/pytests/sli2py_neurons/test_iaf_psc_exp_multisynapse.py new file mode 100644 index 0000000000..eb54482bef --- /dev/null +++ b/testsuite/pytests/sli2py_neurons/test_iaf_psc_exp_multisynapse.py @@ -0,0 +1,157 @@ +# -*- coding: utf-8 -*- +# +# test_iaf_psc_exp_multisynapse.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 . + +""" +Test ``iaf_psc_exp_multisynapse`` recordables and simulated PSCs against expectation. +""" + + +import nest +import numpy as np +import numpy.testing as nptest +import pytest + + +@pytest.fixture(autouse=True) +def reset(): + nest.ResetKernel() + + +def exp_psc_fn(t, tau_syn): + vals = np.zeros_like(t) + nonzero_inds = t > 0.0 + vals[nonzero_inds] = np.exp(-t[nonzero_inds] / tau_syn) + return vals + + +def exp_psc_voltage_response(t, tau_syn, tau_m, C_m, w): + vals = np.zeros_like(t) + nonzero_inds = t > 0.0 + delta_e = np.exp(-t[nonzero_inds] / tau_m) - np.exp(-t[nonzero_inds] / tau_syn) + vals[nonzero_inds] = w / (C_m * (1.0 / tau_syn - 1.0 / tau_m)) * delta_e + return vals + + +def test_set_synaptic_time_constants(): + """Tests that synaptic time constants can be set correctly""" + taus = [2.0, 20.0, 60.0, 100.0] + nrn = nest.Create("iaf_psc_exp_multisynapse") + nrn.set(tau_syn=taus) + nptest.assert_array_almost_equal(nrn.get("tau_syn"), taus) + + +def test_simulation_against_analytical_solution(): + """ + Test simulated PSCs against analytical expectation. + + This test checks that the integration of the exponential currents of inputs + from multiple different synaptic ports are the same as the analytical solution. + """ + + tau_syns = [2.0, 20.0, 60.0, 100.0] + delays = [7.0, 5.0, 2.0, 1.0] + weights = [30.0, 50.0, 20.0, 10.0] + C_m = 250.0 + tau_m = 15.0 + spike_time = 0.1 + simtime = 8.0 + dt = 0.1 + + nest.set(resolution=dt) + + nrn = nest.Create( + "iaf_psc_exp_multisynapse", + params={ + "C_m": C_m, + "E_L": 0.0, + "V_m": 0.0, + "V_th": 1500.0, + "I_e": 0.0, + "tau_m": tau_m, + "tau_syn": tau_syns, + }, + ) + + sg = nest.Create("spike_generator", params={"spike_times": [spike_time]}) + + for syn_idx, (delay, weight) in enumerate(zip(delays, weights)): + syn_spec = { + "synapse_model": "static_synapse", + "delay": delay, + "weight": weight, + "receptor_type": syn_idx + 1, + } + nest.Connect(sg, nrn, conn_spec="one_to_one", syn_spec=syn_spec) + + mm = nest.Create( + "multimeter", + params={"record_from": ["I_syn_1", "I_syn_2", "I_syn_3", "I_syn_4", "V_m", "I_syn"], "interval": dt}, + ) + + nest.Connect(mm, nrn, syn_spec={"delay": 0.1}) + nest.Simulate(simtime) + times = mm.get("events", "times") + + I_syns_analytical = [] + V_m_analytical = np.zeros_like(times) + for weight, delay, tau_s in zip(weights, delays, tau_syns): + I_syns_analytical.append(exp_psc_fn(times - delay - spike_time, tau_s) * weight) + V_m_analytical += exp_psc_voltage_response(times - delay - spike_time, tau_s, tau_m, C_m, weight) + + for idx, I_syn_analytical in enumerate(I_syns_analytical): + nptest.assert_array_almost_equal(mm.get("events", f"I_syn_{idx+1}"), I_syn_analytical) + nptest.assert_array_almost_equal(mm.get("events", "V_m"), V_m_analytical) + + +# The following tests address #800 +# - Test that the default recordables are V_m, w and I_syn_1 +# - Test that the recordable I_syn's change when changing the number of receptor ports + + +def test_default_recordables(): + nrn = nest.Create("iaf_psc_exp_multisynapse") + recordables = nrn.get("recordables") + assert len(recordables) == 3 + assert "I_syn" in recordables + assert "I_syn_1" in recordables + assert "V_m" in recordables + + +def test_resize_recordables(): + """ + Test resizing of recordables. + + This test ensures that recordables are updated correctly when the number + of synaptic ports are changed. + """ + + tau_syn1 = [5.0, 1.0, 25.0] + tau_syn2 = [5.0, 1.0] + tau_syn3 = [5.0, 1.0, 25.0, 50.0] + + nrn = nest.Create("iaf_psc_alpha_multisynapse", params={"tau_syn": tau_syn1}) + assert len(nrn.recordables) == 5 + + nrn.set(tau_syn=tau_syn2) + assert len(nrn.recordables) == 4 + + nrn.set(tau_syn=tau_syn3) + assert len(nrn.recordables) == 6 diff --git a/testsuite/unittests/test_iaf_psc_exp_multisynapse.sli b/testsuite/unittests/test_iaf_psc_exp_multisynapse.sli deleted file mode 100644 index 2284d5693f..0000000000 --- a/testsuite/unittests/test_iaf_psc_exp_multisynapse.sli +++ /dev/null @@ -1,296 +0,0 @@ -/* - * test_iaf_psc_exp_multisynapse.sli - * - * 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 . - * - */ - -/** @BeginDocumentation - Name: testsuite::test_iaf_psc_exp_multisynapse - sli script for the multisynapse model - - Synopsis: (test_iaf_psc_exp_multisynapse) run - - Description: - This test creates a multisynapse neuron and first checks if time constants - can be set correctly. - Afterwards a neuron with four synaptic time constants with different weights and delays - is simulated and the resulting membrane potential values after the income of one spike - are checked against analytically obtained values. - - SeeAlso: testsuite::test_iaf_psc_alpha_multisynapse, iaf_psc_exp_multisynapse - - FirstVersion: April 2013 - Author: Hannah Bos - */ - - -(unittest) run -/unittest using - -0.1 /h Set - -ResetKernel - -<< - /local_num_threads 1 - /resolution h - >> SetKernelStatus - -/taus [2. 20. 60. 100.] def -/delays [7. 5. 2. 1. ] def % ms -/weights [30. 50. 20. 10.] def -/spike_time 0.1 def -/dt 0.1 def - -/iaf_psc_exp_multisynapse Create /npost Set -npost << /tau_syn taus >> SetStatus - -npost 0 get GetStatus /tau_syn get -taus eq assert_or_die - -/voltmeter Create /vm Set -vm << /time_in_steps true /interval h >> SetStatus - -/spike_generator Create /sg Set -sg << /spike_times [spike_time] >> SetStatus - -[ weights delays [ 4 ] Range ] -{ - /receptor Set - /delay Set - /weight Set - sg npost /one_to_one pstack << /weight weight /delay delay /receptor_type receptor >> Connect -} ScanThread - - -vm npost 1.0 h Connect - -8 ms Simulate - - -{ % reference data - dup Transpose First /test_times Set % times of reference - - vm [/events [/times /V_m]] get cva % array of recorded data - 6 ToUnitTestPrecision % to precision of reference - Transpose % all recorded tuples - {First test_times exch MemberQ } Select % those with reference - eq % compare -} - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% Expected output of this program: -% -% time (in steps) voltage (in mV) -[ -[ 1 -70 ] % <----- Spike time -[ 2 -70 ] -[ 3 -70 ] -[ 4 -70 ] -[ 5 -70 ] -[ 6 -70 ] -[ 7 -70 ] -[ 8 -70 ] -[ 9 -70 ] -[ 10 -70 ] -[ 11 -70 ] % <----- The spike arrives at the first synapse but has not yet effected the synaptic current and membrane potential. -[ 12 -69.9960 ] % <- -[ 13 -69.9921 ] % | -[ 14 -69.9882 ] % | -[ 15 -69.9843 ] % --- the effect of the synaptic current is visible in the membrane potential -% -% ... -% -[ 20 -69.9657 ] %<----- The spike arrives at the second synapse but has not yet effected the snynaptic current and membrane potential. -[ 21 -69.9621 ] %<- -[ 22 -69.9506 ] % | -[ 23 -69.9392 ] % | -[ 24 -69.928 ] % --- the effect of the two synaptic currents is visible in the membrane potential -[ 25 -69.9169 ] % -% -% ... -% -[ 50 -69.6771 ] %<----- The spike arrives at the third synapse but has not yet effected the snynaptic current and membrane potential. -[ 51 -69.6689 ] %<- -[ 52 -69.641 ] % | -[ 53 -69.6134 ] % | -[ 54 -69.5863 ] % --- the effect of the three synaptic currents is visible in the membrane potential -[ 55 -69.5595 ] % -% -% ... -% -[ 70 -69.2011 ] %<----- The spike arrives at the fourth synapse but has not yet effected the snynaptic current and membrane potential. -[ 71 -69.18 ] %<- -[ 72 -69.1474 ] % | -[ 73 -69.1159 ] % | -[ 74 -69.0854 ] % --- the effect of the four synaptic currents is visible in the membrane potential -[ 75 -69.0557 ] % -] -% -% - -exch assert_or_die - -% -------------------------------------------- - -% The following tests address #800 -% - Test that the default recordables are V_m, w and I_syn_1 -% - Test that the recordable I_syn's change when changing the number of receptor ports - -% test default recordables include I_syn_1 -{ - << >> begin - ResetKernel - - /nrn /iaf_psc_exp_multisynapse Create def - - /mm /multimeter << /time_in_steps true - /interval 1.0 - /record_from [ /V_m /I_syn /I_syn_1 ] - >> Create def - - mm nrn Connect - - /len nrn /recordables get length def - len 3 eq - end -} -assert_or_die - -% test resize recordables -{ - << >> begin - ResetKernel - - /tau_syn1 [5.0 1.0 25.0] def - /tau_syn2 [5.0 1.0] def - /tau_syn3 [5.0 1.0 25.0 50.] def - - /nrn /iaf_psc_exp_multisynapse << /tau_syn tau_syn1 >> Create def - /len1 nrn /recordables get length def - - nrn << /tau_syn tau_syn2 >> SetStatus - /len2 nrn /recordables get length def - - nrn << /tau_syn tau_syn3 >> SetStatus - /len3 nrn /recordables get length def - - len1 5 eq len2 4 eq and len3 6 eq and - end -} -assert_or_die - -% test record I_syn_i and check for exp function synapse -{ - << >> begin - ResetKernel - /tau_syn [40.0 20.0 30.0 25.] def % synaptic time constants - /weight [1.0 0.5 2.0 1.0] def % synaptic weights - /delays [1.0 3.0 10.0 10.] def % ms - synaptic delays - /spike_time 10. def % time at which the single spike occurs - /dt 0.1 def % time step - /total_t 500. def % total simulation time - /tolerance 1e-7 def % tolerable difference between theoretic and simulated alpha synapse currents - - /exp_function - { - % Call like t0 W tau t alpha_function - << /tau 1.0 /W 1.0 /t0 0. >> - begin - /t exch def - /tau exch def - /W exch def - /t0 exch def - - t - { - /tt exch def - tt t0 geq - { - /tdiff_over_tau tt t0 sub tau div def - tdiff_over_tau neg exp W mul - } - % else - { - 0. - } ifelse - } Map - end - } def - - << /resolution dt >> SetKernelStatus - - % Create a spike generator that generates a single spike - /spike_generator Create /sg Set - sg << /spike_times [spike_time] >> SetStatus % generates a single peak - - % Create the multisynapse neuron - /nrn /iaf_psc_exp_multisynapse - << /I_e 0. /tau_syn tau_syn >> Create def - - % Create an array of synaptic indexes to loop through - delays length 1 arraystore Range /synapses_idx exch def - [delays weights synapses_idx] % loop on synaptic receptor ports - { - /syn_id exch def - /W exch def - /delay exch def - % Connect spike generator to each port - sg nrn /one_to_one << - /synapse_model /static_synapse - /delay delay - /weight W - /receptor_type syn_id >> - Connect - } ScanThread - - % Create the multimeter that will record from the 4 synapse channels - /mm /multimeter << /time_in_steps true - /interval dt - /record_from [ /I_syn_1 /I_syn_2 /I_syn_3 /I_syn_4 ] - >> Create def - - mm nrn Connect - - % Simulate - total_t Simulate - - % Get the conductances measured during the simulation - /t mm /events get /times get cva dt mul def - /sim_I_syn_1 mm /events get /I_syn_1 get cva def - /sim_I_syn_2 mm /events get /I_syn_2 get cva def - /sim_I_syn_3 mm /events get /I_syn_3 get cva def - /sim_I_syn_4 mm /events get /I_syn_4 get cva def - /sim_I_syns [sim_I_syn_1 sim_I_syn_2 sim_I_syn_3 sim_I_syn_4] def - - true - [delays weights tau_syn sim_I_syns] - { - /sim_I_syn exch def - /tau exch def - /W exch def - /t0 exch spike_time add def - /theo_I_syn t0 W tau t exp_function def - sim_I_syn theo_I_syn sub { abs } Map Max tolerance leq and - } ScanThread - end -} -assert_or_die - -endusing