1 /*
2 * aeif_cond_exp.cpp
3 *
4 * This file is part of NEST.
5 *
6 * Copyright (C) 2004 The NEST Initiative
7 *
8 * NEST is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 2 of the License, or
11 * (at your option) any later version.
12 *
13 * NEST is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with NEST. If not, see <http://www.gnu.org/licenses/>.
20 *
21 */
22
23 #include "aeif_cond_exp.h"
24
25 #ifdef HAVE_GSL
26
27 // C++ includes:
28 #include <cmath>
29 #include <cstdio>
30 #include <iomanip>
31 #include <iostream>
32 #include <limits>
33
34 // Includes from libnestutil:
35 #include "dict_util.h"
36 #include "numerics.h"
37
38 // Includes from nestkernel:
39 #include "exceptions.h"
40 #include "kernel_manager.h"
41 #include "nest_names.h"
42 #include "universal_data_logger_impl.h"
43
44 // Includes from sli:
45 #include "dict.h"
46 #include "dictutils.h"
47 #include "doubledatum.h"
48 #include "integerdatum.h"
49
50 /* ----------------------------------------------------------------
51 * Recordables map
52 * ---------------------------------------------------------------- */
53
54 nest::RecordablesMap< nest::aeif_cond_exp > nest::aeif_cond_exp::recordablesMap_;
55
56 namespace nest
57 {
58 /*
59 * template specialization must be placed in namespace
60 *
61 * Override the create() method with one call to RecordablesMap::insert_()
62 * for each quantity to be recorded.
63 */
64 template <>
65 void
create()66 RecordablesMap< aeif_cond_exp >::create()
67 {
68 // use standard names whereever you can for consistency!
69 insert_( names::V_m, &aeif_cond_exp::get_y_elem_< aeif_cond_exp::State_::V_M > );
70 insert_( names::g_ex, &aeif_cond_exp::get_y_elem_< aeif_cond_exp::State_::G_EXC > );
71 insert_( names::g_in, &aeif_cond_exp::get_y_elem_< aeif_cond_exp::State_::G_INH > );
72 insert_( names::w, &aeif_cond_exp::get_y_elem_< aeif_cond_exp::State_::W > );
73 }
74 }
75
76
77 extern "C" int
aeif_cond_exp_dynamics(double,const double y[],double f[],void * pnode)78 nest::aeif_cond_exp_dynamics( double, const double y[], double f[], void* pnode )
79 {
80 // a shorthand
81 typedef nest::aeif_cond_exp::State_ S;
82
83 // get access to node so we can almost work as in a member function
84 assert( pnode );
85 const nest::aeif_cond_exp& node = *( reinterpret_cast< nest::aeif_cond_exp* >( pnode ) );
86
87 const bool is_refractory = node.S_.r_ > 0;
88
89 // y[] here is---and must be---the state vector supplied by the integrator,
90 // not the state vector in the node, node.S_.y[].
91
92 // The following code is verbose for the sake of clarity. We assume that a
93 // good compiler will optimize the verbosity away ...
94
95 // Clamp membrane potential to V_reset while refractory, otherwise bound
96 // it to V_peak. Do not use V_.V_peak_ here, since that is set to V_th if
97 // Delta_T == 0.
98 const double& V = is_refractory ? node.P_.V_reset_ : std::min( y[ S::V_M ], node.P_.V_peak_ );
99 // shorthand for the other state variables
100 const double& g_ex = y[ S::G_EXC ];
101 const double& g_in = y[ S::G_INH ];
102 const double& w = y[ S::W ];
103
104 const double I_syn_exc = g_ex * ( V - node.P_.E_ex );
105 const double I_syn_inh = g_in * ( V - node.P_.E_in );
106
107 const double I_spike =
108 node.P_.Delta_T == 0. ? 0. : ( node.P_.g_L * node.P_.Delta_T * std::exp( ( V - node.P_.V_th ) / node.P_.Delta_T ) );
109
110 // dv/dt
111 f[ S::V_M ] = is_refractory ? 0. : ( -node.P_.g_L * ( V - node.P_.E_L ) + I_spike - I_syn_exc - I_syn_inh - w
112 + node.P_.I_e + node.B_.I_stim_ ) / node.P_.C_m;
113
114 f[ S::G_EXC ] = -g_ex / node.P_.tau_syn_ex; // Synaptic Conductance (nS)
115
116 f[ S::G_INH ] = -g_in / node.P_.tau_syn_in; // Synaptic Conductance (nS)
117
118 // Adaptation current w.
119 f[ S::W ] = ( node.P_.a * ( V - node.P_.E_L ) - w ) / node.P_.tau_w;
120
121 return GSL_SUCCESS;
122 }
123
124
125 /* ----------------------------------------------------------------
126 * Default constructors defining default parameters and state
127 * ---------------------------------------------------------------- */
128
Parameters_()129 nest::aeif_cond_exp::Parameters_::Parameters_()
130 : V_peak_( 0.0 ) // mV
131 , V_reset_( -60.0 ) // mV
132 , t_ref_( 0.0 ) // ms
133 , g_L( 30.0 ) // nS
134 , C_m( 281.0 ) // pF
135 , E_ex( 0.0 ) // mV
136 , E_in( -85.0 ) // mV
137 , E_L( -70.6 ) // mV
138 , Delta_T( 2.0 ) // mV
139 , tau_w( 144.0 ) // ms
140 , a( 4.0 ) // nS
141 , b( 80.5 ) // pA
142 , V_th( -50.4 ) // mV
143 , tau_syn_ex( 0.2 ) // ms
144 , tau_syn_in( 2.0 ) // ms
145 , I_e( 0.0 ) // pA
146 , gsl_error_tol( 1e-6 )
147 {
148 }
149
State_(const Parameters_ & p)150 nest::aeif_cond_exp::State_::State_( const Parameters_& p )
151 : r_( 0 )
152 {
153 y_[ 0 ] = p.E_L;
154 for ( size_t i = 1; i < STATE_VEC_SIZE; ++i )
155 {
156 y_[ i ] = 0;
157 }
158 }
159
State_(const State_ & s)160 nest::aeif_cond_exp::State_::State_( const State_& s )
161 : r_( s.r_ )
162 {
163 for ( size_t i = 0; i < STATE_VEC_SIZE; ++i )
164 {
165 y_[ i ] = s.y_[ i ];
166 }
167 }
168
operator =(const State_ & s)169 nest::aeif_cond_exp::State_& nest::aeif_cond_exp::State_::operator=( const State_& s )
170 {
171 r_ = s.r_;
172 for ( size_t i = 0; i < STATE_VEC_SIZE; ++i )
173 {
174 y_[ i ] = s.y_[ i ];
175 }
176 return *this;
177 }
178
179 /* ----------------------------------------------------------------
180 * Paramater and state extractions and manipulation functions
181 * ---------------------------------------------------------------- */
182
183 void
get(DictionaryDatum & d) const184 nest::aeif_cond_exp::Parameters_::get( DictionaryDatum& d ) const
185 {
186 def< double >( d, names::C_m, C_m );
187 def< double >( d, names::V_th, V_th );
188 def< double >( d, names::t_ref, t_ref_ );
189 def< double >( d, names::g_L, g_L );
190 def< double >( d, names::E_L, E_L );
191 def< double >( d, names::V_reset, V_reset_ );
192 def< double >( d, names::E_ex, E_ex );
193 def< double >( d, names::E_in, E_in );
194 def< double >( d, names::tau_syn_ex, tau_syn_ex );
195 def< double >( d, names::tau_syn_in, tau_syn_in );
196 def< double >( d, names::a, a );
197 def< double >( d, names::b, b );
198 def< double >( d, names::Delta_T, Delta_T );
199 def< double >( d, names::tau_w, tau_w );
200 def< double >( d, names::I_e, I_e );
201 def< double >( d, names::V_peak, V_peak_ );
202 def< double >( d, names::gsl_error_tol, gsl_error_tol );
203 }
204
205 void
set(const DictionaryDatum & d,Node * node)206 nest::aeif_cond_exp::Parameters_::set( const DictionaryDatum& d, Node* node )
207 {
208 updateValueParam< double >( d, names::V_th, V_th, node );
209 updateValueParam< double >( d, names::V_peak, V_peak_, node );
210 updateValueParam< double >( d, names::t_ref, t_ref_, node );
211 updateValueParam< double >( d, names::E_L, E_L, node );
212 updateValueParam< double >( d, names::V_reset, V_reset_, node );
213 updateValueParam< double >( d, names::E_ex, E_ex, node );
214 updateValueParam< double >( d, names::E_in, E_in, node );
215
216 updateValueParam< double >( d, names::C_m, C_m, node );
217 updateValueParam< double >( d, names::g_L, g_L, node );
218
219 updateValueParam< double >( d, names::tau_syn_ex, tau_syn_ex, node );
220 updateValueParam< double >( d, names::tau_syn_in, tau_syn_in, node );
221
222 updateValueParam< double >( d, names::a, a, node );
223 updateValueParam< double >( d, names::b, b, node );
224 updateValueParam< double >( d, names::Delta_T, Delta_T, node );
225 updateValueParam< double >( d, names::tau_w, tau_w, node );
226
227 updateValueParam< double >( d, names::I_e, I_e, node );
228
229 updateValueParam< double >( d, names::gsl_error_tol, gsl_error_tol, node );
230
231 if ( V_peak_ < V_th )
232 {
233 throw BadProperty( "V_peak >= V_th required." );
234 }
235
236 if ( Delta_T < 0. )
237 {
238 throw BadProperty( "Delta_T must be positive." );
239 }
240 else if ( Delta_T > 0. )
241 {
242 // check for possible numerical overflow with the exponential divergence at
243 // spike time, keep a 1e20 margin for the subsequent calculations
244 const double max_exp_arg = std::log( std::numeric_limits< double >::max() / 1e20 );
245 if ( ( V_peak_ - V_th ) / Delta_T >= max_exp_arg )
246 {
247 throw BadProperty(
248 "The current combination of V_peak, V_th and Delta_T"
249 "will lead to numerical overflow at spike time; try"
250 "for instance to increase Delta_T or to reduce V_peak"
251 "to avoid this problem." );
252 }
253 }
254
255 if ( V_reset_ >= V_peak_ )
256 {
257 throw BadProperty( "Ensure that: V_reset < V_peak ." );
258 }
259
260 if ( C_m <= 0 )
261 {
262 throw BadProperty( "Ensure that C_m >0" );
263 }
264
265 if ( t_ref_ < 0 )
266 {
267 throw BadProperty( "Refractory time cannot be negative." );
268 }
269
270 if ( tau_syn_ex <= 0 || tau_syn_in <= 0 || tau_w <= 0 )
271 {
272 throw BadProperty( "All time constants must be strictly positive." );
273 }
274
275 if ( gsl_error_tol <= 0. )
276 {
277 throw BadProperty( "The gsl_error_tol must be strictly positive." );
278 }
279 }
280
281 void
get(DictionaryDatum & d) const282 nest::aeif_cond_exp::State_::get( DictionaryDatum& d ) const
283 {
284 def< double >( d, names::V_m, y_[ V_M ] );
285 def< double >( d, names::g_ex, y_[ G_EXC ] );
286 def< double >( d, names::g_in, y_[ G_INH ] );
287 def< double >( d, names::w, y_[ W ] );
288 }
289
290 void
set(const DictionaryDatum & d,const Parameters_ &,Node * node)291 nest::aeif_cond_exp::State_::set( const DictionaryDatum& d, const Parameters_&, Node* node )
292 {
293 updateValueParam< double >( d, names::V_m, y_[ V_M ], node );
294 updateValueParam< double >( d, names::g_ex, y_[ G_EXC ], node );
295 updateValueParam< double >( d, names::g_in, y_[ G_INH ], node );
296 updateValueParam< double >( d, names::w, y_[ W ], node );
297 if ( y_[ G_EXC ] < 0 || y_[ G_INH ] < 0 )
298 {
299 throw BadProperty( "Conductances must not be negative." );
300 }
301 }
302
Buffers_(aeif_cond_exp & n)303 nest::aeif_cond_exp::Buffers_::Buffers_( aeif_cond_exp& n )
304 : logger_( n )
305 , s_( 0 )
306 , c_( 0 )
307 , e_( 0 )
308 {
309 // Initialization of the remaining members is deferred to
310 // init_buffers_().
311 }
312
Buffers_(const Buffers_ &,aeif_cond_exp & n)313 nest::aeif_cond_exp::Buffers_::Buffers_( const Buffers_&, aeif_cond_exp& n )
314 : logger_( n )
315 , s_( 0 )
316 , c_( 0 )
317 , e_( 0 )
318 {
319 // Initialization of the remaining members is deferred to
320 // init_buffers_().
321 }
322
323 /* ----------------------------------------------------------------
324 * Default and copy constructor for node, and destructor
325 * ---------------------------------------------------------------- */
326
aeif_cond_exp()327 nest::aeif_cond_exp::aeif_cond_exp()
328 : ArchivingNode()
329 , P_()
330 , S_( P_ )
331 , B_( *this )
332 {
333 recordablesMap_.create();
334 }
335
aeif_cond_exp(const aeif_cond_exp & n)336 nest::aeif_cond_exp::aeif_cond_exp( const aeif_cond_exp& n )
337 : ArchivingNode( n )
338 , P_( n.P_ )
339 , S_( n.S_ )
340 , B_( n.B_, *this )
341 {
342 }
343
~aeif_cond_exp()344 nest::aeif_cond_exp::~aeif_cond_exp()
345 {
346 // GSL structs may not have been allocated, so we need to protect destruction
347 if ( B_.s_ )
348 {
349 gsl_odeiv_step_free( B_.s_ );
350 }
351 if ( B_.c_ )
352 {
353 gsl_odeiv_control_free( B_.c_ );
354 }
355 if ( B_.e_ )
356 {
357 gsl_odeiv_evolve_free( B_.e_ );
358 }
359 }
360
361 /* ----------------------------------------------------------------
362 * Node initialization functions
363 * ---------------------------------------------------------------- */
364
365 void
init_buffers_()366 nest::aeif_cond_exp::init_buffers_()
367 {
368 B_.spike_exc_.clear(); // includes resize
369 B_.spike_inh_.clear(); // includes resize
370 B_.currents_.clear(); // includes resize
371 ArchivingNode::clear_history();
372
373 B_.logger_.reset();
374
375 B_.step_ = Time::get_resolution().get_ms();
376
377 // We must integrate this model with high-precision to obtain decent results
378 B_.IntegrationStep_ = std::min( 0.01, B_.step_ );
379
380 if ( B_.s_ == 0 )
381 {
382 B_.s_ = gsl_odeiv_step_alloc( gsl_odeiv_step_rkf45, State_::STATE_VEC_SIZE );
383 }
384 else
385 {
386 gsl_odeiv_step_reset( B_.s_ );
387 }
388
389 if ( B_.c_ == 0 )
390 {
391 B_.c_ = gsl_odeiv_control_yp_new( P_.gsl_error_tol, P_.gsl_error_tol );
392 }
393 else
394 {
395 gsl_odeiv_control_init( B_.c_, P_.gsl_error_tol, P_.gsl_error_tol, 0.0, 1.0 );
396 }
397
398 if ( B_.e_ == 0 )
399 {
400 B_.e_ = gsl_odeiv_evolve_alloc( State_::STATE_VEC_SIZE );
401 }
402 else
403 {
404 gsl_odeiv_evolve_reset( B_.e_ );
405 }
406
407 B_.sys_.jacobian = NULL;
408 B_.sys_.dimension = State_::STATE_VEC_SIZE;
409 B_.sys_.params = reinterpret_cast< void* >( this );
410 B_.sys_.function = aeif_cond_exp_dynamics;
411
412 B_.I_stim_ = 0.0;
413 }
414
415 void
calibrate()416 nest::aeif_cond_exp::calibrate()
417 {
418 // ensures initialization in case mm connected after Simulate
419 B_.logger_.init();
420
421 // set the right threshold and GSL function depending on Delta_T
422 if ( P_.Delta_T > 0. )
423 {
424 V_.V_peak = P_.V_peak_;
425 }
426 else
427 {
428 V_.V_peak = P_.V_th; // same as IAF dynamics for spikes if Delta_T == 0.
429 }
430
431 V_.refractory_counts_ = Time( Time::ms( P_.t_ref_ ) ).get_steps();
432 }
433
434 /* ----------------------------------------------------------------
435 * Update and spike handling functions
436 * ---------------------------------------------------------------- */
437
438 void
update(const Time & origin,const long from,const long to)439 nest::aeif_cond_exp::update( const Time& origin, const long from, const long to )
440 {
441 assert( to >= 0 && ( delay ) from < kernel().connection_manager.get_min_delay() );
442 assert( from < to );
443 assert( State_::V_M == 0 );
444
445 for ( long lag = from; lag < to; ++lag )
446 {
447 double t = 0.0;
448
449 // numerical integration with adaptive step size control:
450 // ------------------------------------------------------
451 // gsl_odeiv_evolve_apply performs only a single numerical
452 // integration step, starting from t and bounded by step;
453 // the while-loop ensures integration over the whole simulation
454 // step (0, step] if more than one integration step is needed due
455 // to a small integration step size;
456 // note that (t+IntegrationStep > step) leads to integration over
457 // (t, step] and afterwards setting t to step, but it does not
458 // enforce setting IntegrationStep to step-t
459 while ( t < B_.step_ )
460 {
461 const int status = gsl_odeiv_evolve_apply( B_.e_,
462 B_.c_,
463 B_.s_,
464 &B_.sys_, // system of ODE
465 &t, // from t
466 B_.step_, // to t <= step
467 &B_.IntegrationStep_, // integration step size
468 S_.y_ ); // neuronal state
469 if ( status != GSL_SUCCESS )
470 {
471 throw GSLSolverFailure( get_name(), status );
472 }
473
474 // check for unreasonable values; we allow V_M to explode
475 if ( S_.y_[ State_::V_M ] < -1e3 || S_.y_[ State_::W ] < -1e6 || S_.y_[ State_::W ] > 1e6 )
476 {
477 throw NumericalInstability( get_name() );
478 }
479
480 // spikes are handled inside the while-loop
481 // due to spike-driven adaptation
482 if ( S_.r_ > 0 )
483 {
484 S_.y_[ State_::V_M ] = P_.V_reset_;
485 }
486 else if ( S_.y_[ State_::V_M ] >= V_.V_peak )
487 {
488 S_.y_[ State_::V_M ] = P_.V_reset_;
489 S_.y_[ State_::W ] += P_.b; // spike-driven adaptation
490
491 /* Initialize refractory step counter.
492 * - We need to add 1 to compensate for count-down immediately after
493 * while loop.
494 * - If neuron has no refractory time, set to 0 to avoid refractory
495 * artifact inside while loop.
496 */
497 S_.r_ = V_.refractory_counts_ > 0 ? V_.refractory_counts_ + 1 : 0;
498
499 set_spiketime( Time::step( origin.get_steps() + lag + 1 ) );
500 SpikeEvent se;
501 kernel().event_delivery_manager.send( *this, se, lag );
502 }
503 }
504
505 // decrement refractory count
506 if ( S_.r_ > 0 )
507 {
508 --S_.r_;
509 }
510
511 // apply spikes
512 S_.y_[ State_::G_EXC ] += B_.spike_exc_.get_value( lag );
513 S_.y_[ State_::G_INH ] += B_.spike_inh_.get_value( lag );
514
515 // set new input current
516 B_.I_stim_ = B_.currents_.get_value( lag );
517
518 // log state data
519 B_.logger_.record_data( origin.get_steps() + lag );
520 }
521 }
522
523 void
handle(SpikeEvent & e)524 nest::aeif_cond_exp::handle( SpikeEvent& e )
525 {
526 assert( e.get_delay_steps() > 0 );
527
528 if ( e.get_weight() > 0.0 )
529 {
530 B_.spike_exc_.add_value( e.get_rel_delivery_steps( kernel().simulation_manager.get_slice_origin() ),
531 e.get_weight() * e.get_multiplicity() );
532 }
533 else
534 {
535 B_.spike_inh_.add_value( e.get_rel_delivery_steps( kernel().simulation_manager.get_slice_origin() ),
536 -e.get_weight() * e.get_multiplicity() );
537 }
538 }
539
540 void
handle(CurrentEvent & e)541 nest::aeif_cond_exp::handle( CurrentEvent& e )
542 {
543 assert( e.get_delay_steps() > 0 );
544
545 const double c = e.get_current();
546 const double w = e.get_weight();
547
548 B_.currents_.add_value( e.get_rel_delivery_steps( kernel().simulation_manager.get_slice_origin() ), w * c );
549 }
550
551 void
handle(DataLoggingRequest & e)552 nest::aeif_cond_exp::handle( DataLoggingRequest& e )
553 {
554 B_.logger_.handle( e );
555 }
556
557 #endif // HAVE_GSL
558