1 /*
2 * ht_neuron.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 HT_NEURON_H
24 #define HT_NEURON_H
25
26 // Generated includes:
27 #include "config.h"
28
29 #ifdef HAVE_GSL
30
31 // C++ includes:
32 #include <string>
33 #include <vector>
34
35 // C includes:
36 #include <gsl/gsl_errno.h>
37 #include <gsl/gsl_matrix.h>
38 #include <gsl/gsl_odeiv.h>
39
40 // Includes from nestkernel:
41 #include "archiving_node.h"
42 #include "connection.h"
43 #include "recordables_map.h"
44 #include "ring_buffer.h"
45 #include "universal_data_logger.h"
46
47 // Includes from
48
49 // Includes from sli:
50 #include "stringdatum.h"
51
52 namespace nest
53 {
54 /**
55 * Function computing right-hand side of ODE for GSL solver.
56 * @note Must be declared here so we can befriend it in class.
57 * @note Must have C-linkage for passing to GSL. Internally, it is
58 * a first-class C++ function, but cannot be a member function
59 * because of the C-linkage.
60 * @note No point in declaring it inline, since it is called
61 * through a function pointer.
62 * @param void* Pointer to model neuron instance.
63 */
64 extern "C" int ht_neuron_dynamics( double, const double*, double*, void* );
65
66 /* BeginUserDocs: neuron, Hill-Tononi plasticity
67
68 Short description
69 +++++++++++++++++
70
71 Neuron model after Hill & Tononi (2005)
72
73 Description
74 +++++++++++
75
76 This model neuron implements a slightly modified version of the
77 neuron model described in [1]_. The most important properties are:
78
79 - Integrate-and-fire with threshold adaptive threshold.
80 - Repolarizing potassium current instead of hard reset.
81 - AMPA, NMDA, GABA_A, and GABA_B conductance-based synapses with
82 beta-function (difference of exponentials) time course.
83 - Voltage-dependent NMDA with instantaneous or two-stage unblocking [1]_, [2]_.
84 - Intrinsic currents I_h, I_T, I_Na(p), and I_KNa.
85 - Synaptic "minis" are not implemented.
86
87 For implementation details see:
88
89 - `HillTononi_model <../model_details/HillTononiModels.ipynb>`_
90
91 For examples, see:
92
93 - :doc:`../auto_examples/intrinsic_currents_spiking`
94 - :doc:`../auto_examples/intrinsic_currents_subthreshold`
95
96 For an example network model using ``ht_neuron`` (based on [1]_), see:
97
98 - `Multiarea Hill-Tononi thalamocortical network model
99 <https://github.com/ricardomurphy/Multiarea-Hill-Tononi-thalamocortical-network-model>`_
100
101 Parameters
102 ++++++++++
103
104 =============== ======= =========================================================
105 V_m mV Membrane potential
106 tau_m ms Membrane time constant applying to all currents except
107 repolarizing K-current (see [1]_, p 1677)
108 t_ref ms Refractory time and duration of post-spike repolarizing
109 potassium current (t_spike in [1]_)
110 tau_spike ms Membrane time constant for post-spike repolarizing
111 potassium current
112 voltage_clamp boolean If true, clamp voltage to value at beginning of
113 simulation
114 (default: false, mainly for testing)
115 theta mV Threshold
116 theta_eq mV Equilibrium value
117 tau_theta ms Time constant
118 g_KL nS Conductance for potassium leak current
119 E_K mV Reversal potential for potassium leak currents
120 g_NaL nS Conductance for sodium leak currents
121 E_Na mV Reversal potential for Na leak currents
122 tau_D_KNa ms Relaxation time constant for I_KNa
123 receptor_types Dictionary mapping synapse names to ports on neuron model
124 recordables List of recordable quantities
125 =============== ======= =========================================================
126
127 +------------------------------------------------------------+
128 |{E_rev,g_peak,tau_rise,tau_decay}_{AMPA,NMDA,GABA_A,GABA_B} |
129 +------------------------------------------------------------+
130 | Reversal potentials, peak conductances and time constants |
131 | for synapses (tau_rise/tau_decay correspond to tau_1/tau_2 |
132 | in the paper) |
133 +------------------------------------------------------------+
134
135 +------------------------+------------------------------------------------+
136 |V_act_NMDA, S_act_NMDA, | Parameters for voltage dependence of NMDA- |
137 |tau_Mg_{fast, slow}_NMDA| conductance, see above |
138 +------------------------+------------------------------------------------+
139
140 ============================ =================================================
141 instant_unblock_NMDA Instantaneous NMDA unblocking (default: false)
142 {E_rev,g_peak}_{h,T,NaP,KNa} Reversal potential and peak conductance for
143 intrinsic currents
144 {N}_{T,NaP} Exponent activation term m_inf, corresponding to
145 number of activation particles
146 equilibrate If given and true, time-dependent activation
147 and inactivation state variables (h, m) of
148 intrinsic currents and NMDA channels are set
149 to their equilibrium values during this
150 SetStatus call; otherwise they retain their
151 present values.
152 ============================ =================================================
153
154 .. note::
155 Conductances are unitless in this model and currents are in mV.
156
157 Sends
158 +++++
159
160 SpikeEvent
161
162 Receives
163 ++++++++
164
165 SpikeEvent, CurrentEvent, DataLoggingRequest
166
167 References
168 ++++++++++
169
170 .. [1] Hill S, Tononi G (2005). Modeling sleep and wakefulness in the
171 thalamocortical system. Journal of Neurophysiology. 93:1671-1698.
172 DOI: https://doi.org/10.1152/jn.00915.2004
173 .. [2] Vargas-Caballero M, Robinson HPC (2003). A slow fraction of Mg2+
174 unblock of NMDA receptors limits their contribution to spike generation
175 in cortical pyramidal neurons. Journal of Neurophysiology 89:2778-2783.
176 DOI: https://doi.org/10.1152/jn.01038.2002
177
178 See also
179 ++++++++
180
181 ht_synapse
182
183 EndUserDocs */
184
185 class ht_neuron : public ArchivingNode
186 {
187 public:
188 ht_neuron();
189 ht_neuron( const ht_neuron& );
190 ~ht_neuron();
191
192 /**
193 * Import sets of overloaded virtual functions.
194 * @see Technical Issues / Virtual Functions: Overriding, Overloading, and
195 * Hiding
196 */
197 using Node::handle;
198 using Node::handles_test_event;
199
200 port send_test_event( Node&, rport, synindex, bool );
201
202 void handle( SpikeEvent& e );
203 void handle( CurrentEvent& e );
204 void handle( DataLoggingRequest& );
205
206 port handles_test_event( SpikeEvent&, rport );
207 port handles_test_event( CurrentEvent&, rport );
208 port handles_test_event( DataLoggingRequest&, rport );
209
210 void get_status( DictionaryDatum& ) const;
211 void set_status( const DictionaryDatum& );
212
213 private:
214 /**
215 * Synapse types to connect to
216 * @note Excluded upper and lower bounds are defined as INF_, SUP_.
217 * Excluding port 0 avoids accidental connections.
218 */
219 enum SynapseTypes
220 {
221 INF_SPIKE_RECEPTOR = 0,
222 AMPA,
223 NMDA,
224 GABA_A,
225 GABA_B,
226 SUP_SPIKE_RECEPTOR
227 };
228
229 void init_buffers_();
230 void calibrate();
231
232 void update( Time const&, const long, const long );
233
234 double get_synapse_constant( double, double, double );
235
236 // END Boilerplate function declarations ----------------------------
237
238 // Friends --------------------------------------------------------
239
240 // make dynamics function quasi-member
241 friend int ht_neuron_dynamics( double, const double*, double*, void* );
242
243 // ----------------------------------------------------------------
244
245 /**
246 * Independent parameters of the model.
247 */
248 struct Parameters_
249 {
250 Parameters_();
251
252 void get( DictionaryDatum& ) const; //!< Store current values in dictionary
253 void set( const DictionaryDatum&, Node* node ); //!< Set values from dicitonary
254
255 // Note: Conductances are unitless
256 // Leaks
257 double E_Na; // mV
258 double E_K; // mV
259 double g_NaL;
260 double g_KL;
261 double tau_m; // ms
262
263 // Dynamic threshold
264 double theta_eq; // mV
265 double tau_theta; // ms
266
267 // Post-spike potassium current
268 double tau_spike; // ms, membrane time constant for this current
269 double t_ref; // ms, refractory time
270
271 // Parameters for synapse of type AMPA, GABA_A, GABA_B and NMDA
272 double g_peak_AMPA;
273 double tau_rise_AMPA; // ms
274 double tau_decay_AMPA; // ms
275 double E_rev_AMPA; // mV
276
277 double g_peak_NMDA;
278 double tau_rise_NMDA; // ms
279 double tau_decay_NMDA; // ms
280 double E_rev_NMDA; // mV
281 double V_act_NMDA; // mV, inactive for V << Vact, inflection of sigmoid
282 double S_act_NMDA; // mV, scale of inactivation
283 double tau_Mg_slow_NMDA; // ms
284 double tau_Mg_fast_NMDA; // ms
285 bool instant_unblock_NMDA;
286
287 double g_peak_GABA_A;
288 double tau_rise_GABA_A; // ms
289 double tau_decay_GABA_A; // ms
290 double E_rev_GABA_A; // mV
291
292 double g_peak_GABA_B;
293 double tau_rise_GABA_B; // ms
294 double tau_decay_GABA_B; // ms
295 double E_rev_GABA_B; // mV
296
297 // parameters for intrinsic currents
298 double g_peak_NaP;
299 double E_rev_NaP; // mV
300 double N_NaP;
301
302 double g_peak_KNa;
303 double E_rev_KNa; // mV
304 double tau_D_KNa; // ms
305
306 double g_peak_T;
307 double E_rev_T; // mV
308 double N_T;
309
310 double g_peak_h;
311 double E_rev_h; // mV
312
313 bool voltage_clamp;
314 };
315
316 // ----------------------------------------------------------------
317
318 /**
319 * State variables of the model.
320 */
321 public:
322 struct State_
323 {
324 // y_ = [V, theta, Synapses]
325 enum StateVecElems_
326 {
327 V_M = 0,
328 THETA,
329 DG_AMPA,
330 G_AMPA,
331 DG_NMDA_TIMECOURSE,
332 G_NMDA_TIMECOURSE,
333 DG_GABA_A,
334 G_GABA_A,
335 DG_GABA_B,
336 G_GABA_B, // DO NOT INSERT ANYTHING UP TO HERE, WILL MIX UP
337 // SPIKE DELIVERY
338 m_fast_NMDA,
339 m_slow_NMDA,
340 m_Ih,
341 D_IKNa,
342 m_IT,
343 h_IT,
344 STATE_VEC_SIZE
345 };
346
347 //! neuron state, must be C-array for GSL solver
348 double y_[ STATE_VEC_SIZE ];
349
350 /** Timer (counter) for spike-activated repolarizing potassium current.
351 * Neuron is absolutely refractory during this period.
352 */
353 long ref_steps_;
354
355 double I_NaP_; //!< Persistent Na current; member only to allow recording
356 double I_KNa_; //!< Depol act. K current; member only to allow recording
357 double I_T_; //!< Low-thresh Ca current; member only to allow recording
358 double I_h_; //!< Pacemaker current; member only to allow recording
359
360 State_( const ht_neuron&, const Parameters_& p );
361 State_( const State_& s );
362
363 State_& operator=( const State_& s );
364
365 ~State_();
366
367 void get( DictionaryDatum& ) const;
368 void set( const DictionaryDatum&, const ht_neuron&, Node* node );
369 };
370
371 private:
372 // These friend declarations must be precisely here.
373 friend class RecordablesMap< ht_neuron >;
374 friend class UniversalDataLogger< ht_neuron >;
375
376
377 // ----------------------------------------------------------------
378
379 /**
380 * Buffers of the model.
381 */
382 struct Buffers_
383 {
384 Buffers_( ht_neuron& );
385 Buffers_( const Buffers_&, ht_neuron& );
386
387 UniversalDataLogger< ht_neuron > logger_;
388
389 /** buffers and sums up incoming spikes/currents */
390 std::vector< RingBuffer > spike_inputs_;
391 RingBuffer currents_;
392
393 /** GSL ODE stuff */
394 gsl_odeiv_step* s_; //!< stepping function
395 gsl_odeiv_control* c_; //!< adaptive stepsize control function
396 gsl_odeiv_evolve* e_; //!< evolution function
397 gsl_odeiv_system sys_; //!< struct describing system
398
399 // Since IntergrationStep_ is initialized with step_, and the resolution
400 // cannot change after nodes have been created, it is safe to place both
401 // here.
402 double step_; //!< step size in ms
403 double integration_step_; //!< current integration time step, updated by GSL
404
405 /**
406 * Input current injected by CurrentEvent.
407 * This variable is used to transport the current applied into the
408 * _dynamics function computing the derivative of the state vector.
409 * It must be a part of Buffers_, since it is initialized once before
410 * the first simulation, but not modified before later Simulate calls.
411 */
412 double I_stim_;
413 };
414
415 // ----------------------------------------------------------------
416
417 /**
418 * Internal variables of the model.
419 */
420 struct Variables_
421 {
422 //! size of conductance steps for arriving spikes
423 std::vector< double > cond_steps_;
424
425 //! Duration of potassium current.
426 int PotassiumRefractoryCounts_;
427
428 //! Voltage at beginning of simulation, for clamping
429 double V_clamp_;
430 };
431
432
433 // readout functions, can use template for vector elements
434 template < State_::StateVecElems_ elem >
435 double
get_y_elem_()436 get_y_elem_() const
437 {
438 return S_.y_[ elem ];
439 }
440 double
get_I_NaP_()441 get_I_NaP_() const
442 {
443 return S_.I_NaP_;
444 }
445 double
get_I_KNa_()446 get_I_KNa_() const
447 {
448 return S_.I_KNa_;
449 }
450 double
get_I_T_()451 get_I_T_() const
452 {
453 return S_.I_T_;
454 }
455 double
get_I_h_()456 get_I_h_() const
457 {
458 return S_.I_h_;
459 }
460
461 double get_g_NMDA_() const;
462
463 /**
464 * NMDA activation for given parameters
465 * Needs to take parameter values explicitly since it is called from
466 * _dynamics.
467 */
468 double m_NMDA_( double V, double m_eq, double m_fast, double m_slow ) const;
469
470 /**
471 * Return equilibrium value of I_h activation
472 *
473 * @param V Membrane potential for which to evaluate
474 * (may differ from y_[V_M] when clamping)
475 */
476 double m_eq_h_( double V ) const;
477
478 /**
479 * Return equilibrium value of I_T activation
480 *
481 * @param V Membrane potential for which to evaluate
482 * (may differ from y_[V_M] when clamping)
483 */
484 double m_eq_T_( double V ) const;
485
486 /**
487 * Return equilibrium value of I_T inactivation
488 *
489 * @param V Membrane potential for which to evaluate
490 * (may differ from y_[V_M] when clamping)
491 */
492 double h_eq_T_( double V ) const;
493
494 /**
495 * Return steady-state magnesium unblock ratio.
496 *
497 * Receives V_m as argument since it is called from ht_neuron_dyamics
498 * with temporary state values.
499 */
500 double m_eq_NMDA_( double V ) const;
501
502 /**
503 * Steady-state "D" value for given voltage.
504 */
505 double D_eq_KNa_( double V ) const;
506
507 static RecordablesMap< ht_neuron > recordablesMap_;
508
509 Parameters_ P_;
510 State_ S_;
511 Variables_ V_;
512 Buffers_ B_;
513 };
514
515
516 inline port
send_test_event(Node & target,rport receptor_type,synindex,bool)517 ht_neuron::send_test_event( Node& target, rport receptor_type, synindex, bool )
518 {
519 SpikeEvent e;
520 e.set_sender( *this );
521
522 return target.handles_test_event( e, receptor_type );
523 }
524
525
526 inline port
handles_test_event(SpikeEvent &,rport receptor_type)527 ht_neuron::handles_test_event( SpikeEvent&, rport receptor_type )
528 {
529 assert( B_.spike_inputs_.size() == 4 );
530
531 if ( not( INF_SPIKE_RECEPTOR < receptor_type && receptor_type < SUP_SPIKE_RECEPTOR ) )
532 {
533 throw UnknownReceptorType( receptor_type, get_name() );
534 return 0;
535 }
536 else
537 {
538 return receptor_type - 1;
539 }
540
541
542 /*
543 if (receptor_type != 0)
544 {
545 throw UnknownReceptorType(receptor_type, get_name());
546 }
547 return 0;*/
548 }
549
550 inline port
handles_test_event(CurrentEvent &,rport receptor_type)551 ht_neuron::handles_test_event( CurrentEvent&, rport receptor_type )
552 {
553 if ( receptor_type != 0 )
554 {
555 throw UnknownReceptorType( receptor_type, get_name() );
556 }
557 return 0;
558 }
559
560 inline port
handles_test_event(DataLoggingRequest & dlr,rport receptor_type)561 ht_neuron::handles_test_event( DataLoggingRequest& dlr, rport receptor_type )
562 {
563 if ( receptor_type != 0 )
564 {
565 throw UnknownReceptorType( receptor_type, get_name() );
566 }
567 return B_.logger_.connect_logging_device( dlr, recordablesMap_ );
568 }
569 }
570
571 #endif // HAVE_GSL
572 #endif // HT_NEURON_H
573