diff --git a/models/CMakeLists.txt b/models/CMakeLists.txt index 24461cbe4b..62a0813b19 100644 --- a/models/CMakeLists.txt +++ b/models/CMakeLists.txt @@ -49,6 +49,7 @@ set( models_sources gif_pop_psc_exp.h gif_pop_psc_exp.cpp ginzburg_neuron.h ginzburg_neuron.cpp hh_cond_exp_traub.h hh_cond_exp_traub.cpp + hh_cond_beta_gap_traub.h hh_cond_beta_gap_traub.cpp hh_psc_alpha.h hh_psc_alpha.cpp hh_psc_alpha_clopath.h hh_psc_alpha_clopath.cpp hh_psc_alpha_gap.h hh_psc_alpha_gap.cpp diff --git a/models/hh_cond_beta_gap_traub.cpp b/models/hh_cond_beta_gap_traub.cpp new file mode 100644 index 0000000000..be02e92463 --- /dev/null +++ b/models/hh_cond_beta_gap_traub.cpp @@ -0,0 +1,784 @@ +/* + * hh_cond_beta_gap_traub.cpp + * + * 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 . + * + */ + + +#include "hh_cond_beta_gap_traub.h" + +#ifdef HAVE_GSL + +// C++ includes: +#include // in case we need isnan() // fabs +#include +#include +#include +#include + +// External includes: +#include +#include +#include + +// Includes from libnestutil: +#include "numerics.h" + +// Includes from nestkernel: +#include "exceptions.h" +#include "kernel_manager.h" +#include "universal_data_logger_impl.h" + +// Includes from sli: +#include "dict.h" +#include "dictutils.h" +#include "doubledatum.h" +#include "integerdatum.h" + +nest::RecordablesMap< nest::hh_cond_beta_gap_traub > + nest::hh_cond_beta_gap_traub::recordablesMap_; + +namespace nest +{ +// Override the create() method with one call to RecordablesMap::insert_() +// for each quantity to be recorded. +template <> +void +RecordablesMap< hh_cond_beta_gap_traub >::create() +{ + // use standard names whereever you can for consistency! + insert_( names::V_m, + &hh_cond_beta_gap_traub:: + get_y_elem_< hh_cond_beta_gap_traub::State_::V_M > ); + insert_( names::g_ex, + &hh_cond_beta_gap_traub:: + get_y_elem_< hh_cond_beta_gap_traub::State_::G_EXC > ); + insert_( names::g_in, + &hh_cond_beta_gap_traub:: + get_y_elem_< hh_cond_beta_gap_traub::State_::G_INH > ); + insert_( names::Act_m, + &hh_cond_beta_gap_traub:: + get_y_elem_< hh_cond_beta_gap_traub::State_::HH_M > ); + insert_( names::Act_h, + &hh_cond_beta_gap_traub:: + get_y_elem_< hh_cond_beta_gap_traub::State_::HH_H > ); + insert_( names::Inact_n, + &hh_cond_beta_gap_traub:: + get_y_elem_< hh_cond_beta_gap_traub::State_::HH_N > ); +} + +extern "C" int +hh_cond_beta_gap_traub_dynamics( double time, + const double y[], + double f[], + void* pnode ) +{ + // a shorthand + typedef nest::hh_cond_beta_gap_traub::State_ S; + + // get access to node so we can almost work as in a member function + assert( pnode ); + const nest::hh_cond_beta_gap_traub& node = + *( reinterpret_cast< nest::hh_cond_beta_gap_traub* >( pnode ) ); + + // y[] here is---and must be---the state vector supplied by the integrator, + // not the state vector in the node, node.S_.y[]. + + // The following code is verbose for the sake of clarity. We assume that a + // good compiler will optimize the verbosity away ... + + // ionic currents + const double I_Na = node.P_.g_Na * y[ S::HH_M ] * y[ S::HH_M ] * y[ S::HH_M ] + * y[ S::HH_H ] * ( y[ S::V_M ] - node.P_.E_Na ); + const double I_K = node.P_.g_K * y[ S::HH_N ] * y[ S::HH_N ] * y[ S::HH_N ] + * y[ S::HH_N ] * ( y[ S::V_M ] - node.P_.E_K ); + const double I_L = node.P_.g_L * ( y[ S::V_M ] - node.P_.E_L ); + + // chemical synaptic currents + const double I_syn_exc = y[ S::G_EXC ] * ( y[ S::V_M ] - node.P_.E_ex ); + const double I_syn_inh = y[ S::G_INH ] * ( y[ S::V_M ] - node.P_.E_in ); + + // gap junction currents + // set I_gap depending on interpolation order + double gap = 0.0; + + const double t = time / node.B_.step_; + + switch ( kernel().simulation_manager.get_wfr_interpolation_order() ) + { + case 0: + gap = -node.B_.sumj_g_ij_ * y[ S::V_M ] + + node.B_.interpolation_coefficients[ node.B_.lag_ ]; + break; + + case 1: + gap = -node.B_.sumj_g_ij_ * y[ S::V_M ] + + node.B_.interpolation_coefficients[ node.B_.lag_ * 2 + 0 ] + + node.B_.interpolation_coefficients[ node.B_.lag_ * 2 + 1 ] * t; + break; + + case 3: + gap = -node.B_.sumj_g_ij_ * y[ S::V_M ] + + node.B_.interpolation_coefficients[ node.B_.lag_ * 4 + 0 ] + + node.B_.interpolation_coefficients[ node.B_.lag_ * 4 + 1 ] * t + + node.B_.interpolation_coefficients[ node.B_.lag_ * 4 + 2 ] * t * t + + node.B_.interpolation_coefficients[ node.B_.lag_ * 4 + 3 ] * t * t * t; + break; + + default: + throw BadProperty( "Interpolation order must be 0, 1, or 3." ); + } + + const double I_gap = gap; + + // membrane potential + f[ S::V_M ] = ( -I_Na - I_K - I_L - I_syn_exc - I_syn_inh + node.B_.I_stim_ + + I_gap + node.P_.I_e ) / node.P_.C_m; + + // channel dynamics + const double V = y[ S::V_M ] - node.P_.V_T; + + const double alpha_n = + 0.032 * ( 15. - V ) / ( std::exp( ( 15. - V ) / 5. ) - 1. ); + const double beta_n = 0.5 * std::exp( ( 10. - V ) / 40. ); + const double alpha_m = + 0.32 * ( 13. - V ) / ( std::exp( ( 13. - V ) / 4. ) - 1. ); + const double beta_m = + 0.28 * ( V - 40. ) / ( std::exp( ( V - 40. ) / 5. ) - 1. ); + const double alpha_h = 0.128 * std::exp( ( 17. - V ) / 18. ); + const double beta_h = 4. / ( 1. + std::exp( ( 40. - V ) / 5. ) ); + + f[ S::HH_M ] = alpha_m - ( alpha_m + beta_m ) * y[ S::HH_M ]; // m-variable + f[ S::HH_H ] = alpha_h - ( alpha_h + beta_h ) * y[ S::HH_H ]; // h-variable + f[ S::HH_N ] = alpha_n - ( alpha_n + beta_n ) * y[ S::HH_N ]; // n-variable + + // synapses: beta function + // d^2g_exc/dt^2, dg_exc/dt + f[ S::DG_EXC ] = -y[ S::DG_EXC ] / node.P_.tau_decay_ex; + f[ S::G_EXC ] = y[ S::DG_EXC ] - ( y[ S::G_EXC ] / node.P_.tau_rise_ex ); + + // d^2g_inh/dt^2, dg_inh/dt + f[ S::DG_INH ] = -y[ S::DG_INH ] / node.P_.tau_decay_in; + f[ S::G_INH ] = y[ S::DG_INH ] - ( y[ S::G_INH ] / node.P_.tau_rise_in ); + + return GSL_SUCCESS; +} + +/* ---------------------------------------------------------------- + * Default constructors defining default parameters and state + * ---------------------------------------------------------------- */ + +nest::hh_cond_beta_gap_traub::Parameters_::Parameters_() + : g_Na( 20000.0 ) // Sodium Conductance (nS) + , g_K( 6000.0 ) // Potassium Conductance (nS) + , g_L( 10.0 ) // Leak Conductance (nS) + , C_m( 200.0 ) // Membrane Capacitance (pF) + , E_Na( 50.0 ) // Sodium Reversal potential (mV) + , E_K( -90.0 ) // Potassium Reversal potential (mV) + , E_L( -60.0 ) // Leak Reversal potential (mV) + , V_T( -50.0 ) // adjusts firing threshold (mV) + , E_ex( 0.0 ) // Excitatory reversal potential (mV) + , E_in( -80.0 ) // Inhibitory reversal potential (mV) + , tau_rise_ex( 0.5 ) // Excitatory Synaptic Rise Time Constant (ms) + , tau_decay_ex( 5.0 ) // Excitatory Synaptic Decay Time Constant (ms) + , tau_rise_in( 0.5 ) // Inhibitory Synaptic Rise Time Constant (ms) + , tau_decay_in( 10.0 ) // Inhibitory Synaptic Decay Time Constant (ms) + , t_ref_( 2.0 ) // Refractory time in ms (ms) + , I_e( 0.0 ) // Stimulus Current (pA) +{ +} + +nest::hh_cond_beta_gap_traub::State_::State_( const Parameters_& p ) + : r_( 0 ) +{ + y_[ 0 ] = p.E_L; + for ( size_t i = 1; i < STATE_VEC_SIZE; ++i ) + { + y_[ i ] = 0.0; + } + + // equilibrium values for (in)activation variables + const double alpha_n = + 0.032 * ( 15. - y_[ 0 ] ) / ( std::exp( ( 15. - y_[ 0 ] ) / 5. ) - 1. ); + const double beta_n = 0.5 * std::exp( ( 10. - y_[ 0 ] ) / 40. ); + const double alpha_m = + 0.32 * ( 13. - y_[ 0 ] ) / ( std::exp( ( 13. - y_[ 0 ] ) / 4. ) - 1. ); + const double beta_m = + 0.28 * ( y_[ 0 ] - 40. ) / ( std::exp( ( y_[ 0 ] - 40. ) / 5. ) - 1. ); + const double alpha_h = 0.128 * std::exp( ( 17. - y_[ 0 ] ) / 18. ); + const double beta_h = 4. / ( 1. + std::exp( ( 40. - y_[ 0 ] ) / 5. ) ); + + y_[ HH_H ] = alpha_h / ( alpha_h + beta_h ); + y_[ HH_N ] = alpha_n / ( alpha_n + beta_n ); + y_[ HH_M ] = alpha_m / ( alpha_m + beta_m ); +} + +nest::hh_cond_beta_gap_traub::State_::State_( const State_& s ) + : r_( s.r_ ) +{ + for ( size_t i = 0; i < STATE_VEC_SIZE; ++i ) + { + y_[ i ] = s.y_[ i ]; + } +} + +nest::hh_cond_beta_gap_traub::State_& nest::hh_cond_beta_gap_traub::State_:: +operator=( const State_& s ) +{ + assert( this != &s ); // would be bad logical error in program + for ( size_t i = 0; i < STATE_VEC_SIZE; ++i ) + { + y_[ i ] = s.y_[ i ]; + } + r_ = s.r_; + return *this; +} + +/* ---------------------------------------------------------------- + * Parameter and state extractions and manipulation functions + * ---------------------------------------------------------------- */ + +void +nest::hh_cond_beta_gap_traub::Parameters_::get( DictionaryDatum& d ) const +{ + def< double >( d, names::g_Na, g_Na ); + def< double >( d, names::g_K, g_K ); + def< double >( d, names::g_L, g_L ); + def< double >( d, names::C_m, C_m ); + def< double >( d, names::E_Na, E_Na ); + def< double >( d, names::E_K, E_K ); + def< double >( d, names::E_L, E_L ); + def< double >( d, names::V_T, V_T ); + def< double >( d, names::E_ex, E_ex ); + def< double >( d, names::E_in, E_in ); + def< double >( d, names::tau_rise_ex, tau_rise_ex ); + def< double >( d, names::tau_decay_ex, tau_decay_ex ); + def< double >( d, names::tau_rise_in, tau_rise_in ); + def< double >( d, names::tau_decay_in, tau_decay_in ); + def< double >( d, names::t_ref, t_ref_ ); + def< double >( d, names::I_e, I_e ); +} + +void +nest::hh_cond_beta_gap_traub::Parameters_::set( const DictionaryDatum& d ) +{ + updateValue< double >( d, names::g_Na, g_Na ); + updateValue< double >( d, names::g_K, g_K ); + updateValue< double >( d, names::g_L, g_L ); + updateValue< double >( d, names::C_m, C_m ); + updateValue< double >( d, names::E_Na, E_Na ); + updateValue< double >( d, names::E_K, E_K ); + updateValue< double >( d, names::E_L, E_L ); + updateValue< double >( d, names::V_T, V_T ); + updateValue< double >( d, names::E_ex, E_ex ); + updateValue< double >( d, names::E_in, E_in ); + updateValue< double >( d, names::tau_rise_ex, tau_rise_ex ); + updateValue< double >( d, names::tau_decay_ex, tau_decay_ex ); + updateValue< double >( d, names::tau_rise_in, tau_rise_in ); + updateValue< double >( d, names::tau_decay_in, tau_decay_in ); + updateValue< double >( d, names::t_ref, t_ref_ ); + updateValue< double >( d, names::I_e, I_e ); + + if ( C_m <= 0 ) + { + throw BadProperty( "Capacitance must be strictly positive." ); + } + + if ( t_ref_ < 0 ) + { + throw BadProperty( "Refractory time cannot be negative." ); + } + + if ( tau_rise_ex <= 0 || tau_decay_ex <= 0 || tau_rise_in <= 0 + || tau_decay_in <= 0 ) + { + throw BadProperty( "All time constants must be strictly positive." ); + } + + if ( g_K < 0 || g_Na < 0 || g_L < 0 ) + { + throw BadProperty( "All conductances must be non-negative." ); + } +} + +void +nest::hh_cond_beta_gap_traub::State_::get( DictionaryDatum& d ) const +{ + def< double >( d, names::V_m, y_[ V_M ] ); // Membrane potential + def< double >( d, names::Act_m, y_[ HH_M ] ); + def< double >( d, names::Act_h, y_[ HH_H ] ); + def< double >( d, names::Inact_n, y_[ HH_N ] ); +} + +void +nest::hh_cond_beta_gap_traub::State_::set( const DictionaryDatum& d, + const Parameters_& ) +{ + updateValue< double >( d, names::V_m, y_[ V_M ] ); + updateValue< double >( d, names::Act_m, y_[ HH_M ] ); + updateValue< double >( d, names::Act_h, y_[ HH_H ] ); + updateValue< double >( d, names::Inact_n, y_[ HH_N ] ); + if ( y_[ HH_M ] < 0 || y_[ HH_H ] < 0 || y_[ HH_N ] < 0 ) + { + throw BadProperty( "All (in)activation variables must be non-negative." ); + } +} + +nest::hh_cond_beta_gap_traub::Buffers_::Buffers_( hh_cond_beta_gap_traub& n ) + : logger_( n ) + , s_( 0 ) + , c_( 0 ) + , e_( 0 ) +{ + // Initialization of the remaining members is deferred to + // init_buffers_(). +} + +nest::hh_cond_beta_gap_traub::Buffers_::Buffers_( const Buffers_&, + hh_cond_beta_gap_traub& n ) + : logger_( n ) + , s_( 0 ) + , c_( 0 ) + , e_( 0 ) +{ + // Initialization of the remaining members is deferred to + // init_buffers_(). +} + +/* ---------------------------------------------------------------- + * Default and copy constructor for node, and destructor + * ---------------------------------------------------------------- */ + +nest::hh_cond_beta_gap_traub::hh_cond_beta_gap_traub() + : Archiving_Node() + , P_() + , S_( P_ ) + , B_( *this ) +{ + recordablesMap_.create(); + Node::set_node_uses_wfr( kernel().simulation_manager.use_wfr() ); +} + +nest::hh_cond_beta_gap_traub::hh_cond_beta_gap_traub( + const hh_cond_beta_gap_traub& n ) + : Archiving_Node( n ) + , P_( n.P_ ) + , S_( n.S_ ) + , B_( n.B_, *this ) +{ + Node::set_node_uses_wfr( kernel().simulation_manager.use_wfr() ); +} + +nest::hh_cond_beta_gap_traub::~hh_cond_beta_gap_traub() +{ + // GSL structs may not have been allocated, so we need to protect destruction + if ( B_.s_ ) + { + gsl_odeiv_step_free( B_.s_ ); + } + if ( B_.c_ ) + { + gsl_odeiv_control_free( B_.c_ ); + } + if ( B_.e_ ) + { + gsl_odeiv_evolve_free( B_.e_ ); + } +} + +/* ---------------------------------------------------------------- + * Node initialization functions + * ---------------------------------------------------------------- */ + +void +nest::hh_cond_beta_gap_traub::init_state_( const Node& proto ) +{ + const hh_cond_beta_gap_traub& pr = + downcast< hh_cond_beta_gap_traub >( proto ); + S_ = pr.S_; +} + +void +nest::hh_cond_beta_gap_traub::init_buffers_() +{ + B_.spike_exc_.clear(); // includes resize + B_.spike_inh_.clear(); // includes resize + B_.currents_.clear(); // includes resize + + // allocate strucure for gap events here + // function is called from Scheduler::prepare_nodes() before the + // first call to update + // so we already know which interpolation scheme to use according + // to the properties of this neurons + // determine size of structure depending on interpolation scheme + // and unsigned int Scheduler::min_delay() (number of simulation time steps + // per min_delay step) + + // resize interpolation_coefficients depending on interpolation order + const size_t buffer_size = kernel().connection_manager.get_min_delay() + * ( kernel().simulation_manager.get_wfr_interpolation_order() + 1 ); + + B_.interpolation_coefficients.resize( buffer_size, 0.0 ); + + B_.last_y_values.resize( kernel().connection_manager.get_min_delay(), 0.0 ); + + B_.sumj_g_ij_ = 0.0; + + Archiving_Node::clear_history(); + + B_.logger_.reset(); + + B_.step_ = Time::get_resolution().get_ms(); + B_.IntegrationStep_ = B_.step_; + + if ( B_.s_ == 0 ) + { + B_.s_ = + gsl_odeiv_step_alloc( gsl_odeiv_step_rkf45, State_::STATE_VEC_SIZE ); + } + else + { + gsl_odeiv_step_reset( B_.s_ ); + } + + if ( B_.c_ == 0 ) + { + B_.c_ = gsl_odeiv_control_y_new( 1e-3, 0.0 ); + } + else + { + gsl_odeiv_control_init( B_.c_, 1e-3, 0.0, 1.0, 0.0 ); + } + + if ( B_.e_ == 0 ) + { + B_.e_ = gsl_odeiv_evolve_alloc( State_::STATE_VEC_SIZE ); + } + else + { + gsl_odeiv_evolve_reset( B_.e_ ); + } + + B_.sys_.function = hh_cond_beta_gap_traub_dynamics; + B_.sys_.jacobian = 0; + B_.sys_.dimension = State_::STATE_VEC_SIZE; + B_.sys_.params = reinterpret_cast< void* >( this ); + + B_.I_stim_ = 0.0; +} + +double +nest::hh_cond_beta_gap_traub::get_normalisation_factor( double tau_rise, + double tau_decay ) +{ + // Factor used to normalise the synaptic conductance such that + // incoming spike causes a peak conductance of 1 nS. + // The denominator (denom1) that appears in the expression of the peak time + // is computed here to check that it is != 0 + // another denominator denom2 appears in the expression of the + // normalization factor g0 + // Both denom1 and denom2 are null if tau_decay = tau_rise, but they + // can also be null if tau_decay and tau_rise are not equal but very + // close to each other, due to the numerical precision limits. + // In such case the beta function reduces to the alpha function, + // and the normalization factor for the alpha function should be used. + double denom1 = tau_decay - tau_rise; + double normalisation_factor = 0; + if ( std::abs( denom1 ) > std::numeric_limits< double >::epsilon() ) + // if rise time != decay time use beta function + { + // peak time + const double t_p = + tau_decay * tau_rise * std::log( tau_decay / tau_rise ) / denom1; + // another denominator is computed here to check that it is != 0 + double denom2 = std::exp( -t_p / tau_decay ) - std::exp( -t_p / tau_rise ); + normalisation_factor = ( 1. / tau_rise - 1. / tau_decay ) / denom2; + } + else // if rise time == decay time use alpha function + { + normalisation_factor = 1. * numerics::e / tau_decay; + } + return normalisation_factor; +} + +void +nest::hh_cond_beta_gap_traub::calibrate() +{ + // ensures initialization in case mm connected after Simulate + B_.logger_.init(); + + V_.PSConInit_E = nest::hh_cond_beta_gap_traub::get_normalisation_factor( + P_.tau_rise_ex, P_.tau_decay_ex ); + V_.PSConInit_I = nest::hh_cond_beta_gap_traub::get_normalisation_factor( + P_.tau_rise_in, P_.tau_decay_in ); + + V_.refractory_counts_ = Time( Time::ms( P_.t_ref_ ) ).get_steps(); + V_.U_old_ = S_.y_[ State_::V_M ]; + + // since t_ref_ >= 0, this can only fail in error + assert( V_.refractory_counts_ >= 0 ); +} + +/* ---------------------------------------------------------------- + * Update and spike handling functions + * ---------------------------------------------------------------- */ + +bool +nest::hh_cond_beta_gap_traub::update_( Time const& origin, + const long from, + const long to, + const bool called_from_wfr_update ) +{ + + assert( + to >= 0 && ( delay ) from < kernel().connection_manager.get_min_delay() ); + assert( from < to ); + + const size_t interpolation_order = + kernel().simulation_manager.get_wfr_interpolation_order(); + const double wfr_tol = kernel().simulation_manager.get_wfr_tol(); + bool wfr_tol_exceeded = false; + + // allocate memory to store the new interpolation coefficients + // to be sent by gap event + const size_t buffer_size = + kernel().connection_manager.get_min_delay() * ( interpolation_order + 1 ); + std::vector< double > new_coefficients( buffer_size, 0.0 ); + + // parameters needed for piecewise interpolation + double y_i = 0.0, y_ip1 = 0.0, hf_i = 0.0, hf_ip1 = 0.0; + double f_temp[ State_::STATE_VEC_SIZE ]; + + for ( long lag = from; lag < to; ++lag ) + { + + // B_.lag is needed by hh_cond_beta_gap_traub_dynamics to + // determine the current section + B_.lag_ = lag; + + if ( called_from_wfr_update ) + { + y_i = S_.y_[ State_::V_M ]; + if ( interpolation_order == 3 ) + { + hh_cond_beta_gap_traub_dynamics( + 0, S_.y_, f_temp, reinterpret_cast< void* >( this ) ); + hf_i = B_.step_ * f_temp[ State_::V_M ]; + } + } + + double t = 0.0; + const double U_old = S_.y_[ State_::V_M ]; + + // numerical integration with adaptive step size control: + // ------------------------------------------------------ + // gsl_odeiv_evolve_apply performs only a single numerical + // integration step, starting from t and bounded by step; + // the while-loop ensures integration over the whole simulation + // step (0, step] if more than one integration step is needed due + // to a small integration step size; + // note that (t+IntegrationStep > step) leads to integration over + // (t, step] and afterwards setting t to step, but it does not + // enforce setting IntegrationStep to step-t; this is of advantage + // for a consistent and efficient integration across subsequent + // simulation intervals + while ( t < B_.step_ ) + { + const int status = gsl_odeiv_evolve_apply( B_.e_, + B_.c_, + B_.s_, + &B_.sys_, // system of ODE + &t, // from t + B_.step_, // to t <= step + &B_.IntegrationStep_, // integration step size + S_.y_ ); // neuronal state + if ( status != GSL_SUCCESS ) + { + throw GSLSolverFailure( get_name(), status ); + } + } + + if ( not called_from_wfr_update ) + { + S_.y_[ State_::DG_EXC ] += + B_.spike_exc_.get_value( lag ) * V_.PSConInit_E; + S_.y_[ State_::DG_INH ] += + B_.spike_inh_.get_value( lag ) * V_.PSConInit_I; + // sending spikes: crossing 0 mV, pseudo-refractoriness and local + // maximum... + // refractory? + if ( S_.r_ > 0 ) + { + --S_.r_; + } + else + // ( threshold && maximum ) + if ( S_.y_[ State_::V_M ] >= P_.V_T + 30. + && U_old > S_.y_[ State_::V_M ] ) + { + S_.r_ = V_.refractory_counts_; + + set_spiketime( Time::step( origin.get_steps() + lag + 1 ) ); + + SpikeEvent se; + kernel().event_delivery_manager.send( *this, se, lag ); + } + + // log state data + B_.logger_.record_data( origin.get_steps() + lag ); + + // set new input current + B_.I_stim_ = B_.currents_.get_value( lag ); + } + else // if(called_from_wfr_update) + { + S_.y_[ State_::DG_EXC ] += + B_.spike_exc_.get_value_wfr_update( lag ) * V_.PSConInit_E; + S_.y_[ State_::DG_INH ] += + B_.spike_inh_.get_value_wfr_update( lag ) * V_.PSConInit_I; + // check if deviation from last iteration exceeds wfr_tol + wfr_tol_exceeded = wfr_tol_exceeded + or fabs( S_.y_[ State_::V_M ] - B_.last_y_values[ lag ] ) > wfr_tol; + B_.last_y_values[ lag ] = S_.y_[ State_::V_M ]; + + // update different interpolations + + // constant term is the same for each interpolation order + new_coefficients[ lag * ( interpolation_order + 1 ) + 0 ] = y_i; + + switch ( interpolation_order ) + { + case 0: + break; + + case 1: + y_ip1 = S_.y_[ State_::V_M ]; + + new_coefficients[ lag * ( interpolation_order + 1 ) + 1 ] = y_ip1 - y_i; + break; + + case 3: + y_ip1 = S_.y_[ State_::V_M ]; + hh_cond_beta_gap_traub_dynamics( + B_.step_, S_.y_, f_temp, reinterpret_cast< void* >( this ) ); + hf_ip1 = B_.step_ * f_temp[ State_::V_M ]; + + new_coefficients[ lag * ( interpolation_order + 1 ) + 1 ] = hf_i; + new_coefficients[ lag * ( interpolation_order + 1 ) + 2 ] = + -3 * y_i + 3 * y_ip1 - 2 * hf_i - hf_ip1; + new_coefficients[ lag * ( interpolation_order + 1 ) + 3 ] = + 2 * y_i - 2 * y_ip1 + hf_i + hf_ip1; + break; + + default: + throw BadProperty( "Interpolation order must be 0, 1, or 3." ); + } + } + + + } // end for-loop + + // if not called_from_wfr_update perform constant extrapolation + // and reset last_y_values + if ( not called_from_wfr_update ) + { + for ( long temp = from; temp < to; ++temp ) + { + new_coefficients[ temp * ( interpolation_order + 1 ) + 0 ] = + S_.y_[ State_::V_M ]; + } + + std::vector< double >( kernel().connection_manager.get_min_delay(), 0.0 ) + .swap( B_.last_y_values ); + } + + // Send gap-event + GapJunctionEvent ge; + ge.set_coeffarray( new_coefficients ); + kernel().event_delivery_manager.send_secondary( *this, ge ); + + // Reset variables + B_.sumj_g_ij_ = 0.0; + std::vector< double >( buffer_size, 0.0 ) + .swap( B_.interpolation_coefficients ); + + return wfr_tol_exceeded; +} + +void +nest::hh_cond_beta_gap_traub::handle( SpikeEvent& e ) +{ + assert( e.get_delay_steps() > 0 ); + + if ( e.get_weight() > 0.0 ) + { + B_.spike_exc_.add_value( e.get_rel_delivery_steps( + kernel().simulation_manager.get_slice_origin() ), + e.get_weight() * e.get_multiplicity() ); + } + else + { + // add with negative weight, ie positive value, since we are changing a + // conductance + B_.spike_inh_.add_value( e.get_rel_delivery_steps( + kernel().simulation_manager.get_slice_origin() ), + -e.get_weight() * e.get_multiplicity() ); + } +} + +void +nest::hh_cond_beta_gap_traub::handle( CurrentEvent& e ) +{ + assert( e.get_delay_steps() > 0 ); + + const double c = e.get_current(); + const double w = e.get_weight(); + + // add weighted current; HEP 2002-10-04 + B_.currents_.add_value( + e.get_rel_delivery_steps( kernel().simulation_manager.get_slice_origin() ), + w * c ); +} + +void +nest::hh_cond_beta_gap_traub::handle( DataLoggingRequest& e ) +{ + B_.logger_.handle( e ); +} + +void +nest::hh_cond_beta_gap_traub::handle( GapJunctionEvent& e ) +{ + const double weight = e.get_weight(); + + B_.sumj_g_ij_ += weight; + + size_t i = 0; + std::vector< unsigned int >::iterator it = e.begin(); + // The call to get_coeffvalue( it ) in this loop also advances the iterator it + while ( it != e.end() ) + { + B_.interpolation_coefficients[ i ] += weight * e.get_coeffvalue( it ); + ++i; + } +} + +} // namespace nest + +#endif // HAVE_GSL diff --git a/models/hh_cond_beta_gap_traub.h b/models/hh_cond_beta_gap_traub.h new file mode 100644 index 0000000000..56ffb7744e --- /dev/null +++ b/models/hh_cond_beta_gap_traub.h @@ -0,0 +1,501 @@ +/* + * hh_cond_beta_gap_traub.h + * + * 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 . + * + */ + +#ifndef HH_COND_BETA_GAP_TRAUB_H +#define HH_COND_BETA_GAP_TRAUB_H + +// Generated includes: +#include "config.h" + +#ifdef HAVE_GSL + +// C includes: +#include +#include +#include +#include + +// Includes from nestkernel: +#include "archiving_node.h" +#include "connection.h" +#include "event.h" +#include "nest_types.h" +#include "node.h" +#include "recordables_map.h" +#include "ring_buffer.h" +#include "universal_data_logger.h" + +namespace nest +{ + +/** + * Function computing right-hand side of ODE for GSL solver. + * @note Must be declared here so we can befriend it in class. + * @note Must have C-linkage for passing to GSL. Internally, it is + * a first-class C++ function, but cannot be a member function + * because of the C-linkage. + * @note No point in declaring it inline, since it is called + * through a function pointer. + * @param void* Pointer to model neuron instance. + */ +extern "C" int +hh_cond_beta_gap_traub_dynamics( double, const double*, double*, void* ); + +/** @BeginDocumentation +Name: hh_cond_beta_gap_traub - modified Hodgkin-Huxley neuron as featured in +Brette et al (2007) review with added gap junction support and beta function +synaptic conductance. + +Description: + +hh_cond_beta_gap_traub is an implementation of a modified Hodgkin-Huxley model +that also supports gap junctions. + +This model was specifically developed for a major review of simulators [1], +based on a model of hippocampal pyramidal cells by Traub and Miles[2]. +The key differences between the current model and the model in [2] are: + +- This model is a point neuron, not a compartmental model. +- This model includes only I_Na and I_K, with simpler I_K dynamics than + in [2], so it has only three instead of eight gating variables; + in particular, all Ca dynamics have been removed. +- Incoming spikes induce an instantaneous conductance change followed by + exponential decay instead of activation over time. + +This model is primarily provided as reference implementation for hh_coba +example of the Brette et al (2007) review. Default parameter values are chosen +to match those used with NEST 1.9.10 when preparing data for [1]. Code for all +simulators covered is available from ModelDB [3]. + +Note: +In this model, a spike is emitted if + + V_m >= V_T + 30 mV and V_m has fallen during the current time step + +To avoid that this leads to multiple spikes during the falling flank of a +spike, it is essential to chose a sufficiently long refractory period. +Traub and Miles used t_ref = 3 ms [2, p 118], while we used t_ref = 2 ms +in [2]. + +Post-synaptic currents +Incoming spike events induce a post-synaptic change of conductance modelled by a +beta function as outlined in [4,5]. The beta function is normalised such that an +event of weight 1.0 results in a peak current of 1 nS at t = tau_rise_xx where +xx is ex or in. + +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. + +Gap Junctions +Gap Junctions are implemented by a gap current of the form g_ij( V_i - V_j). + +Parameters: + +The following parameters can be set in the status dictionary. + +V_m double - Membrane potential in mV +V_T double - Voltage offset that controls dynamics. For default + parameters, V_T = -63mV results in a threshold around + -50mV. +E_L double - Leak reversal potential in mV. +C_m double - Capacity of the membrane in pF. +g_L double - Leak conductance in nS. +tau_rise_ex double - Excitatory synaptic beta function rise time in ms. +tau_decay_ex double - Excitatory synaptic beta function decay time in ms. +tau_rise_in double - Inhibitory synaptic beta function rise time in ms. +tau_decay_in double - Inhibitory synaptic beta function decay time in ms. +t_ref double - Duration of refractory period in ms (see Note). +E_ex double - Excitatory synaptic reversal potential in mV. +E_in double - Inhibitory synaptic reversal potential in mV. +E_Na double - Sodium reversal potential in mV. +g_Na double - Sodium peak conductance in nS. +E_K double - Potassium reversal potential in mV. +g_K double - Potassium peak conductance in nS. +I_e double - External input current in pA. + +References: + +[1] Brette R et al (2007) Simulation of networks of spiking neurons: A review + of tools and strategies. J Comp Neurosci 23:349-98. + doi 10.1007/s10827-007-0038-6 +[2] Traub RD and Miles R (1991) Neuronal Networks of the Hippocampus. + Cambridge University Press, Cambridge UK. +[3] http://modeldb.yale.edu/83319 +[4] Rotter & Diesmann, Biol Cybern 81:381 (1999) +[5] Roth and van Rossum, + Ch 6, in De Schutter, Computational Modeling Methods for Neuroscientists, + MIT Press, 2010. + +Sends: SpikeEvent + +Receives: SpikeEvent, CurrentEvent, DataLoggingRequest + +Author: Daniel Naoumenko (modified hh_cond_exp_traub by Schrader and +hh_psc_alpha_gap by Jan Hahne, Moritz Helias and Susanne Kunkel) + +SeeAlso: hh_psc_alpha_gap, hh_cond_exp_traub, gap_junction, iaf_cond_beta +*/ +class hh_cond_beta_gap_traub : public Archiving_Node +{ + +public: + typedef Node base; + + hh_cond_beta_gap_traub(); + hh_cond_beta_gap_traub( const hh_cond_beta_gap_traub& ); + ~hh_cond_beta_gap_traub(); + + /** + * Import sets of overloaded virtual functions. + * @see Technical Issues / Virtual Functions: Overriding, Overloading, and + * Hiding + */ + using Node::handle; + using Node::handles_test_event; + using Node::sends_secondary_event; + + port send_test_event( Node& target, rport receptor_type, synindex, bool ); + + void handle( SpikeEvent& ); + void handle( CurrentEvent& ); + void handle( DataLoggingRequest& ); + void handle( GapJunctionEvent& ); + + port handles_test_event( SpikeEvent&, rport ); + port handles_test_event( CurrentEvent&, rport ); + port handles_test_event( DataLoggingRequest&, rport ); + port handles_test_event( GapJunctionEvent&, rport ); + + void + sends_secondary_event( GapJunctionEvent& ) + { + } + + void get_status( DictionaryDatum& ) const; + void set_status( const DictionaryDatum& ); + +private: + void init_state_( const Node& proto ); + void init_buffers_(); + double get_normalisation_factor( double, double ); + void calibrate(); + + /** This is the actual update function. The additional boolean parameter + * determines if the function is called by update (false) or wfr_update (true) + */ + bool update_( Time const&, const long, const long, const bool ); + + void update( Time const&, const long, const long ); + bool wfr_update( Time const&, const long, const long ); + + // END Boilerplate function declarations ---------------------------- + + // Friends -------------------------------------------------------- + + // make dynamics function quasi-member + friend int + hh_cond_beta_gap_traub_dynamics( double, const double*, double*, void* ); + + // The next two classes need to be friends to access the State_ class/member + friend class RecordablesMap< hh_cond_beta_gap_traub >; + friend class UniversalDataLogger< hh_cond_beta_gap_traub >; + +private: + // ---------------------------------------------------------------- + + /** + * Independent parameters of the model. + */ + struct Parameters_ + { + double g_Na; //!< Sodium Conductance in nS + double g_K; //!< Potassium Conductance in nS + double g_L; //!< Leak Conductance in nS + double C_m; //!< Membrane Capacitance in pF + double E_Na; //!< Sodium Reversal Potential in mV + double E_K; //!< Potassium Reversal Potential in mV + double E_L; //!< Leak Reversal Potential in mV + double V_T; //!< Voltage offset for dynamics in mV + double E_ex; //!< Excitatory reversal Potential in mV + double E_in; //!< Inhibitory reversal Potential in mV + double tau_rise_ex; //!< Excitatory Synaptic Rise Time Constant in ms + double tau_decay_ex; //!< Excitatory Synaptic Decay Time Constant in ms + double tau_rise_in; //!< Inhibitory Synaptic Rise Time Constant in ms + double tau_decay_in; //!< Inhibitory Synaptic Decay Time Constant in ms + double t_ref_; //!< Refractory time in ms + double I_e; //!< External Current in pA + + Parameters_(); + + void get( DictionaryDatum& ) const; //!< Store current values in dictionary + void set( const DictionaryDatum& ); //!< Set values from dicitonary + }; + +public: + // ---------------------------------------------------------------- + + /** + * State variables of the model. + */ + struct State_ + { + + //! Symbolic indices to the elements of the state vector y + enum StateVecElems + { + V_M = 0, + HH_M, // 1 + HH_H, // 2 + HH_N, // 3 + DG_EXC, // 4 + G_EXC, // 5 + DG_INH, // 6 + G_INH, // 7 + STATE_VEC_SIZE + }; + + //! neuron state, must be C-array for GSL solver + double y_[ STATE_VEC_SIZE ]; + int r_; //!< number of refractory steps remaining + + State_( const Parameters_& p ); + State_( const State_& s ); + + State_& operator=( const State_& s ); + + void get( DictionaryDatum& ) const; + void set( const DictionaryDatum&, const Parameters_& ); + }; + + // Variables class ------------------------------------------------------- + + /** + * Internal variables of the model. + * Variables are re-initialized upon each call to Simulate. + */ + struct Variables_ + { + /** + * Impulse to add to DG_EXC on spike arrival to evoke unit-amplitude + * conductance excursion. + */ + double PSConInit_E; + + /** + * Impulse to add to DG_INH on spike arrival to evoke unit-amplitude + * conductance excursion. + */ + double PSConInit_I; + + //! refractory time in steps + int refractory_counts_; + double U_old_; // for spike-detection + }; + + // ---------------------------------------------------------------- + + /** + * Buffers of the model. + */ + struct Buffers_ + { + Buffers_( hh_cond_beta_gap_traub& ); //!< Sets buffer pointers to 0 + //! Sets buffer pointers to 0 + Buffers_( const Buffers_&, hh_cond_beta_gap_traub& ); + + //! Logger for all analog data + UniversalDataLogger< hh_cond_beta_gap_traub > logger_; + + /** buffers and sums up incoming spikes/currents */ + RingBuffer spike_exc_; + RingBuffer spike_inh_; + RingBuffer currents_; + + /** GSL ODE stuff */ + gsl_odeiv_step* s_; //!< stepping function + gsl_odeiv_control* c_; //!< adaptive stepsize control function + gsl_odeiv_evolve* e_; //!< evolution function + gsl_odeiv_system sys_; //!< struct describing system + + // IntergrationStep_ should be reset with the neuron on ResetNetwork, + // but remain unchanged during calibration. Since it is initialized with + // step_, and the resolution cannot change after nodes have been created, + // it is safe to place both here. + double step_; //!< step size in ms + double IntegrationStep_; //!< current integration time step, updated by GSL + + // remembers current lag for piecewise interpolation + long lag_; + + // remembers y_values from last wfr_update + std::vector< double > last_y_values; + + // summarized gap weight + double sumj_g_ij_; + + // summarized coefficients of the interpolation polynomial + std::vector< double > interpolation_coefficients; + + /** + * Input current injected by CurrentEvent. + * This variable is used to transport the current applied into the + * _dynamics function computing the derivative of the state vector. + * It must be a part of Buffers_, since it is initialized once before + * the first simulation, but not modified before later Simulate calls. + */ + double I_stim_; + }; + + // Access functions for UniversalDataLogger ------------------------------- + + //! Read out state vector elements, used by UniversalDataLogger + template < State_::StateVecElems elem > + double + get_y_elem_() const + { + return S_.y_[ elem ]; + } + + Parameters_ P_; + State_ S_; + Variables_ V_; + Buffers_ B_; + + //! Mapping of recordables names to access functions + static RecordablesMap< hh_cond_beta_gap_traub > recordablesMap_; +}; + +inline void +hh_cond_beta_gap_traub::update( Time const& origin, + const long from, + const long to ) +{ + update_( origin, from, to, false ); +} + +inline bool +hh_cond_beta_gap_traub::wfr_update( Time const& origin, + const long from, + const long to ) +{ + State_ old_state = S_; // save state before wfr_update + const bool wfr_tol_exceeded = update_( origin, from, to, true ); + S_ = old_state; // restore old state + + return not wfr_tol_exceeded; +} + +inline port +hh_cond_beta_gap_traub::send_test_event( Node& target, + rport receptor_type, + synindex, + bool ) +{ + SpikeEvent e; + e.set_sender( *this ); + + return target.handles_test_event( e, receptor_type ); +} + + +inline port +hh_cond_beta_gap_traub::handles_test_event( SpikeEvent&, rport receptor_type ) +{ + if ( receptor_type != 0 ) + { + throw UnknownReceptorType( receptor_type, get_name() ); + } + return 0; +} + +inline port +hh_cond_beta_gap_traub::handles_test_event( CurrentEvent&, rport receptor_type ) +{ + if ( receptor_type != 0 ) + { + throw UnknownReceptorType( receptor_type, get_name() ); + } + return 0; +} + +inline port +hh_cond_beta_gap_traub::handles_test_event( DataLoggingRequest& dlr, + rport receptor_type ) +{ + if ( receptor_type != 0 ) + { + throw UnknownReceptorType( receptor_type, get_name() ); + } + return B_.logger_.connect_logging_device( dlr, recordablesMap_ ); +} + +inline port +hh_cond_beta_gap_traub::handles_test_event( GapJunctionEvent&, + rport receptor_type ) +{ + if ( receptor_type != 0 ) + { + throw UnknownReceptorType( receptor_type, get_name() ); + } + return 0; +} + +inline void +hh_cond_beta_gap_traub::get_status( DictionaryDatum& d ) const +{ + P_.get( d ); + S_.get( d ); + Archiving_Node::get_status( d ); + + ( *d )[ names::recordables ] = recordablesMap_.get_list(); + + def< double >( d, names::t_spike, get_spiketime_ms() ); +} + +inline void +hh_cond_beta_gap_traub::set_status( const DictionaryDatum& d ) +{ + Parameters_ ptmp = P_; // temporary copy in case of errors + ptmp.set( d ); // throws if BadProperty + State_ stmp = S_; // temporary copy in case of errors + stmp.set( d, ptmp ); // throws if BadProperty + + // We now know that (ptmp, stmp) are consistent. We do not + // write them back to (P_, S_) before we are also sure that + // the properties to be set in the parent class are internally + // consistent. + Archiving_Node::set_status( d ); + + // if we get here, temporaries contain consistent set of properties + P_ = ptmp; + S_ = stmp; + + calibrate(); +} + +} // namespace + + +#endif // HAVE_GSL +#endif // HH_COND_BETA_GAP_TRAUB_H diff --git a/models/modelsmodule.cpp b/models/modelsmodule.cpp index c34f011f6c..2887affd82 100644 --- a/models/modelsmodule.cpp +++ b/models/modelsmodule.cpp @@ -57,6 +57,7 @@ #include "gauss_rate.h" #include "ginzburg_neuron.h" #include "hh_cond_exp_traub.h" +#include "hh_cond_beta_gap_traub.h" #include "hh_psc_alpha.h" #include "hh_psc_alpha_clopath.h" #include "hh_psc_alpha_gap.h" @@ -390,6 +391,8 @@ ModelsModule::init( SLIInterpreter* ) "iaf_cond_exp_sfa_rr" ); kernel().model_manager.register_node_model< iaf_cond_alpha_mc >( "iaf_cond_alpha_mc" ); + kernel().model_manager.register_node_model< hh_cond_beta_gap_traub >( + "hh_cond_beta_gap_traub" ); kernel().model_manager.register_node_model< hh_psc_alpha >( "hh_psc_alpha" ); kernel().model_manager.register_node_model< hh_psc_alpha_clopath >( "hh_psc_alpha_clopath" ); diff --git a/pynest/nest/tests/test_refractory.py b/pynest/nest/tests/test_refractory.py index 62a3187dea..56caa8c9f4 100644 --- a/pynest/nest/tests/test_refractory.py +++ b/pynest/nest/tests/test_refractory.py @@ -45,6 +45,7 @@ * ``aeif_cond_alpha_RK5`` * ``gif_pop_psc_exp`` * ``hh_cond_exp_traub`` +* ``hh_cond_beta_gap_traub`` * ``hh_psc_alpha`` * ``hh_psc_alpha_gap`` * ``iaf_psc_exp_ps_lossless`` @@ -82,6 +83,7 @@ "aeif_cond_alpha_RK5", # This one is faulty and will be removed "gif_pop_psc_exp", # This one commits spikes at same time "hh_cond_exp_traub", # This one does not support V_reset + "hh_cond_beta_gap_traub", # This one does not support V_reset "hh_psc_alpha", # This one does not support V_reset "hh_psc_alpha_clopath", # This one does not support V_reset "hh_psc_alpha_gap", # This one does not support V_reset diff --git a/testsuite/regressiontests/ticket-421.sli b/testsuite/regressiontests/ticket-421.sli index b6547118cb..de156ba39a 100644 --- a/testsuite/regressiontests/ticket-421.sli +++ b/testsuite/regressiontests/ticket-421.sli @@ -51,7 +51,7 @@ Author: Hans Ekkehard Plesser, 2010-05-05 /exclude_models [/aeif_cond_exp /aeif_cond_alpha /aeif_cond_alpha_RK5 /a2eif_cond_exp /a2eif_cond_exp_HW /aeif_cond_alpha_multisynapse /aeif_psc_delta_clopath /aeif_psc_exp /aeif_psc_alpha /aeif_psc_delta /aeif_cond_beta_multisynapse /hh_cond_exp_traub - /hh_psc_alpha /hh_psc_alpha_clopath /hh_psc_alpha_gap /ht_neuron /ht_neuron_fs + /hh_cond_beta_gap_traub /hh_psc_alpha /hh_psc_alpha_clopath /hh_psc_alpha_gap /ht_neuron /ht_neuron_fs /iaf_cond_exp_sfa_rr /izhikevich] def % use power-of-two resolution to avoid roundof problems diff --git a/testsuite/unittests/test_hh_cond_beta_gap_traub.sli b/testsuite/unittests/test_hh_cond_beta_gap_traub.sli new file mode 100644 index 0000000000..829cfbaf98 --- /dev/null +++ b/testsuite/unittests/test_hh_cond_beta_gap_traub.sli @@ -0,0 +1,160 @@ +/* + * test_hh_cond_beta_gap_traub.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_hh_cond_beta_gap_traub - sli script for overall test of hh_cond_beta_gap_traub model + with gap_junction connection + + Synopsis: (test_hh_cond_beta_gap_traub) run -> compares response to current step with reference data + + Description: + test_hh_cond_beta_gap_traub.sli is an overall test of the hh_cond_beta_gap_traub model connected + by gap_junction connection. + + Two neurons of whom one receives an constant input current of 200 pA are connected + by gap_junction with an (unrealistic) high gap weight. The accurate functionality + of the gap_junction feature is tested by measuring the membrane potential of the + neuron without input current. + + Although 0.1 cannot be represented in the IEEE double data type, it + is safe to simulate with a resolution (computation step size) of 0.1 + ms because by default nest is built with a timebase enabling exact + representation of 0.1 ms. + + The expected output is documented at the end of + the script. The textual output of the voltmeter documented in this file + can be regenerated by setting adding /to_screen true to the SetStatus + call of vm below. + + Author: October 2015, Hahne + SeeAlso: testsuite::test_gap_junction, hh_cond_beta_gap_traub, gap_junction +*/ + +(unittest) run +/unittest using + +% The following test needs the model hh_cond_beta_gap_traub, so +% this test should only run if we have GSL +skip_if_without_gsl + +0.1 /h Set + +ResetKernel + +0 << + /local_num_threads 1 + /resolution h + /use_wfr true + /wfr_tol 0.0001 + /wfr_interpolation_order 3 + /wfr_max_iterations 10 + /wfr_comm_interval 1.0 + >> SetStatus + +/hh_cond_beta_gap_traub Create /neuron1 Set +/hh_cond_beta_gap_traub Create /neuron2 Set + +neuron1 +<< + /I_e 200. +>> SetStatus + +/voltmeter Create /vm Set +vm << /withtime true /time_in_steps true /interval h >> SetStatus + +[neuron1] [neuron2] +<< /rule /one_to_one /make_symmetric true >> +<< /model /gap_junction /weight 20.0 >> Connect + +vm neuron2 1.0 h Connect + + +20 Simulate + +{ % reference data + dup Transpose First /test_times Set % times of reference + + vm [/events [/times /V_m]] get cva % array of recorded data + 5 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: +% +% The output send to std::cout is a superposition of +% the output of the voltmeter. +% +% time (in steps) voltage (in mV) +[ +[ 1 -60 ] +[ 2 -5.999800e+01 ] +[ 3 -5.999600e+01 ] +[ 4 -5.999200e+01 ] +[ 5 -5.998800e+01 ] +[ 6 -5.998300e+01 ] +[ 7 -5.997700e+01 ] +[ 8 -5.997000e+01 ] +[ 9 -5.996300e+01 ] +[ 10 -5.995500e+01 ] +[ 11 -5.994600e+01 ] +[ 12 -5.993600e+01 ] +[ 13 -5.992600e+01 ] +[ 14 -5.991400e+01 ] +[ 15 -5.990300e+01 ] +% +% ... +% +[ 117 -5.746100e+01 ] +[ 118 -5.743600e+01 ] +[ 119 -5.741100e+01 ] +[ 120 -5.738600e+01 ] +[ 121 -5.736100e+01 ] +[ 122 -5.733600e+01 ] +[ 123 -5.731100e+01 ] +[ 124 -5.728600e+01 ] +[ 125 -5.726100e+01 ] +[ 126 -5.723700e+01 ] +[ 127 -5.721200e+01 ] +[ 128 -5.718800e+01 ] +[ 129 -5.716300e+01 ] +[ 130 -5.713900e+01 ] +[ 131 -5.711500e+01 ] +% +% ... +% +[ 190 -5.583700e+01 ] +[ 191 -5.581800e+01 ] +[ 192 -5.579900e+01 ] +[ 193 -5.578000e+01 ] +[ 194 -5.576100e+01 ] +[ 195 -5.574200e+01 ] +[ 196 -5.572400e+01 ] +[ 197 -5.570500e+01 ] +[ 198 -5.568600e+01 ] +[ 199 -5.566800e+01 ] +] + +exch assert_or_die \ No newline at end of file