Spike-Time Dependent Plasticity in a parallel network

Moderator: wwlytton

Post Reply
Posts: 28
Joined: Tue Jan 13, 2015 11:25 am

Spike-Time Dependent Plasticity in a parallel network

Post by Grado » Thu Feb 05, 2015 5:05 pm

As the title suggests, I am trying to add spike-time dependent plasticity to a parallel network.

Basically what I need to do is modify the synaptic weights every time an action potential occurs based upon the relative time of the presynaptic and postsynaptic spikes of the spiking neuron in question.

My intuition tells me that this needs to be done in the NET_RECEIVE block of my synapse code, but I'm unsure as to how. Assume we have five neurons A1, A2, B, C1, and C2 connected in such that A1 and A2 synapse onto B and B onto C1 and C2, using NetCon objects A1B, A2B, BC1, and BC2. When neuron B spikes, I need to check the most recent spike times of all cells that synapse onto it as well as all cells it synapses onto. If A1 or A2 recently fired, their synaptic weights will be increased, while if C1 or C2 recently fired, their synaptic weights will be decreased.

Am I on the right track?

Code: Select all

TITLE simple AMPA receptors (discrete connections)


        Simple model for glutamate AMPA receptors


    Whole-cell recorded postsynaptic currents mediated by AMPA/Kainate
    receptors (Xiang et al., J. Neurophysiol. 71: 2552-2556, 1994) were used
    to estimate the parameters of the present model; the fit was performed
    using a simplex algorithm (see Destexhe et al., J. Computational Neurosci.
    1: 195-230, 1994).


    The simplified model was obtained from a detailed synaptic model that
    included the release of transmitter in adjacent terminals, its lateral
    diffusion and uptake, and its binding on postsynaptic receptors (Destexhe
    and Sejnowski, 1995).  Short pulses of transmitter with first-order
    kinetics were found to be the best fast alternative to represent the more
    detailed models.


    The first-order model can be solved analytically, leading to a very fast
    mechanism for simulating synapses, since no differential equation must be
    solved (see references below).


   Destexhe, A., Mainen, Z.F. and Sejnowski, T.J.  An efficient method for
   computing synaptic conductances based on a kinetic model of receptor binding
   Neural Computation 6: 10-14, 1994.

   Destexhe, A., Mainen, Z.F. and Sejnowski, T.J. Synthesis of models for
   excitable membranes, synaptic transmission and neuromodulation using a
   common kinetic formalism, Journal of Computational Neuroscience 1:
   195-230, 1994.

See also:




        RANGE C, R, R0, R1, g, Cmax
        GLOBAL Cdur, Alpha, Beta, Erev, blockTime, Rinf, Rtau

        (nA) = (nanoamp)
        (mV) = (millivolt)
        (umho) = (micromho)
        (mM) = (milli/liter)


        Cmax    = 0.5   (mM)            : max transmitter concentration (set = 1 to match ~/netcon/ampa.hoc)
        Cdur    = 0.3   (ms)            : transmitter duration (rising phase)
        Alpha   = 0.94  (/ms mM)        : forward (binding) rate
        Beta    = 0.18  (/ms)           : backward (unbinding) rate
        Erev    = 0     (mV)            : reversal potential
        blockTime = 2   (ms)            : time window following dbs event during which non-dbs events are blocked

        v               (mV)            : postsynaptic voltage
        i               (nA)            : current = g*(v - Erev)
        g               (umho)          : conductance
        C               (mM)            : transmitter concentration
        R                               : fraction of open channels
        R0                              : open channels at start of time period
        Rinf                            : steady state channels open
        Rtau            (ms)            : time constant of channel binding
        on                              : rising phase of PSC
        gmax                            : max conductance

        R = 0
        C = 0
        Rinf = Cmax*Alpha / (Cmax*Alpha + Beta)
        Rtau = 1 / ((Alpha * Cmax) + Beta)
        on = 0
        R0 = 0
        nspike = 0
        collisionBlock = 0

        SOLVE release
        i = R*(v - Erev)

PROCEDURE release() {
        if (on) {                               : transmitter being released?
           R = gmax*Rinf + (R0 - gmax*Rinf) * exptable (- (t - tLast) / Rtau)
        } else {                                : no release occuring
           R = R0 * exptable (- Beta * (t - tLast))
        return 0;

: following supports both saturation from single input and
: summation from multiple inputs
: if spike occurs during CDur then new off time is t + CDur
: ie. transmitter concatenates but does not summate
: Note: automatic initialization of all reference args to 0 except first

NET_RECEIVE(weight, ncType, ncPrb) {LOCAL ok, tmp :ncType 0=presyn cell, 1=dbs activated axon
    : flag is an implicit argument of NET_RECEIVE and  normally 0
    if (flag == 0) { : a spike, so turn on if not already in a Cdur pulse
        ok = 0
        if (ncType == 1) {
            collisionBlock = collisionBlock + 1
            net_send(blockTime, -1)
            ok = 1
        if (collisionBlock == 0) {
            ok = 1
        if (ok) {
            if (!on) {
                on = 1
                tLast = t
                R0 = R
                gmax = weight   :weight not additive from separate sources as in original ampa.mod
            nspike = nspike + 1 : come again in Cdur with flag = current value of nspike
            net_send(Cdur, nspike)
    if (flag == nspike) { : if this associated with last spike then turn off
        if (on) {
            on = 0
            tLast = t
            R0 = R
            gmax = 0
    if (flag == -1) {
        collisionBlock = collisionBlock - 1

FUNCTION exptable(x) {
        TABLE  FROM -10 TO 10 WITH 2000

        if ((x > -10) && (x < 10)) {
                exptable = exp(x)
        } else {
                exptable = 0.

Site Admin
Posts: 5053
Joined: Wed May 18, 2005 3:50 pm
Location: Yale University School of Medicine

Re: Spike-Time Dependent Plasticity in a parallel network

Post by ted » Fri Feb 06, 2015 7:50 pm

Your question was how to implement STDP, and the answer does involve the NET_RECEIVE block, but first here are a few comments about the model of a saturating synapse that you're starting with.

1. It is not compatible with adaptive integration ("cvode"). This is fixable, but I wouldn't bother with that if you only plan to use fixed time steps.

2. Each instance of this mechanism can handle only one input stream, that is, it can only be the target of a single NetCon. The comment just before the NET_RECEIVE block says

Code: Select all

  : following supports both saturation from single input and
  : summation from multiple inputs
but that's incorrect, and could be a trap for future users of this code who do not read and understand it in detail. Inside the nested conditional statement one finds

Code: Select all

  gmax = weight   :weight not additive from separate sources as in original ampa.mod
which is true.

Neither of these limitations applies to the saturating synaptic transmission mechanism AMPA_S that is presented in chapter 10 of The NEURON Book.

Next I should say something about the implementation of DBS-induced axonal block. It looks like you're doing something like this

Code: Select all

       nc0           dbsc0
pre0 -------> syn0 <------- dbs0

       nc1           dbsc1
pre1 -------> syn1 <------- dbs1

       nc2           dbsc2
pre2 -------> syn2 <------- dbs2
each pre is a presynaptic spike source
each nc is a NetCon
each syn is a separate synaptic instance
each dbs is a source of special events so that when syn receives a dbs event, it becomes temporarily unresponsive to new events from pre

Instead of implementing DBS-induced axonal block via NMODL code in the synaptic mechanism's mod file (which will be complicated enough without it), I'd suggest actually gating the spikes conveyed from each pre to the syn, before they even get to the syn. That is, the connection from one pre to the syn that it targets would be

Code: Select all

      nca         ncb
pre -------> ax -------> syn
dbs ---------+
where ax is an artificial spiking cell used as a gate that passes pre's events except for a user-specified time interval after it receives an event from dbs. nca's weight is 2 and its latency is 0, ncb's weight is whatever you want it to have, and its latency is whatever the total latency should be from pre to syn.
ncd's weight is -1 (any value < 0 would do) and its latency is 0.

The nice thing about this implementation of axonal block is that you don't have to anything to syn--syn can be any synaptic mechanism, even one that can handle multiple input streams.

I'm going to call the class that implements ax the EGate class ("event gate"). Here's the NMODL code for egate.mod

Code: Select all

: passes events from input to output
: until some time t' at which it receives an event with w < 0
: after which it blocks events until after t'+toff ms
: arrival of another event at time t'' with w < 0 starts a new "blocking interval"
: or prolongs the current one until t''+toff

  RANGE toff : how long it stays off after it is turned off
  RANGE on : 1 means passes events, 0 means blocks events

  toff = 5 (ms)


  on = 1 : default is to pass all events

  if (on == 1) {
    if (w>=0) {
    } else { : w<0, turn off
      on = 0
      net_send(toff, 1) : to turn back on
  } else { : off
    if (flag == 1) { : turn back on
      on = 1
    } else if (w<0) { : prolong block
I have put this mod file and a demo program here
Unzip the zip file, compile the mod file, then use NEURON to execute toynet.hoc.
Pre is a NetStim that generates an event stream at fixed intervals.
DBS is another NetStim that produces events that are used to emulate the effects of DBS.
EGate is an instance of the EGate mechanism.
The SpikePlot window shows the events generated by these three artificial spiking cells.
Click on Init & Run and see what happens.
Use the ArtCellGUI tool to change the EGate's toff, or the interval between events generated by DBS.

Now, what about STDP?

Here's an implementation of a fast activating, nonsaturating excitatory synapse (think AMPAergic synapse) that has STDP with multiplicative potentiation and depression.
https://www.neuron.yale.edu/ftp/ted/neu ... ynstdp.zip

The relationship between the tpre-tpost interval (tpre<tpost means the presynaptic spike precedes the postsynaptic spike) and the degree of potentiation or depression per spike pair follows the description by Bi and Poo in their 1998 and 2001 papers (see Fig. 1 in Bi2001).
expsynstdp.mod contains the NMODL code that defines the synaptic mechanism.
init.hoc generates a plot of synaptic potentiation or depression vs. pre-post latency.
rig.ses creates a RunControl, a couple of Graphs,
a VariableTimeStep panel,
and a Grapher. The Grapher plots the relationship between the pre-post latency and the
resulting degree of potentiation or depression induced by a single spike pair
which is specified by FUNCTION factor in expsynstdp.mod. Click on its "Plot" button to see the graph. You may have to use the graph's
View = plot
to rescale for best visualization.

"Well, nice, but I really want a synapse that saturates."

Doable, but the only way you're going to potentiate or depress such a synapse is by changing the duration of transmitter persistence in the synaptic cleft. That introduces a nonlinearity, the severity of which depends on whether the "standard 1 ms duration" is long enough for a significant fraction (say 30% or more) of channels to be activated. The nonlinearity would affect potentiation more than it affects depression. And because of that nonlinearity, the synaptic mechanism will only be able to handle a single input stream. Let me know how you want to proceed.

"Well, can I at least have a synapse with a conductance change that is biexponential, that is, described by two time constants?"

If that would be useful, I can create one based on Exp2Syn.

Posts: 28
Joined: Tue Jan 13, 2015 11:25 am

Re: Spike-Time Dependent Plasticity in a parallel network

Post by Grado » Fri May 22, 2015 10:26 am

I did end up getting STDP to work in my code, so if anyone comes here later and is interested in how I did it, message me and I'll be happy to help.

Post Reply