How to save synaptic changes?

[behii]

Hello,
I use a network similar to network of Cutsuridis et al in his paper ‘’Encoding and
retrieval in a model of the hippocampal CA1 microcircuit’’ . There are
100 CA1 pyramidal cells in this network, each of which receives input
from 100 CA3 cells. Synaptic weights are modifiable through STDP
learning rule and a mod file is available for implementing STDP.

I am trying to find out how synaptic weights change (or synaptic conductance) and save these
changes. I wonder if I should adjust the mod file of STDP or add commands to the
hoc file? I truly appreciate if you would please help me to figure it
out.

[ted]

NetCon weights are visible to hoc--see the documentation of the NetCon class.
Depending on the NMODL code used to implement your synaptic mechanisms, synaptic conductance may or may not be visible to hoc, but even if visible to hoc it may not be the most convenient indicator of synaptic weight. It would be more convenient if the code is written in a way that exposes the variables that govern the degree of potentiation/depotentiation. If these are part of the NetCon's weight vector they might be useful for inferring the degree of potentiation/depotentiation at the time of the most recent synaptic activation.

In order to provide more specific comments, I would have to know the details of how plasticity was implepemted. If you are referring to code that exists somewhere in ModelDB, it might help if you mention the accession number and the specific file, or provide the URL of that file.

[behii]

Thank you for the explanation.

I use this model:
http://senselab.med.yale.edu/modeldb/Sh ... del=123815

And here is the NMODL code for implementing STDP:

Code: Select all

: STDP by Hines, changed to dual exponential (BPG 6-1-09)
: Modified by BPG 13-12-08
: Limited weights: max weight is wmax and min weight is wmin
: (initial weight is specified by netconn - usually set to wmin)
: Rhythmic GABAB suppresses conductance and promotes plasticity.
: When GABAB is low, conductance is high and plasticity is off.

NEURON {
	POINT_PROCESS STDPE2
	RANGE tau1, tau2, e, i, d, p, dtau, ptau, thresh, wmax, wmin
	RANGE g, gbdel, gblen, gbint, gscale
	NONSPECIFIC_CURRENT i
}

UNITS {
	(nA) = (nanoamp)
	(mV) = (millivolt)
	(uS) = (microsiemens)
}

PARAMETER {
	tau1=.1 (ms) <1e-9,1e9>
	tau2 = 10 (ms) <1e-9,1e9>
	e = 0	(mV)
	wmax = 0 (uS)
	wmin = 0 (uS)	: not used - use netconn weight instead (BPG)
	d = 0 <0,1>: depression factor (multiplicative to prevent < 0)
	p = 0.5 : potentiation factor (additive, non-saturating)
	dtau = 34 (ms) : depression effectiveness time constant
	ptau = 17 (ms) : Bi & Poo (1998, 2001)
	thresh = -20 (mV)	: postsynaptic voltage threshold
	gbdel = 100 (ms) <1e-9,1e9> : initial GABAB off interval (ms)
	gbint = 100 (ms) <1e-9,1e9> : GABAB off interval (ms)
	gblen = 100 (ms) <1e-9,1e9> : GABAB on length (ms)
	gscale = 0.1	: relative suppression by GABAB
}

ASSIGNED {
	v (mV)
	i (nA)
	tpost (ms)
	on
	g (uS)
	gs
	factor
}

STATE {
	C (uS)
	B (uS)
}

INITIAL {
	LOCAL tp
	if (tau1/tau2 > .9999) {
		tau1 = .9999*tau2
	}
	C = 0
	B = 0
	tp = (tau1*tau2)/(tau2 - tau1) * log(tau2/tau1)
	factor = -exp(-tp/tau1) + exp(-tp/tau2)
	factor = 1/factor
	gs=1
	on=0	: initially not plastic
	tpost = -1e9
	net_send(0, 1)
	net_send(gbdel, 3)	: initial GABAB off period
}

BREAKPOINT {
	SOLVE state METHOD cnexp
	g = B - C
	i = g*gs*(v - e)
}

DERIVATIVE state {
	C' = -C/tau1
	B' = -B/tau2
}

NET_RECEIVE(w (uS), A, tpre (ms)) {
	INITIAL { A = 0  tpre = -1e9 }
	if (flag == 0) { : presynaptic spike  (after last post so depress)
:		printf("entry flag=%g t=%g w=%g A=%g tpre=%g tpost=%g\n", flag, t, w, A, tpre, tpost)
:		g = g + w + A	: only for single exp (BPG)
		C = C + (w + A)*factor
		B = B + (w + A)*factor
		tpre = t
		if (on == 1) {
			A = A * (1 - d*exp((tpost - t)/dtau))
                        print "dep" , A
		}
	}else if (flag == 2 && on == 1) { : postsynaptic spike
:		printf("entry flag=%g t=%g tpost=%g\n", flag, t, tpost)
		tpost = t
		FOR_NETCONS(w1, A1, tp) { : also can hide NET_RECEIVE args
:			printf("entry FOR_NETCONS w1=%g A1=%g tp=%g\n", w1, A1, tp)
			A1 = A1 + (wmax-w1-A1)*p*exp((tp - t)/ptau)
                        print "pot" , A1
		}
	} else if (flag == 1) { : flag == 1 from INITIAL block
:		printf("entry flag=%g t=%g\n", flag, t)
		WATCH (v > thresh) 2
	}
	else if (flag == 3) { : plasticity control
		if (on == 0) { : start plasticity
			on = 1
			gs = gscale
			net_send(gblen, 3)
		}
		else { : end burst
			on = 0
			gs = 1
			net_send(gbint, 3)
		}
	}
}

In fact, I want to extract the 100x100 weight matrix between CA3 and CA1 cells after training the network.At the begining of training this matrix is developed through clipped Hebbian learning rule and during training, weights are modified through STDP learning rule.

I am trying to find out:

- how synaptic conductances change after given time of training.
- how to save conductance or potentiation degree, A, -as you suggested- that exist in a NMODL code and in NET_RECEIVE block with such a model.
- what is the synaptic conductance between tow given cells.
- how to save the time that each conductance or potentiation degree changes.

It is important for my research to extract each connection matrix entry. I am confused with that!

Thank you so much for your time.

[ted]

I am trying to find out:

- how synaptic conductances change after given time of training.
- how to save conductance or potentiation degree, A, -as you suggested- that exist in a NMODL code and in NET_RECEIVE block with such a model.
- what is the synaptic conductance between tow given cells.
- how to save the time that each conductance or potentiation degree changes.

The synaptic conductances are described by continuous differential equations. Each will change continuously with time. Once any synaptic conductance becomes nonzero, it will change at every time step thereafter. You really need to focus on the degree of potentiation of each synapse.

The potentiation variable of any synaptic connection implemented with STDPE2 will be the second element of the weight vector of the NetCon that makes that connection. In other words, every NetCon attached to an instance of STDPE2 will have a weight vector whose second element reflects the degree of potentiation. So if you have a NetCon called nc, its potentiation will be called
nc.weight[1]
(see documentation of the NetCon class in the Programmer's Reference).

From the NET_RECEIVE block of STDPE2 note that
(1) the second element is initialized to 0, and
(2) every time a new synaptic event is delivered to an instance of STDPE2, the second element of every NetCon that attaches to that instance is updated.
(3) the second element is NOT updated at any other time

This has the following consequences:

Suppose you have a NetCon called nc attached to an instance of STDPE2 called syn.
If you stop a simulation at some time t and examine the second element of nc, you will get the state of potentiation of that particular connection as of the last time the second element of nc was evaluated. If tinput is the last time that nc delivered an event to syn, and tspike is the last time that a spike occurred in the segment to which nc is attached, then the second element of nc will have the value it was assigned at max(tinput, tspike) where max(a, b) = a if a>b, b if b>a. If you want to know when that was, you'll have to compare tpre (the third element of nc.weight, i.e. nc.weight[2]) and tpost (time of most recent spike at the location where syn is attached). Currently, STDPE2 does not allow you to discover tpost with hoc statements because tpost is only an ASSIGNED; if that is important to you, you will want to include the statement
RANGE tpost
in STDPE2's NEURON block.

[behii]

Thank you for your reply. Your previous comments were very helpful.

I still have some problems! As I have mentioned, I need to specify potentiation degree of each synapse in a network of 100 CA3 cells, each of which connecting to 100 CA1 cells.
There is a procedure named “connect CA3” that connects CA3 cells to CA1 and other cells in the network. Here is the procedure:

Code: Select all

proc connectCA3() {local i, j, cp, gid  localobj src, syn, synN, nc, fc, rs, conns, rc
  cp = $2	// connection probability
  mcell_ran4_init(connect_random_low_start_)
  conns = new Vector(nCA3)  // connection weights
  rc = new Vector(nCA3)  // random physical connectivity
  ncslist = new List()
  //ncilist = new List()
  
  // inputs to PCs determined by weight matrix
  for i=0, cells.count-1 {	// loop over possible target cells
    gid = gidvec.x[i]	// id of cell
    if (gid >= iPC && gid < npcell+iPC) {	// appropriate target cell
    print "cell ", gid
    syn = cells.object(i).pre_list.object($3)	// AMPA synapse with STDP
    syn.wmax = CHWGT
    syn.wmin = CLWGT
    syn.d = STDPDFAC	// depression
    syn.p = STDPPFAC	// potentiation
    syn.gscale = AMPASUPP	// fraction of AMPA during storage
    syn.thresh = STDPTHRESH	// threshold for postsynaptic voltage detection
    syn.gbdel = STDPSTART
    syn.gbint = STDPINT
    syn.gblen = STDPLEN
    synN = cells.object(i).pre_list.object($4)	// NMDA synapse
    rs = ranlist.object(i)  // the corresponding RandomStream
    rs.start()
    rs.r.uniform(0, 1)  // return integer in range 0..1
    rc.setrand(rs.r)	// generate random connectivity
    // open connections file
    fc = new File($s1)
    fc.ropen()
    conns.scanf(fc, gid+1, nCA3)	// read incoming weights for cell gid
    fc.close()
    for j=0, nCA3-1 {
      // only connection if physical connection exists
      if (rc.x[j] <= cp) {
        //print "   src ", j
        src = cells.object(j+iCA3).stim
        // set up connection from source to target NMDA synapse
//        nc = pc.gid_connect(j+iCA3, synN)
        nc = new NetCon(src, synN)
        ncslist.append(nc)
        nc.delay = CDEL
        nc.weight = CNWGT	// NMDA weight same for all connections
        // high AMPA if weight is 1
        if (conns.x[j] == 1) {
          // set up connection from source to target
          //nc = pc.gid_connect(j+iCA3, syn)
          nc = new NetCon(src, syn)
          ncslist.append(nc)
          nc.delay = CDEL
          nc.weight = CHWGT
        } else {
          // set up connection from source to target
          //nc = pc.gid_connect(j+iCA3, syn)
          nc = new NetCon(src, syn)
          ncslist.append(nc)
          nc.delay = CDEL
          nc.weight = CLWGT	// unlearned weight
        }
      }
    }
    }
  }
}

There is also another procedure to make CA3 input. Now, I have added a procedure named “savepotdeg” to the hoc file to save potentiation degree as:

Code: Select all

objref recA, fa, 
strdef fna,
recA = new Vector()
strdef potdeg
potdeg="Results/pdeg"

proc savepotdeg() { local i, j,k  localobj src, syn, nc
  for i=0, 99 {
   syn = cells.object(i).pre_list.object(EM_CA3)
   for j=0, nCA3-1 {
     src = cells.object(j+iCA3).stim
     nc= new NetCon(src, syn)
     a=nc.weight[1]
     recA.record(&nc.weight[1])
     print "potentiation degree   " , nc.weight[1]
     sprint(fna,"%s_matrixA.dat", potdeg)
     fa = new File(fna)
     fa.wopen()
     print "recA.size   " , recA.size
     for k=0, recA.size-1 {
         fa.printf("%g\n", recA.x[k])
     }
     fa.close()
     }
     }
     }

where src and syn refer to pre and post synaptic cells of the synapse, respectively, that I want to save corresponding potentiation degree. But after running the hoc file, there is no data saved in “potdeg_A.dat”.

Is it the right procedure to save potentiation degree?

Is netcon.weight[1] a vector or a single value?

How should I define source and target of Netcon to reach a specified synaptic weight vector?

Where should this procedure (savepotdeg) be called in hoc file?

Thank you so much,

[ted]

Presumably the network architecture has already been built before proc savepotdeg() is called. Then savepotdeg() comes along and tries to set up a new bunch of connections between cells. The result will be a different network than the one you started with. This is not useful. Get rid of it.

From proc connectCA3(), it appears that there is a List, called ncslist, to which will be appended every NetCon that connects two cells.

Code: Select all

for i=0,ncslist.count()-1 {
  . . . statements . . .
}

will iterate over every element of this list and do whatever the "statements" are. For example, to discover the weight and potentiation for each NetCon in your network,

Code: Select all

for i=0,ncslist.count()-1 {
  print i, ncslist.o(i).weight[0], ncslist.o(i).weight[2]
}

[behii]

I tried to work with potentiation degree but, it is not good enough for my purpose. I really need to find out conductance. When I save “ncslist.o(i).weight[0]”, this generates the same initial conductances that was set in “connectCA3”. It seems that “ncslist.o(i).weight[0]” cannot track conductances' changes! Is it correct?

In your last post, you mentioned I cannot derive conductance because it changes continuously with time. How should I change mod file of STDP so that I can record conductance?

In addition, I found an object named “STDPE2-“ in “plot what” that has 100 components from 0 to 99, each of which has other objects like B,C,e,g, etc. I ploted some conductances in this way, but I don’t know each of these conductances are between which cells. Do you have any idea how can I know that? and why these objects are 100 instead of 100*100=10000, because there are 10000 synapses? Do you think if I can use these objects to find out conductances' changes? How should I adjust them to have 100*100 connection matrix?

Thank you very much again for your help.

[ted]

From proc connectCA3(), it appears that there is a List, called ncslist, to which will be appended every NetCon that connects two cells.

Code: Select all

for i=0,ncslist.count()-1 {
  . . . statements . . .
}

will iterate over every element of this list and do whatever the "statements" are. For example, to discover the weight and potentiation for each NetCon in your network,

Code: Select all

for i=0,ncslist.count()-1 {
  print i, ncslist.o(i).weight[0], ncslist.o(i).weight[2]
}

My mistake--the potentiation term is the 2nd element in the weight vector. The correct statement would be

Code: Select all

for i=0,ncslist.count()-1 {
  print i, ncslist.o(i).weight[0], ncslist.o(i).weight[1]
}

[ted]

“ncslist.o(i).weight[0]” cannot track conductances' changes!

True. There is no assignment statement that has w on the left hand side. w remains unchanged throughout the simulation. The effect of any synaptic input event is manifest by what it does to the state variables B and C. These statements
C = C + (w + A)*factor
B = B + (w + A)*factor
in the NET_RECEIVE block, the calculation of factor in the INITIAL block, and the statements
g = B - C
i = g*gs*(v - e)
in the BREAKPOINT block, make me suspect that an input event with weight w will cause a synaptic conductance transient with peak amplitude given by
(w + A)*gs

You should test this for yourself with a very simple toy model: a single compartment with the Hodgkin-Huxley mechanism hh, to which you have attached an IClamp. Set the IClamp to deliver a 0.1 ms current pulse at 1 ms, and adjust its amplitude so that it is large enough to elicit a spike. Next change the IClamp's amplitude to 0, and attach a an instance of the STDPE2 mechanism driven by a NetStim via a NetCon. Set the NetStim to generate a single event at 1 ms, and set the NetCon's delay to 1 ms and its weight to 0.001. Run a simulation and verify that the STDPE2's g starts out at 0 ms, begins to increase at 2 ms, and peaks at 0.001 uS.

Next increase the number of events delivered by the NetStim to 2, and set the interval between them to 10 ms. Check the peak amplitudes of the two synaptic conductance transients. See if these agree with what you would predict from the values of w and A (the [0] and [1] elements of the NetCon's weight vector) and gs.

Now restore the IClamp's amplitude to a nonzero value. Vary the time at which the current pulse is delivered and see what happens to the synaptic conductance transients. Compare with what you would expect based on the values of w and A.

How should I change mod file of STDP so that I can record conductance?

You don't want to do that. Conductance changes continuously throughout the simulation. You really want the peak ampitude of each condutance transient that is elicited by a synaptic input--that's much less data to have to deal with. And you can calculate that from w, A, and gs. You will only need to record the value of gs whenever an input arrives, and of each NetCon's A when it changes, i.e. when it delivers an input event. I see that the NMODL code does not expose gs to hoc. This is not good--gs may change in the course of a simulation. You will want to declare gs as a RANGE variable in the NEURON block.

In addition, I found an object named “STDPE2-“ in “plot what” that has 100 components from 0 to 99, each of which has other objects like B,C,e,g, etc.

Good for you. It can be very difficult to understand someone else's code (hard enough for me to understand my own code). What you found was 100 instances of the STDPE2 class. These are not 100 synapses. They are 100 synaptic mechanisms. They may be used to represent the effects of many more than 100 synapses.

"How?" you might well ask.

By superposition. Suppose you have a biological cell with two synapses of the same kind (ampaergic, gabaergic, whatever) that are anatomically close to each other, so that they are also "electrically" close to each other. Activate one and it produces a conductance change g1. Activate the other and it also produces a conductance change g2. The total conductance produced by both is the sum of their individual conductances. They both see the same local membrane potential v, and the ions that flow through their channels have the same reversal potential es. So the total current is
istotal = (g1 + g2)*(v - es)
Now, the time courses of g1 and g2 are governed by one or more differential equations (DEs) that have identical form--in this case,
B' = -B/tau
C' = -C/tau
Furthermore, these DEs are linear.

So in building a computational model of these two biological synapses, we can take advantage of superposition, and use a single set of DEs to represent their response to synaptic activation. And if there were 100 synapses of the same type attached to an anatomically small region of a cell, and they were described by a set of linear DEs, we could represent all 100 of them with a single set of DEs.

So there may be only 100 instances of some synaptic mechanism X, but each instance may receive events from as many presynaptic spike sources as you like.

I ploted some conductances in this way

Well, the actual conductance associated with any STDPE2 is actually g*gs. In its present form STDPE2 only allows you to know the value of g, not the value of gs. gs may be either 1 or gscale, depending on what's happening in the network. Furthermore, the value of g represents the sum of all inputs that converge onto the STDPE2. So you are completely in the dark about the conductance produced by any individual synaptic connection between two cells. You can only discover the peak conductance change elicited by an individual synaptic activation, and you can only do that if you know gs, w, and A (the former is a property of the STDPE2 and the latter two belong to a particular NetCon).

but I don’t know each of these conductances are between which cells. Do you have any idea how can I know that?

You need to know which NetCon belongs to which presynaptic-postsynaptic cell pair. Keeping track of what is connected to what is something that is typically done during model setup, although sometimes for diagnostic purposes there is code that traces connections after setup has been completed (you may find it interesting to look at the code for the second model described in
Hines, M.L. and Carnevale, N.T.
Translating network models to parallel hardware in NEURON.
J. Neurosci. Methods 169:425-455, 2008
). I wouldn't be surprised if the code in ModelDB entry 123815 includes one or more procedures that take care of that.

[behii]

Thank you for your detailed explanations.

I built a toy model. There is a little difference between what is expected and what the results from STDPE2.g show.


Firstly,

[“Set the NetStim to generate a single event at 1 ms, and set the NetCon's delay to 1 ms and its weight to 0.001. Run a simulation and verify that the STDPE2's g starts out at 0 ms, begins to increase at 2 ms, and peaks at 0.001 uS.”*]

the exact value of peak is 0.000999953 instead of 0.001!


Secondly, I found in every step of computing g, gnew is near g´old + w +A. and g´old is gold * gs. When I computed g in this way, again, there is a little difference between what I computed from w and A and old g, and what STDPE2.g shows. For example, I calculated g to be 0.00175764 but STDPE2.g is equal to 0.00166772.



What is the cause of these differences?

And how should I estimate these differences to get true values of g by A,W and gs?