1 /*
2 * aeif_psc_exp.h
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 #ifndef AEIF_PSC_EXP_H
24 #define AEIF_PSC_EXP_H
25
26 // Generated includes:
27 #include "config.h"
28
29 #ifdef HAVE_GSL
30
31 // External includes:
32 #include <gsl/gsl_errno.h>
33 #include <gsl/gsl_matrix.h>
34 #include <gsl/gsl_odeiv.h>
35
36 // Includes from nestkernel:
37 #include "archiving_node.h"
38 #include "connection.h"
39 #include "event.h"
40 #include "nest_types.h"
41 #include "recordables_map.h"
42 #include "ring_buffer.h"
43 #include "universal_data_logger.h"
44
45
46 namespace nest
47 {
48 /**
49 * Function computing right-hand side of ODE for GSL solver.
50 * @note Must be declared here so we can befriend it in class.
51 * @note Must have C-linkage for passing to GSL. Internally, it is
52 * a first-class C++ function, but cannot be a member function
53 * because of the C-linkage.
54 * @note No point in declaring it inline, since it is called
55 * through a function pointer.
56 * @param void* Pointer to model neuron instance.
57 */
58 extern "C" int aeif_psc_exp_dynamics( double, const double*, double*, void* );
59
60 /* BeginUserDocs: neuron, integrate-and-fire, adaptive threshold, current-based
61
62 Short description
63 +++++++++++++++++
64
65 Current-based exponential integrate-and-fire neuron model
66
67 Description
68 +++++++++++
69
70 aeif_psc_exp is the adaptive exponential integrate and fire neuron
71 according to Brette and Gerstner (2005), with postsynaptic currents
72 in the form of truncated exponentials.
73
74 This implementation uses the embedded 4th order Runge-Kutta-Fehlberg
75 solver with adaptive stepsize to integrate the differential equation.
76
77 The membrane potential is given by the following differential equation:
78
79 .. math::
80
81 C dV/dt= -g_L(V-E_L)+g_L*\Delta_T*\exp((V-V_T)/\Delta_T)-g_e(t)(V-E_e) \\
82 -g_i(t)(V-E_i)-w +I_e
83
84 and
85
86 .. math::
87
88 \tau_w * dw/dt= a(V-E_L) -W
89
90
91 Note that the spike detection threshold V_peak is automatically set to
92 :math:`V_th+10` mV to avoid numerical instabilites that may result from
93 setting V_peak too high.
94
95 For implementation details see the
96 `aeif_models_implementation <../model_details/aeif_models_implementation.ipynb>`_ notebook.
97
98 Parameters
99 ++++++++++
100
101 The following parameters can be set in the status dictionary.
102
103 ======== ======= =======================================
104 **Dynamic state variables:**
105 --------------------------------------------------------
106 V_m mV Membrane potential
107 I_ex pA Excitatory synaptic current
108 I_in pA Inhibitory synaptic current
109 w pA Spike-adaptation current
110 ======== ======= =======================================
111
112 ======== ======= =======================================
113 **Membrane Parameters**
114 --------------------------------------------------------
115 C_m pF Capacity of the membrane
116 t_ref ms Duration of refractory period
117 V_reset mV Reset value for V_m after a spike
118 E_L mV Leak reversal potential
119 g_L nS Leak conductance
120 I_e pA Constant external input current
121 ======== ======= =======================================
122
123 ======== ======= ==================================
124 **Spike adaptation parameters**
125 ---------------------------------------------------
126 a ns Subthreshold adaptation
127 b pA Spike-triggered adaptation
128 Delta_T mV Slope factor
129 tau_w ms Adaptation time constant
130 V_t mV Spike initiation threshold
131 V_peak mV Spike detection threshold
132 ======== ======= ==================================
133
134 =========== ======= ===========================================================
135 **Synaptic parameters**
136 -------------------------------------------------------------------------------
137 tau_syn_ex ms Exponential decay time constant of excitatory synaptic
138 conductance kernel
139 tau_syn_in ms Exponential decay time constant of inhibitory synaptic
140 conductance kernel
141 =========== ======= ===========================================================
142
143 ============= ======= =========================================================
144 **Integration parameters**
145 -------------------------------------------------------------------------------
146 gsl_error_tol real This parameter controls the admissible error of the
147 GSL integrator. Reduce it if NEST complains about
148 numerical instabilities
149 ============= ======= =========================================================
150
151 Sends
152 +++++
153
154 SpikeEvent
155
156 Receives
157 ++++++++
158
159 SpikeEvent, CurrentEvent, DataLoggingRequest
160
161 References
162 ++++++++++
163
164 .. [1] Brette R and Gerstner W (2005). Adaptive Exponential
165 Integrate-and-Fire Model as an Effective Description of Neuronal
166 Activity. J Neurophysiol 94:3637-3642.
167 DOI: https://doi.org/10.1152/jn.00686.2005
168
169 See also
170 ++++++++
171
172 iaf_psc_exp, aeif_cond_exp
173
174 EndUserDocs */
175
176 class aeif_psc_exp : public ArchivingNode
177 {
178
179 public:
180 aeif_psc_exp();
181 aeif_psc_exp( const aeif_psc_exp& );
182 ~aeif_psc_exp();
183
184 /**
185 * Import sets of overloaded virtual functions.
186 * @see Technical Issues / Virtual Functions: Overriding, Overloading, and
187 * Hiding
188 */
189 using Node::handle;
190 using Node::handles_test_event;
191
192 port send_test_event( Node&, rport, synindex, bool );
193
194 void handle( SpikeEvent& );
195 void handle( CurrentEvent& );
196 void handle( DataLoggingRequest& );
197
198 port handles_test_event( SpikeEvent&, rport );
199 port handles_test_event( CurrentEvent&, rport );
200 port handles_test_event( DataLoggingRequest&, rport );
201
202 void get_status( DictionaryDatum& ) const;
203 void set_status( const DictionaryDatum& );
204
205 private:
206 void init_buffers_();
207 void calibrate();
208 void update( const Time&, const long, const long );
209
210 // END Boilerplate function declarations ----------------------------
211
212 // Friends --------------------------------------------------------
213
214 // make dynamics function quasi-member
215 friend int aeif_psc_exp_dynamics( double, const double*, double*, void* );
216
217 // The next two classes need to be friends to access the State_ class/member
218 friend class RecordablesMap< aeif_psc_exp >;
219 friend class UniversalDataLogger< aeif_psc_exp >;
220
221 private:
222 // ----------------------------------------------------------------
223
224 //! Independent parameters
225 struct Parameters_
226 {
227 double V_peak_; //!< Spike detection threshold in mV
228 double V_reset_; //!< Reset Potential in mV
229 double t_ref_; //!< Refractory period in ms
230
231 double g_L; //!< Leak Conductance in nS
232 double C_m; //!< Membrane Capacitance in pF
233 double E_L; //!< Leak reversal Potential (aka resting potential) in mV
234 double Delta_T; //!< Slope factor in ms
235 double tau_w; //!< Adaptation time-constant in ms
236 double a; //!< Subthreshold adaptation in nS
237 double b; //!< Spike-triggered adaptation in pA
238 double V_th; //!< Spike threshold in mV
239 double tau_syn_ex; //!< Excitatory synaptic kernel decay time in ms
240 double tau_syn_in; //!< Inhibitory synaptic kernel decay time in ms
241 double I_e; //!< Intrinsic current in pA
242
243 double gsl_error_tol; //!< Error bound for GSL integrator
244
245 Parameters_(); //!< Sets default parameter values
246
247 void get( DictionaryDatum& ) const; //!< Store current values in dictionary
248 void set( const DictionaryDatum&, Node* node ); //!< Set values from dicitonary
249 };
250
251 public:
252 // ----------------------------------------------------------------
253
254 /**
255 * State variables of the model.
256 * @note Copy constructor required because of C-style array.
257 */
258 struct State_
259 {
260 /**
261 * Enumeration identifying elements in state array State_::y_.
262 * The state vector must be passed to GSL as a C array. This enum
263 * identifies the elements of the vector. It must be public to be
264 * accessible from the iteration function.
265 */
266 enum StateVecElems
267 {
268 V_M = 0,
269 I_EXC, // 1
270 I_INH, // 2
271 W, // 3
272 STATE_VEC_SIZE
273 };
274
275 //! neuron state, must be C-array for GSL solver
276 double y_[ STATE_VEC_SIZE ];
277 unsigned int r_; //!< number of refractory steps remaining
278
279 State_( const Parameters_& ); //!< Default initialization
280 State_( const State_& );
281
282 State_& operator=( const State_& );
283
284 void get( DictionaryDatum& ) const;
285 void set( const DictionaryDatum&, const Parameters_&, Node* );
286 };
287
288 // ----------------------------------------------------------------
289
290 /**
291 * Buffers of the model.
292 */
293 struct Buffers_
294 {
295 Buffers_( aeif_psc_exp& ); //!<Sets buffer pointers to 0
296 Buffers_( const Buffers_&, aeif_psc_exp& ); //!<Sets buffer pointers to 0
297
298 //! Logger for all analog data
299 UniversalDataLogger< aeif_psc_exp > logger_;
300
301 /** buffers and sums up incoming spikes/currents */
302 RingBuffer spike_exc_;
303 RingBuffer spike_inh_;
304 RingBuffer currents_;
305
306 /** GSL ODE stuff */
307 gsl_odeiv_step* s_; //!< stepping function
308 gsl_odeiv_control* c_; //!< adaptive stepsize control function
309 gsl_odeiv_evolve* e_; //!< evolution function
310 gsl_odeiv_system sys_; //!< struct describing the GSL system
311
312 // Since IntergrationStep_ is initialized with step_, and the resolution
313 // cannot change after nodes have been created, it is safe to place both
314 // here.
315 double step_; //!< step size in ms
316 double IntegrationStep_; //!< current integration time step, updated by GSL
317
318 /**
319 * Input current injected by CurrentEvent.
320 * This variable is used to transport the current applied into the
321 * _dynamics function computing the derivative of the state vector.
322 * It must be a part of Buffers_, since it is initialized once before
323 * the first simulation, but not modified before later Simulate calls.
324 */
325 double I_stim_;
326 };
327
328 // ----------------------------------------------------------------
329
330 /**
331 * Internal variables of the model.
332 */
333 struct Variables_
334 {
335 /**
336 * Threshold detection for spike events: P.V_peak if Delta_T > 0.,
337 * P.V_th if Delta_T == 0.
338 */
339 double V_peak;
340
341 unsigned int refractory_counts_;
342 };
343
344 // Access functions for UniversalDataLogger -------------------------------
345
346 //! Read out state vector elements, used by UniversalDataLogger
347 template < State_::StateVecElems elem >
348 double
get_y_elem_()349 get_y_elem_() const
350 {
351 return S_.y_[ elem ];
352 }
353
354 // ----------------------------------------------------------------
355
356 Parameters_ P_;
357 State_ S_;
358 Variables_ V_;
359 Buffers_ B_;
360
361 //! Mapping of recordables names to access functions
362 static RecordablesMap< aeif_psc_exp > recordablesMap_;
363 };
364
365 inline port
send_test_event(Node & target,rport receptor_type,synindex,bool)366 aeif_psc_exp::send_test_event( Node& target, rport receptor_type, synindex, bool )
367 {
368 SpikeEvent e;
369 e.set_sender( *this );
370
371 return target.handles_test_event( e, receptor_type );
372 }
373
374 inline port
handles_test_event(SpikeEvent &,rport receptor_type)375 aeif_psc_exp::handles_test_event( SpikeEvent&, rport receptor_type )
376 {
377 if ( receptor_type != 0 )
378 {
379 throw UnknownReceptorType( receptor_type, get_name() );
380 }
381 return 0;
382 }
383
384 inline port
handles_test_event(CurrentEvent &,rport receptor_type)385 aeif_psc_exp::handles_test_event( CurrentEvent&, rport receptor_type )
386 {
387 if ( receptor_type != 0 )
388 {
389 throw UnknownReceptorType( receptor_type, get_name() );
390 }
391 return 0;
392 }
393
394 inline port
handles_test_event(DataLoggingRequest & dlr,rport receptor_type)395 aeif_psc_exp::handles_test_event( DataLoggingRequest& dlr, rport receptor_type )
396 {
397 if ( receptor_type != 0 )
398 {
399 throw UnknownReceptorType( receptor_type, get_name() );
400 }
401 return B_.logger_.connect_logging_device( dlr, recordablesMap_ );
402 }
403
404 inline void
get_status(DictionaryDatum & d)405 aeif_psc_exp::get_status( DictionaryDatum& d ) const
406 {
407 P_.get( d );
408 S_.get( d );
409 ArchivingNode::get_status( d );
410
411 ( *d )[ names::recordables ] = recordablesMap_.get_list();
412 }
413
414 inline void
set_status(const DictionaryDatum & d)415 aeif_psc_exp::set_status( const DictionaryDatum& d )
416 {
417 Parameters_ ptmp = P_; // temporary copy in case of errors
418 ptmp.set( d, this ); // throws if BadProperty
419 State_ stmp = S_; // temporary copy in case of errors
420 stmp.set( d, ptmp, this ); // throws if BadProperty
421
422 // We now know that (ptmp, stmp) are consistent. We do not
423 // write them back to (P_, S_) before we are also sure that
424 // the properties to be set in the parent class are internally
425 // consistent.
426 ArchivingNode::set_status( d );
427
428 // if we get here, temporaries contain consistent set of properties
429 P_ = ptmp;
430 S_ = stmp;
431 }
432
433 } // namespace
434
435 #endif // HAVE_GSL
436 #endif // AEIF_PSC_EXP_H
437