MOOSE - Multiscale Object Oriented Simulation Environment
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
GraupnerBrunel2012CaPlasticitySynHandler.cpp
Go to the documentation of this file.
1 /**********************************************************************
2 ** This program is part of 'MOOSE', the
3 ** Messaging Object Oriented Simulation Environment.
4 ** Copyright (C) 2011 Upinder S. Bhalla. and NCBS
5 ** It is made available under the terms of the
6 ** GNU Lesser General Public License version 2.1
7 ** See the file COPYING.LIB for the full notice.
8 **********************************************************************/
9 
10 #include "../basecode/header.h"
11 #include "Synapse.h"
12 #include "SynEvent.h" // only using the SynEvent class from this
13 #include "SynHandlerBase.h"
15 
16 #include <queue>
17 
19 {
20  static string doc[] =
21  {
22  "Name", "GraupnerBrunel2012CaPlasticitySynHandler",
23  "Author", "Aditya Gilra",
24  "Description",
25  "The GraupnerBrunel2012CaPlasticitySynHandler handles synapses"
26  "with Ca-based plasticity as per Higgins et al. 2014 and Graupner and Brunel 2012."
27  "Note 1:"
28  " Here, Ca ('chemical Ca') is updated only at each pre-spike, pre-spike+delayD and post-spike!"
29  " So it is inaccurate to use it for say Ca-dependent K channels in the electrical compartment,"
30  " for which you use are advised to use the CaPool i.e. 'electrical Ca'."
31  "Note 2:"
32  " Ca here is post-synaptic 'chemical Ca' common for all synapses in this SynHandler,"
33  " so weights of all pre-synapses connected to this SynHandler get updated"
34  " at each pre-spike, pre-spike+delayD and post-spike!"
35  " So if all pre-synaptic weights start out the same, they remain the same!!"
36  " If you want to consider each pre-synapse independently,"
37  " have independent SynHandlers for each synapse."
38  " If these SynHandlers are in the same electrical compartment,"
39  " you're essentially assuming these are on different spines,"
40  " with their own 'chemical Ca' which won't match the"
41  " 'electrical Ca' of the compartment (=dendrite)."
42  " If you put each SynHandler with a single synapse"
43  " in its own electrical compartment (=spine),"
44  " only then can you have an 'electrical Ca'"
45  " corresponding to the 'chemical Ca'."
46  "Three priority queues are used to manage pre, post, and pre+delayD spikes."
47  };
48 
50  "Ca",
51  "Ca is a post-synaptic decaying variable as a proxy for Ca concentration"
52  "and receives an impulse whenever a pre- or post- spike occurs."
53  "Caution: Ca is updated via an event-based rule, so it is only updated and valid"
54  "when a pre- or post- spike has occured, or at time delayD after a pre-spike."
55  "Do not use it to control a Ca dependent current, etc."
56  "See notes in the class Description: all pre-synapses get updated via the same post-synaptic Ca.",
59  );
60 
62  "CaInit",
63  "CaInit is the initial value for Ca",
66  );
67 
69  "tauCa",
70  "tauCa is the time constant for decay of Ca",
73  );
74 
76  "tauSyn",
77  "tauSyn is the time constant for synaptic weight evolution equation",
80  );
81 
83  "CaPre",
84  "CaPre is added to Ca on every pre-spike",
87  );
88 
90  "CaPost",
91  "CaPost is added to Ca on every post-spike",
94  );
95 
97  "delayD",
98  "Time delay D after pre-spike, when Ca is increased by Capre."
99  " delayD represents NMDA rise time.",
102  );
103 
105  "gammaP",
106  "gammaP is the potentiation factor for synaptic weight increase if Ca>thetaP",
109  );
110 
112  "gammaD",
113  "gammaD is the depression factor for synaptic weight decrease if Ca>thetaD",
116  );
117 
119  "thetaP",
120  "Potentiation threshold for Ca"
121  "User must ensure thetaP>thetaD, else simulation results will be wrong.",
124  );
125 
127  "thetaD",
128  "Depression threshold for Ca"
129  "User must ensure thetaP>thetaD, else simulation results will be wrong.",
132  );
133 
135  "bistable",
136  "If true, the synapse is bistable as in GraupnerBrunel2012 paper."
137  "The effect of potential on the weight update is usually ignorable"
138  " if Ca is above thetaP and thetaD most of the time.",
141  );
142 
144  "noisy",
145  "If true, turn noise on as per noiseSD",
148  );
149 
151  "noiseSD",
152  "Standard deviation of noise added to Ca",
155  );
156 
158  "weightMax",
159  "An upper bound on the weight",
162  );
163 
165  "weightMin",
166  "A lower bound on the weight",
169  );
170 
172  "weightScale",
173  "Scale all pre-synaptic weights by weightScale before adding to activation (default 1.0)"
174  "In the terminology of the paper Higgins et al 2012, weight is synaptic efficacy,"
175  "while weightScale*weight is what finally is added to activation variable.",
178  );
179 
180  static DestFinfo addPostSpike(
181  "addPostSpike",
182  "Handles arriving spike messages from post-synaptic neuron, inserts into postEvent queue.",
185 
187  "synapse",
188  "Sets up field Elements for synapse",
193  );
194 
195  static Finfo* GraupnerBrunel2012CaPlasticitySynHandlerFinfos[] =
196  {
197  &synFinfo, // FieldElement
198  &addPostSpike, // DestFinfo
199  &Ca, // Field
200  &CaInit, // Field
201  &tauCa, // Field
202  &tauSyn, // Field
203  &CaPre, // Field
204  &CaPost, // Field
205  &delayD, // Field
206  &thetaP, // Field
207  &thetaD, // Field
208  &gammaP, // Field
209  &gammaD, // Field
210  &weightMax, // Field
211  &weightMin, // Field
212  &weightScale, // Field
213  &noisy, // Field
214  &noiseSD, // Field
215  &bistable // Field
216  };
217 
219  static Cinfo synHandlerCinfo (
220  "GraupnerBrunel2012CaPlasticitySynHandler",
222  GraupnerBrunel2012CaPlasticitySynHandlerFinfos,
223  sizeof( GraupnerBrunel2012CaPlasticitySynHandlerFinfos ) / sizeof ( Finfo* ),
224  &dinfo,
225  doc,
226  sizeof( doc ) / sizeof( string )
227  );
228 
229  return &synHandlerCinfo;
230 }
231 
234 
236 {
237  Ca_ = 0.0;
238  CaInit_ = 0.0;
239  tauCa_ = 1.0;
240  tauSyn_ = 1.0;
241  CaPre_ = 0.0;
242  CaPost_ = 0.0;
243  thetaD_ = 0.0;
244  thetaP_ = 0.0;
245  gammaD_ = 0.0;
246  gammaP_ = 0.0;
247  delayD_ = 0.0;
248  weightMin_ = 0.0;
249  weightMax_ = 0.0;
250  weightScale_ = 1.0;
251  noisy_ = false;
252  noiseSD_ = 0.0;
253  bistable_ = true;
254  seed_ = 0;
256  reinitSeed();
257 
258  // std::cout << " Mean " << dist_.mean() << " std " << dist_.stddev() << std::endl;
259 }
260 
262 {
263 }
264 
266 {
267  if( 0 == seed_ )
269 
270  if( 0 == seed_ )
271  seed_ = rd_();
272 
273  MOOSE_DEBUG( "Seed is set to " << seed_ );
274  rng_.seed( seed_ );
275 }
276 
280  )
281 {
282  synapses_ = ssh.synapses_;
283  for ( vector< Synapse >::iterator
284  i = synapses_.begin(); i != synapses_.end(); ++i )
285  i->setHandler( this );
286 
287  // For no apparent reason, priority queues don't have a clear operation.
288  while( !events_.empty() )
289  events_.pop();
290  while( !delayDPreEvents_.empty() )
291  events_.pop();
292  while( !postEvents_.empty() )
293  postEvents_.pop();
294 
295  return *this;
296 }
297 
299 {
300  unsigned int prevSize = synapses_.size();
301  synapses_.resize( v );
302  for ( unsigned int i = prevSize; i < v; ++i )
303  synapses_[i].setHandler( this );
304 }
305 
307 {
308  return synapses_.size();
309 }
310 
312 {
313  static Synapse dummy;
314  if ( i < synapses_.size() )
315  return &synapses_[i];
316  cout << "Warning: GraupnerBrunel2012CaPlasticitySynHandler::getSynapse: index: " << i <<
317  " is out of range: " << synapses_.size() << endl;
318  return &dummy;
319 }
320 
322  unsigned int index, double time, double weight )
323 {
324  assert( index < synapses_.size() );
325  events_.push( PreSynEvent( index, time, weight ) );
326  delayDPreEvents_.push( PreSynEvent( index, time+delayD_, weight ) );
327 }
328 
330  unsigned int index ) const
331 {
332  if ( events_.empty() )
333  return 0.0;
334  return events_.top().time;
335 }
336 
337 
339 {
340  postEvents_.push( PostSynEvent( time ) );
341 }
342 
344 {
345  double CaOld = Ca_;
346  double deltaT = currTime-lastCaUpdateTime_;
347  // Update Ca, but CaPre and CaPost are added in vProcess()
348  Ca_ *= exp(-deltaT/tauCa_);
349  lastCaUpdateTime_ = currTime;
350  weightFactors wUp; // by default all are set to 0.0
351 
352  // calculate/approximate time spent above potentiation and depression thresholds
353  // see pg 13 of Higgins et al | October 2014 | Volume 10 | Issue 10 | e1003834 | PLOS Comp Biol
354  // starting from bottom condition, going upwards in the algorithm given in above paper
355  if (CaOld <= thetaD_)
356  {
357  }
358  else if (CaOld <= thetaP_)
359  {
360  //cout << "tD<Caold<tP" << "\n";
361  if (Ca_ <= thetaD_)
362  {
363  wUp.tD = tauCa_*log(CaOld/thetaD_);
364  //cout << "Ca<tD" << "\n";
365  }
366  else
367  {
368  wUp.tD = deltaT;
369  //cout << "Ca>tD" << "\n";
370  }
371  }
372  else
373  {
374  //cout << "Caold>tP" << "\n";
375  if (Ca_ <= thetaD_)
376  {
377  wUp.tP = tauCa_*log(CaOld/thetaP_);
378  wUp.tD = tauCa_*log(thetaP_/thetaD_);
379  //cout << "Ca<tD" << "\n";
380  //cout << "Caold = " << CaOld << "thetaP = " << thetaP_ << "\n";
381  }
382  else if (Ca_ <= thetaP_)
383  {
384  wUp.tP = tauCa_*log(CaOld/thetaP_);
385  wUp.tD = deltaT - wUp.tP;
386  //cout << "Ca<tP" << "\n";
387  }
388  else
389  {
390  wUp.tP = deltaT;
391  //cout << "Ca>tP" << "\n";
392  }
393  }
394  wUp.t0 = deltaT - wUp.tP - wUp.tD;
395 
396  // Depending on tP and tD, I return A,B,C factors for weight update
397  // (see page 13 of Higgins et al 2014).
398  // A,B,C,D,E are used to compute the weight change for one or multiple synapses
399  // A is weight independent, B,D are weight dependent and C,E are accumulated noise
400  if (wUp.tP > 0)
401  {
402  double gPgD = gammaP_+gammaD_;
403  wUp.A = gammaP_/gPgD*(1.0-exp(-wUp.tP*gPgD/tauSyn_));
404  wUp.B = exp(-wUp.tP*gPgD/tauSyn_);
405  if (noisy_)
406  {
407  wUp.C = noiseSD_ * dist_(rng_) *
408  sqrt( ( 1.0-exp(-2*gPgD*wUp.tP/tauSyn_) ) / gPgD );
409  // cout << " A = " << wUp.A << " B = " << wUp.B << " C = " << wUp.C << "\n";
410  }
411  else
412  {
413  wUp.C = 0.0;
414  }
415  }
416 
417  if (wUp.tD > 0)
418  {
419  wUp.D = exp(-wUp.tD*gammaD_/tauSyn_);
420  if (noisy_)
421  {
422  wUp.E = noiseSD_ * dist_( rng_ ) *
423  sqrt( ( 1.0-exp(-2*gammaD_*wUp.tD/tauSyn_) ) / 2.0/gammaD_ );
424  //cout << "D = " << wUp.D << " E = " << wUp.E << "\n";
425  }
426  else
427  {
428  wUp.E = 0.0;
429  }
430  }
431 
432  //cout << currTime << " tD = " << wUp.tD << " tP = " << wUp.tP << "\n";
433  // tP, tD, A, B, C, D, E of wUp
434  // are set to 0.0 by default in struct constructor
435  return wUp; // return by value, i.e. copies the struct, so wUp going out of scope doesn't matter
436  // malloc and returning pointer is more expensive for small structs
437  // see: http://stackoverflow.com/questions/9590827/is-it-safe-to-return-a-struct-in-c-c
438 }
439 
441 {
442  double newWeight = synPtr->getWeight();
443  //cout << " oldweight = " << newWeight << "\n";
444  if (wFacPtr->tP > 0.0)
445  {
446  newWeight = wFacPtr->A + wFacPtr->B*newWeight + wFacPtr->C;
447  //cout << " midweight = " << newWeight << "\n";
448  }
449  if (wFacPtr->tD > 0.0)
450  {
451  newWeight = wFacPtr->D*newWeight + wFacPtr->E; // update the weight again
452  //cout << " newweight = " << newWeight << "\n";
453  }
454 
455  // potential is usually ignorable when t0 is small,
456  // i.e. Ca is mostly above thetaD or thetaP
457  //cout << "before bistable newWeight = " << newWeight << "\n";
458  if (bistable_)
459  {
460  double chi0 = pow((newWeight-0.5),2.0) / (newWeight*(newWeight-1));
461  double weightDeviation = 0.5*sqrt( 1.0 + 1.0/(chi0*exp(wFacPtr->t0/2.0/tauSyn_)-1.0) );
462  if (newWeight<0.5)
463  {
464  newWeight = 0.5 - weightDeviation;
465  }
466  else
467  {
468  newWeight = 0.5 + weightDeviation;
469  }
470  }
471  //cout << "after bistable newWeight = " << newWeight << "\n";
472 
473  // clip weight within [weightMin,weightMax]
474  newWeight = std::max(weightMin_, std::min(newWeight, weightMax_));
475  //cout << " clipweight = " << newWeight << "\n";
476  //cout << " A = " << wFacPtr->A << " B = " << wFacPtr->B << " C = " << wFacPtr->C
477  // << " D = " << wFacPtr->D << " E = " << wFacPtr->E << " newW = "<< newWeight << "\n";
478  synPtr->setWeight( newWeight );
479 }
480 
482 {
483  double activation = 0.0;
484  double currTime = p->currTime;
485  bool CaFactorsUpdated = false; // Ca-decay and weight Factors updation is done only once
486  // if any pre-spike, pre-spike+delayD or post-spike occurs
487  // Ca is bumped by CaPre and CaPost on
488  // pre-spike+delayD and post-spike respectively
489  // At the end, if any event above occurred,
490  // change all pre-synaptic weights based on Ca
491  weightFactors wFacs;
492 
493  // process pre-synaptic spike events for activation, Ca and weight update
494  while( !events_.empty() && events_.top().time <= currTime )
495  {
496  PreSynEvent currEvent = events_.top();
497 
498  unsigned int synIndex = currEvent.synIndex;
499  // Warning, coder! 'STDPSynapse currSyn = synapses_[synIndex];' is wrong,
500  // it creates a new, shallow-copied object.
501  // We want only to point to the same object.
502  Synapse* currSynPtr = &synapses_[synIndex];
503 
504  // activate the synapse for every pre-spike
505  // If the synapse has a delay, the weight could be updated during the delay!
506  // currEvent.weight is the weight before the delay!
507  // Might be better to use currSynPtr->getWeight()
508  // or even the latest updated weight below?!
509  // Using currSynPtr->getWeight().
510  // Still, weight update is event driven and not continuous
511  // so using currSynPtr->getWeight() takes care of
512  // weight updates due to events/spikes during the delay
513  // but doesn't take care of continuous weight evolution during the delay.
514  // The weight update is done
515  // after sending the current weight to activation
516  // Maybe I should use the updated weight to send to activation?!
517 
518  // See: http://www.genesis-sim.org/GENESIS/Hyperdoc/Manual-26.html#synchan
519  // Send out weight / dt for every spike
520  // Since it is an impulse active only for one dt,
521  // need to send it divided by dt.
522  // Can connect activation to SynChan (double exp)
523  // or to LIF as an impulse to voltage.
524  //activation += currEvent.weight * weightScale_ / p->dt;
525  activation += currSynPtr->getWeight() * weightScale_ / p->dt;
526 
527  // update only once for this time-step if an event occurs
528  if (!CaFactorsUpdated)
529  {
530  // update Ca and weightFactors
531  wFacs = updateCaWeightFactors( currTime );
532  CaFactorsUpdated = true;
533  }
534 
535  events_.pop();
536  }
537  if ( activation != 0.0 )
538  SynHandlerBase::activationOut()->send( e, activation );
539 
540  // process delayed pre-synaptic spike events for Ca and weight update
541  // delayD after pre-spike accounts for NMDA rise time
542  while( !delayDPreEvents_.empty() && delayDPreEvents_.top().time <= currTime )
543  {
544  // Update Ca, and add CaPre
545  // update only once for this time-step if an event occurs
546  if (!CaFactorsUpdated)
547  {
548  // update Ca and weightFactors
549  wFacs = updateCaWeightFactors( currTime );
550  CaFactorsUpdated = true;
551  }
552  Ca_ += CaPre_;
553 
554  delayDPreEvents_.pop();
555  }
556 
557  // process post-synaptic spike events for Ca and weight update
558  while( !postEvents_.empty() && postEvents_.top().time <= currTime )
559  {
560  // update Ca, then add CaPost
561  // update only once for this time-step if an event occurs
562  if (!CaFactorsUpdated)
563  {
564  // update Ca and weightFactors
565  wFacs = updateCaWeightFactors( currTime );
566  CaFactorsUpdated = true;
567  }
568  Ca_ += CaPost_;
569 
570  postEvents_.pop();
571  }
572 
573  // If any event has happened, update all pre-synaptic weights
574  // If you want individual Ca for each pre-synapse
575  // create individual SynHandlers for each
576  if (CaFactorsUpdated)
577  {
578  // Change weight of all synapses
579  for (unsigned int i=0; i<synapses_.size(); i++)
580  {
581  // Warning, coder! 'Synapse currSyn = synapses_[i];' is wrong,
582  // it creates a new, shallow-copied object.
583  // We want only to point to the same object.
584  Synapse* currSynPtrHere = &synapses_[i];
585  updateWeight(currSynPtrHere,&wFacs);
586  }
587  }
588 
589 }
590 
592 {
593  // For no apparent reason, priority queues don't have a clear operation.
594  while( !events_.empty() )
595  events_.pop();
596  while( !delayDPreEvents_.empty() )
597  events_.pop();
598  while( !postEvents_.empty() )
599  postEvents_.pop();
600  Ca_ = CaInit_;
601 }
602 
604 {
605  unsigned int newSynIndex = synapses_.size();
606  synapses_.resize( newSynIndex + 1 );
607  synapses_[newSynIndex].setHandler( this );
608  return newSynIndex;
609 }
610 
611 
613 {
614  assert( msgLookup < synapses_.size() );
615  synapses_[msgLookup].setWeight( -1.0 );
616 }
617 
619 {
620  Ca_ = v;
621 }
622 
624 {
625  return Ca_;
626 }
627 
629 {
630  CaInit_ = v;
631 }
632 
634 {
635  return CaInit_;
636 }
637 
639 {
640  CaPre_ = v;
641 }
642 
644 {
645  return CaPre_;
646 }
647 
649 {
650  CaPost_ = v;
651 }
652 
654 {
655  return CaPost_;
656 }
657 
659 {
660  thetaP_ = v;
661 }
662 
664 {
665  return thetaP_;
666 }
667 
669 {
670  thetaD_ = v;
671 }
672 
674 {
675  return thetaD_;
676 }
677 
679 {
680  gammaD_ = v;
681 }
682 
684 {
685  return gammaD_;
686 }
687 
689 {
690  gammaP_ = v;
691 }
692 
694 {
695  return gammaP_;
696 }
697 
699 {
700  delayD_ = v;
701 }
702 
704 {
705  return delayD_;
706 }
707 
709 {
710  if ( rangeWarning( "tauCa", v ) ) return;
711  tauCa_ = v;
712 }
713 
715 {
716  return tauCa_;
717 }
718 
720 {
721  if ( rangeWarning( "tauSyn", v ) ) return;
722  tauSyn_ = v;
723 }
724 
726 {
727  return tauSyn_;
728 }
729 
731 {
732  noisy_ = v;
733 }
734 
736 {
737  return noisy_;
738 }
739 
741 {
742  noiseSD_ = v;
743 }
744 
746 {
747  return noiseSD_;
748 }
749 
751 {
752  bistable_ = v;
753 }
754 
756 {
757  return bistable_;
758 }
759 
761 {
762  weightMax_ = v;
763 }
764 
766 {
767  return weightMax_;
768 }
769 
771 {
772  weightMin_ = v;
773 }
774 
776 {
777  return weightMin_;
778 }
779 
781 {
782  weightScale_ = v;
783 }
784 
786 {
787  return weightScale_;
788 }
GraupnerBrunel2012CaPlasticitySynHandler & operator=(const GraupnerBrunel2012CaPlasticitySynHandler &other)
void setNumSynapses(unsigned int num)
void addSpike(unsigned int index, double time, double weight)
double currTime
Definition: ProcInfo.h:19
void setWeight(double v)
Definition: Synapse.cpp:78
Definition: Dinfo.h:60
Definition: EpFunc.h:64
static DestFinfo dummy("dummy","This Finfo is a dummy. If you are reading this you have used an invalid index", 0)
double getWeight() const
Definition: Synapse.cpp:88
unsigned int getNumSynapses() const
void log(string msg, serverity_level_ type=debug, bool redirectToConsole=true, bool removeTicks=true)
Log to console (and to a log-file)
static const Cinfo * GraupnerBrunel2012CaPlasticitySynHandlerCinfo
static const Cinfo * initCinfo()
void updateWeight(Synapse *synPtr, weightFactors *wFacPtr)
Synapse * getSynapse(unsigned int i)
priority_queue< PreSynEvent, vector< PreSynEvent >, CompareSynEvent > delayDPreEvents_
static const Cinfo * initCinfo()
Definition: Synapse.cpp:15
priority_queue< PreSynEvent, vector< PreSynEvent >, CompareSynEvent > events_
double dt
Definition: ProcInfo.h:18
Definition: Eref.h:26
int getGlobalSeed()
Definition: global.cpp:206
bool rangeWarning(const string &field, double value)
static SrcFinfo1< double > * activationOut()
unsigned int synIndex
Definition: SynEvent.h:52
priority_queue< PostSynEvent, vector< PostSynEvent >, ComparePostSynEvent > postEvents_
unsigned int addSynapse()
Adds a new synapse, returns its index.
static const Cinfo * synHandlerCinfo
Definition: Cinfo.h:18
Definition: Finfo.h:12