1 /*
2 * iaf_psc_alpha.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 IAF_PSC_ALPHA_H
24 #define IAF_PSC_ALPHA_H
25
26 // Includes from nestkernel:
27 #include "archiving_node.h"
28 #include "connection.h"
29 #include "event.h"
30 #include "nest_types.h"
31 #include "recordables_map.h"
32 #include "ring_buffer.h"
33 #include "universal_data_logger.h"
34
35 namespace nest
36 {
37
38 /* BeginUserDocs: neuron, integrate-and-fire, current-based
39
40 Short description
41 +++++++++++++++++
42
43 Leaky integrate-and-fire neuron model
44
45 Description
46 +++++++++++
47
48 iaf_psc_alpha is an implementation of a leaky integrate-and-fire model
49 with alpha-function shaped synaptic currents. Thus, synaptic currents
50 and the resulting postsynaptic potentials have a finite rise time.
51
52 The threshold crossing is followed by an absolute refractory period
53 during which the membrane potential is clamped to the resting potential.
54
55 The linear subthreshold dynamics is integrated by the Exact
56 Integration scheme [1]_. The neuron dynamics is solved on the time
57 grid given by the computation step size. Incoming as well as emitted
58 spikes are forced to that grid.
59
60 An additional state variable and the corresponding differential
61 equation represents a piecewise constant external current.
62
63 The general framework for the consistent formulation of systems with
64 neuron like dynamics interacting by point events is described in
65 [1]_. A flow chart can be found in [2]_.
66
67 Critical tests for the formulation of the neuron model are the
68 comparisons of simulation results for different computation step
69 sizes. sli/testsuite/nest contains a number of such tests.
70
71 The iaf_psc_alpha is the standard model used to check the consistency
72 of the nest simulation kernel because it is at the same time complex
73 enough to exhibit non-trivial dynamics and simple enough compute
74 relevant measures analytically.
75
76 .. note::
77 The present implementation uses individual variables for the
78 components of the state vector and the non-zero matrix elements of
79 the propagator. Because the propagator is a lower triangular matrix,
80 no full matrix multiplication needs to be carried out and the
81 computation can be done "in place", i.e. no temporary state vector
82 object is required.
83
84 The template support of recent C++ compilers enables a more succinct
85 formulation without loss of runtime performance already at minimal
86 optimization levels. A future version of iaf_psc_alpha will probably
87 address the problem of efficient usage of appropriate vector and
88 matrix objects.
89
90 .. note::
91
92 If `tau_m` is very close to `tau_syn_ex` or `tau_syn_in`, the model
93 will numerically behave as if `tau_m` is equal to `tau_syn_ex` or
94 `tau_syn_in`, respectively, to avoid numerical instabilities.
95
96 For implementation details see the
97 `IAF_neurons_singularity <../model_details/IAF_neurons_singularity.ipynb>`_ notebook.
98
99 Parameters
100 ++++++++++
101
102 The following parameters can be set in the status dictionary.
103
104 =========== ====== ==========================================================
105 V_m mV Membrane potential
106 E_L mV Resting membrane potential
107 C_m pF Capacity of the membrane
108 tau_m ms Membrane time constant
109 t_ref ms Duration of refractory period
110 V_th mV Spike threshold
111 V_reset mV Reset potential of the membrane
112 tau_syn_ex ms Rise time of the excitatory synaptic alpha function
113 tau_syn_in ms Rise time of the inhibitory synaptic alpha function
114 I_e pA Constant input current
115 V_min mV Absolute lower value for the membrane potenial
116 =========== ====== ==========================================================
117
118 References
119 ++++++++++
120
121 .. [1] Rotter S, Diesmann M (1999). Exact simulation of
122 time-invariant linear systems with applications to neuronal
123 modeling. Biologial Cybernetics 81:381-402.
124 DOI: https://doi.org/10.1007/s004220050570
125 .. [2] Diesmann M, Gewaltig M-O, Rotter S, & Aertsen A (2001). State
126 space analysis of synchronous spiking in cortical neural
127 networks. Neurocomputing 38-40:565-571.
128 DOI: https://doi.org/10.1016/S0925-2312(01)00409-X
129 .. [3] Morrison A, Straube S, Plesser H E, Diesmann M (2006). Exact
130 subthreshold integration with continuous spike times in discrete time
131 neural network simulations. Neural Computation, in press
132 DOI: https://doi.org/10.1162/neco.2007.19.1.47
133
134 Sends
135 +++++
136
137 SpikeEvent
138
139 Receives
140 ++++++++
141
142 SpikeEvent, CurrentEvent, DataLoggingRequest
143
144 See also
145 ++++++++
146
147 iaf_psc_delta, iaf_psc_exp, iaf_cond_exp
148
149 EndUserDocs */
150
151 class iaf_psc_alpha : public ArchivingNode
152 {
153
154 public:
155 iaf_psc_alpha();
156 iaf_psc_alpha( const iaf_psc_alpha& );
157
158 /**
159 * Import sets of overloaded virtual functions.
160 * @see Technical Issues / Virtual Functions: Overriding, Overloading, and
161 * Hiding
162 */
163 using Node::handle;
164 using Node::handles_test_event;
165
166 port send_test_event( Node&, rport, synindex, bool );
167
168 void handle( SpikeEvent& );
169 void handle( CurrentEvent& );
170 void handle( DataLoggingRequest& );
171
172 port handles_test_event( SpikeEvent&, rport );
173 port handles_test_event( CurrentEvent&, rport );
174 port handles_test_event( DataLoggingRequest&, rport );
175
176 void get_status( DictionaryDatum& ) const;
177 void set_status( const DictionaryDatum& );
178
179 private:
180 void init_buffers_();
181 void calibrate();
182
183 void update( Time const&, const long, const long );
184
185 // The next two classes need to be friends to access the State_ class/member
186 friend class RecordablesMap< iaf_psc_alpha >;
187 friend class UniversalDataLogger< iaf_psc_alpha >;
188
189 // ----------------------------------------------------------------
190
191 struct Parameters_
192 {
193 /** Membrane time constant in ms. */
194 double Tau_;
195
196 /** Membrane capacitance in pF. */
197 double C_;
198
199 /** Refractory period in ms. */
200 double TauR_;
201
202 /** Resting potential in mV. */
203 double E_L_;
204
205 /** External current in pA */
206 double I_e_;
207
208 /** Reset value of the membrane potential */
209 double V_reset_;
210
211 /** Threshold, RELATIVE TO RESTING POTENTIAL(!).
212 I.e. the real threshold is (E_L_+Theta_). */
213 double Theta_;
214
215 /** Lower bound, RELATIVE TO RESTING POTENTIAL(!).
216 I.e. the real lower bound is (LowerBound_+E_L_). */
217 double LowerBound_;
218
219 /** Time constant of excitatory synaptic current in ms. */
220 double tau_ex_;
221
222 /** Time constant of inhibitory synaptic current in ms. */
223 double tau_in_;
224
225 Parameters_(); //!< Sets default parameter values
226
227 void get( DictionaryDatum& ) const; //!< Store current values in dictionary
228
229 /** Set values from dictionary.
230 * @returns Change in reversal potential E_L, to be passed to State_::set()
231 */
232 double set( const DictionaryDatum&, Node* node );
233 };
234
235 // ----------------------------------------------------------------
236
237 struct State_
238 {
239 double y0_; //!< Constant current
240 double dI_ex_;
241 double I_ex_;
242 double dI_in_;
243 double I_in_;
244 //! This is the membrane potential RELATIVE TO RESTING POTENTIAL.
245 double y3_;
246
247 int r_; //!< Number of refractory steps remaining
248
249 State_(); //!< Default initialization
250
251 void get( DictionaryDatum&, const Parameters_& ) const;
252
253 /** Set values from dictionary.
254 * @param dictionary to take data from
255 * @param current parameters
256 * @param Change in reversal potential E_L specified by this dict
257 */
258 void set( const DictionaryDatum&, const Parameters_&, double, Node* node );
259 };
260
261 // ----------------------------------------------------------------
262
263 struct Buffers_
264 {
265
266 Buffers_( iaf_psc_alpha& );
267 Buffers_( const Buffers_&, iaf_psc_alpha& );
268
269 //! Indices for access to different channels of input_buffer_
270 enum
271 {
272 SYN_IN = 0,
273 SYN_EX,
274 I0,
275 NUM_INPUT_CHANNELS
276 };
277
278 /** buffers and sums up incoming spikes/currents */
279 MultiChannelInputBuffer< NUM_INPUT_CHANNELS > input_buffer_;
280
281 //! Logger for all analog data
282 UniversalDataLogger< iaf_psc_alpha > logger_;
283 };
284
285 // ----------------------------------------------------------------
286
287 struct Variables_
288 {
289
290 /** Amplitude of the synaptic current.
291 This value is chosen such that a postsynaptic potential with
292 weight one has an amplitude of 1 mV.
293 */
294 double EPSCInitialValue_;
295 double IPSCInitialValue_;
296 int RefractoryCounts_;
297
298 double P11_ex_;
299 double P21_ex_;
300 double P22_ex_;
301 double P31_ex_;
302 double P32_ex_;
303 double P11_in_;
304 double P21_in_;
305 double P22_in_;
306 double P31_in_;
307 double P32_in_;
308 double P30_;
309 double P33_;
310 double expm1_tau_m_;
311
312 double weighted_spikes_ex_;
313 double weighted_spikes_in_;
314 };
315
316 // Access functions for UniversalDataLogger -------------------------------
317
318 //! Read out the real membrane potential
319 inline double
get_V_m_()320 get_V_m_() const
321 {
322 return S_.y3_ + P_.E_L_;
323 }
324
325 inline double
get_I_syn_ex_()326 get_I_syn_ex_() const
327 {
328 return S_.I_ex_;
329 }
330
331 inline double
get_I_syn_in_()332 get_I_syn_in_() const
333 {
334 return S_.I_in_;
335 }
336
337 // Data members -----------------------------------------------------------
338
339 /**
340 * @defgroup iaf_psc_alpha_data
341 * Instances of private data structures for the different types
342 * of data pertaining to the model.
343 * @note The order of definitions is important for speed.
344 * @{
345 */
346 Parameters_ P_;
347 State_ S_;
348 Variables_ V_;
349 Buffers_ B_;
350 /** @} */
351
352 //! Mapping of recordables names to access functions
353 static RecordablesMap< iaf_psc_alpha > recordablesMap_;
354 };
355
356 inline port
send_test_event(Node & target,rport receptor_type,synindex,bool)357 nest::iaf_psc_alpha::send_test_event( Node& target, rport receptor_type, synindex, bool )
358 {
359 SpikeEvent e;
360 e.set_sender( *this );
361 return target.handles_test_event( e, receptor_type );
362 }
363
364 inline port
handles_test_event(SpikeEvent &,rport receptor_type)365 iaf_psc_alpha::handles_test_event( SpikeEvent&, rport receptor_type )
366 {
367 if ( receptor_type != 0 )
368 {
369 throw UnknownReceptorType( receptor_type, get_name() );
370 }
371 return 0;
372 }
373
374 inline port
handles_test_event(CurrentEvent &,rport receptor_type)375 iaf_psc_alpha::handles_test_event( CurrentEvent&, 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(DataLoggingRequest & dlr,rport receptor_type)385 iaf_psc_alpha::handles_test_event( DataLoggingRequest& dlr, rport receptor_type )
386 {
387 if ( receptor_type != 0 )
388 {
389 throw UnknownReceptorType( receptor_type, get_name() );
390 }
391 return B_.logger_.connect_logging_device( dlr, recordablesMap_ );
392 }
393
394 inline void
get_status(DictionaryDatum & d)395 iaf_psc_alpha::get_status( DictionaryDatum& d ) const
396 {
397 P_.get( d );
398 S_.get( d, P_ );
399 ArchivingNode::get_status( d );
400
401 ( *d )[ names::recordables ] = recordablesMap_.get_list();
402 }
403
404 inline void
set_status(const DictionaryDatum & d)405 iaf_psc_alpha::set_status( const DictionaryDatum& d )
406 {
407 Parameters_ ptmp = P_; // temporary copy in case of errors
408 const double delta_EL = ptmp.set( d, this ); // throws if BadProperty
409 State_ stmp = S_; // temporary copy in case of errors
410 stmp.set( d, ptmp, delta_EL, this ); // throws if BadProperty
411
412 // We now know that (ptmp, stmp) are consistent. We do not
413 // write them back to (P_, S_) before we are also sure that
414 // the properties to be set in the parent class are internally
415 // consistent.
416 ArchivingNode::set_status( d );
417
418 // if we get here, temporaries contain consistent set of properties
419 P_ = ptmp;
420 S_ = stmp;
421 }
422
423 } // namespace
424
425 #endif /* #ifndef IAF_PSC_ALPHA_H */
426