How to implement short term depression in AMPA/NMDA synapses

[maria diamantaki]

Hi all,

I want to implement the short term depression mechanism described in
Tsodyks M, Pawelzik K, Markram H (1998) Neural networks with dynamic synapses. Neural Comput 10:821-35
https://senselab.med.yale.edu/modeldb/s ... model=3815

Here is the mod file

Code: Select all

COMMENT
Revised 12/15/2000 in light of a personal communication 
from Misha Tsodyks that u is incremented _before_ x is 
converted to y--a point that was not clear in the paper.
If u is incremented _after_ x is converted to y, then 
the first synaptic activation after a long interval of 
silence will produce smaller and smaller postsynaptic 
effect as the length of the silent interval increases, 
eventually becoming vanishingly small.

Implementation of a model of short-term facilitation and depression 
based on the kinetics described in
  Tsodyks et al.
  Synchrony generation in recurrent networks 
  with frequency-dependent synapses
  Journal of Neuroscience 20:RC50:1-5, 2000.
Their mechanism represented synapses as current sources.
The mechanism implemented here uses a conductance change instead.

The basic scheme is

x -------> y    Instantaneous, spike triggered.
                Increment is u*x (see discussion of u below).
                x == fraction of "synaptic resources" that have 
                     "recovered" (fraction of xmtr pool that is 
                     ready for release, or fraction of postsynaptic 
                     channels that are ready to be opened, or some 
                     joint function of these two factors)
                y == fraction of "synaptic resources" that are in the 
                     "active state."  This is proportional to the 
                     number of channels that are open, or the 
                     fraction of max synaptic current that is 
                     being delivered. 
  tau_1
y -------> z    z == fraction of "synaptic resources" that are 
                     in the "inactive state"

  tau_rec
z -------> x

where x + y + z = 1

The active state y is multiplied by a synaptic weight to compute
the actual synaptic conductance (or current, in the original form 
of the model).

In addition, there is a "facilition" term u that 
governs the fraction of x that is converted to y 
on each synaptic activation.

  -------> u    Instantaneous, spike triggered, 
                happens _BEFORE_ x is converted to y.
                Increment is U*(1-u) where U and u both 
                lie in the range 0 - 1.
  tau_facil
u ------->      decay of facilitation

This implementation for NEURON offers the user a parameter 
u0 that has a default value of 0 but can be used to specify 
a nonzero initial value for u.

When tau_facil = 0, u is supposed to equal U.

Note that the synaptic conductance in this mechanism 
has the same kinetics as y, i.e. decays with time 
constant tau_1.

This mechanism can receive multiple streams of 
synaptic input via NetCon objects.  
Each stream keeps track of its own 
weight and activation history.

The printf() statements are for testing purposes only.
ENDCOMMENT


NEURON {
	POINT_PROCESS tmgsyn
	RANGE e, i
	RANGE tau_1, tau_rec, tau_facil, U, u0
	NONSPECIFIC_CURRENT i
}

UNITS {
	(nA) = (nanoamp)
	(mV) = (millivolt)
	(umho) = (micromho)
}

PARAMETER {
	: e = -90 mV for inhibitory synapses,
	:     0 mV for excitatory
	e = -90	(mV)
	: tau_1 was the same for inhibitory and excitatory synapses
	: in the models used by T et al.
	tau_1 = 3 (ms) < 1e-9, 1e9 >
	: tau_rec = 100 ms for inhibitory synapses,
	:           800 ms for excitatory
	tau_rec = 100 (ms) < 1e-9, 1e9 >
	: tau_facil = 1000 ms for inhibitory synapses,
	:             0 ms for excitatory
	tau_facil = 1000 (ms) < 0, 1e9 >
	: U = 0.04 for inhibitory synapses, 
	:     0.5 for excitatory
	: the (1) is needed for the < 0, 1 > to be effective
	:   in limiting the values of U and u0
	U = 0.04 (1) < 0, 1 >
	: initial value for the "facilitation variable"
	u0 = 0 (1) < 0, 1 >
}

ASSIGNED {
	v (mV)
	i (nA)
	x
}

STATE {
	g (umho)
}

INITIAL {
	g=0
}

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

DERIVATIVE state {
	g' = -g/tau_1
}

NET_RECEIVE(weight (umho), y, z, u, tsyn (ms)) {
INITIAL {
: these are in NET_RECEIVE to be per-stream
	y = 0
	z = 0
:	u = 0
	u = u0
	tsyn = t
: this header will appear once per stream
: printf("t\t t-tsyn\t y\t z\t u\t newu\t g\t dg\t newg\t newy\n")
}

	: first calculate z at event-
	:   based on prior y and z
	z = z*exp(-(t - tsyn)/tau_rec)
	z = z + ( y*(exp(-(t - tsyn)/tau_1) - exp(-(t - tsyn)/tau_rec)) / ((tau_1/tau_rec)-1) )
	: now calc y at event-
	y = y*exp(-(t - tsyn)/tau_1)

	x = 1-y-z

	: calc u at event--
	if (tau_facil > 0) {
		u = u*exp(-(t - tsyn)/tau_facil)
	} else {
		u = U
	}

: printf("%g\t%g\t%g\t%g\t%g", t, t-tsyn, y, z, u)

	if (tau_facil > 0) {
		state_discontinuity(u, u + U*(1-u))
	}

: printf("\t%g\t%g\t%g", u, g, weight*x*u)

	state_discontinuity(g, g + weight*x*u)
	state_discontinuity(y, y + x*u)

	tsyn = t

: printf("\t%g\t%g\n", g, y)
}

I was thinking to add that plasticity mechanism in the synaptic mechanisms I already have an AMPA and an NMDA in the NET RECEIVE block. Here is the AMPA mod file (the NMDA is almost the same with an Mg block).

Code: Select all

: AMPA mod

UNITS {
	(nA) 	= (nanoamp)
	(mV) 	= (millivolt)
	:(umho) = (micromho)
	(pS) 	= (picosiemens)
}

NEURON {
	POINT_PROCESS AMPA_S
	NONSPECIFIC_CURRENT i
	RANGE g, gmax, i
	GLOBAL Cdur, Alpha, Beta, Erev, Rinf, Rtau
}

PARAMETER {
	Cdur = 0.05 (ms)  	:1.0 (ms) :rising phase-transmitter duration
	Alpha = 0.2(/ms)  	:1.1 (/ms)
	Beta = 0.4(/ms)  	:0.19 (/ms)
	Erev = 0 (mV)
	gmax = 30(pS)  	:700(pS) 	:10e4  (pS)  
}

ASSIGNED {
	v 	(mV)
	i 	(nA)
	g 	(pS)
	Rtau 	(ms) 	: time const of channel binding
	Rinf 		: fraction of open channels at steady state
	synon 		: sum of weights of all synapses in the "onset" state
}

STATE {
	Ron Roff 	:total conductances of all synapses
} 	

INITIAL {
	synon = 0
	Rtau = 1 / (Alpha + Beta)
	Rinf = Alpha / (Alpha + Beta)
}

BREAKPOINT {
	SOLVE release METHOD cnexp
	g = (Ron + Roff)*gmax
	i = g*(v - Erev)
}

DERIVATIVE release {
	Ron' = (synon*Rinf - Ron)/Rtau
	Roff' = -Beta*Roff
}

NET_RECEIVE(weight, on, r0, t0 (ms)) {
	if (flag == 0) {
		:a spike arrived, start onset state if not already on
		if (!on) {
			synon = synon + weight
			r0 = r0*exp(-Beta*(t-t0)) 
			Ron = Ron + r0
			Roff = Roff - r0
			t0 = t
			on = 1
			net_send(Cdur, 1)
		} else {
			net_move(t + Cdur)
		}
	}

	if (flag == 1) {
		synon = synon - weight
		r0 = weight*Rinf + (r0-weight*Rinf)*exp(-(t-t0)/Rtau)
		Ron = Ron - r0
		Roff = Roff + r0
		t0 = t
		on = 0
	}
} 

So, I thought that the x, y, z parameters of Tsodyks code correspond to the r0, Ron, Roff of the ampa code and I tried to change the ampa mod file according to that, but the results I got were far from what I was expecting.. Probably I haven't understood correctly these codes.

Any hints on how I can implement that?

Thanks
maria

[ted]

I thought that the x, y, z parameters of Tsodyks code correspond to the r0, Ron, Roff of the ampa code

They're not remotely similar. Suggest you read
Tsodyks M, Pawelzik K, Markram H (1998) Neural networks with dynamic synapses. Neural Computation 10:821-35
and the discussion of how to model synaptic transmission with receptor saturation in
Destexhe, A.; Mainen, Z. & Sejnowski, T. (1994) An efficient method for computing synaptic conductances based on a kinetic model of receptor binding. Neural Computation 6:14-18

Any hints on how I can implement that?

The Tsodyks et al. model addresses transmitter release, and the Destexhe et al. model addresses receptor saturation. These are mediated by separate physical processes. If you want a synaptic model that includes both of these processes, both processes will have to be represented in the code.

[maria diamantaki]

Thanks for the reply Ted!

If you want a synaptic model that includes both of these processes, both processes will have to be represented in the code.

So it is possible to combine these two different processes, but I am sure it's not as simple as just adding some parts of one code to the other. How are they connected? I was trying to find that connection, when I ended up assuming (wrongly) that x, y, z could correspond to ro, Ron, Roff, mainly because in Tsodyks mechanism the synaptic conductance has the same kinetics as y, so I assumed (wrong again?) that the decaying time constant tau_1 is the same as Rtau.

Actually, all this is too confusing for me, since I just have a couple of months experience with NEURON. The point is that I want to add depression in the model and also try to keep the kinetics model of the synapses as a more biophysically realistic approach. What do you suggest to do? Is there a not too advanced approach I could follow?

In case I just keep Tsodyks model, how can I modify it for the NMDA synapse? Add the Mg block like in the kinetics model?

Any information and advice is more than welcome!
maria

[ted]

Before we go any further with this discussion, it should be noted that all conceptual models of biological phenomena are simplifications and abstractions.

So it is possible to combine these two different processes, but I am sure it's not as simple as just adding some parts of one code to the other. How are they connected?

You have to set up the connection yourself. What you propose is to combine a transmitter release mechanism with a receptor saturation mechanism. The transmitter release mechanism determines the amount of transmitter that is released into the synaptic cleft by a presynaptic spike. The receptor saturation mechanism assumes that the transmitter interaction with the receptor follows a simple first order reaction scheme.

in Tsodyks mechanism the synaptic conductance has the same kinetics as y

The Tsodyks et al. mechanism can be interpreted in several ways. One possible interpretation is:
y describes the time course of transmitter in the synaptic cleft
the interaction between transmitter and receptor is instantaneous, and produces a synaptic conductance that is proportional to y

This may or may not have been part of the thinking that led to the original conceptual model by Tsodyks et al., but it is an interpretation that allows one to start thinking about receptor saturation.

What you propose to do is to change the relationship between y and the synaptic conductance. You could continue to imagine that the transmitter-receptor interaction is instantaneous, but apply a squashing function to the relationship between synaptic conductance and y. Or you could assume that the transmitter-receptor interaction is described by a first order reaction with a saturable receptor. The choice is up to you, but it should be based on experimental evidence. Otherwise you're just adding arbitrary complications to something that is already artificial. And it is faulty logic to argue that mere complexity makes one model more scientifically useful (or somehow more "realistic") than another.

In case I just keep Tsodyks model, how can I modify it for the NMDA synapse? Add the Mg block like in the kinetics model?

That would work.

[maria diamantaki]

Hi Ted thanks for the feedback, it's really helpful. I'll try to figure out how i'll implement that and I'll come back if I have more questions. For sure I want to keep it simple, so I'll give it a try just with Tsodyks model.

Thanks a lot again!
maria

[ted]

As I think about it, I realize that the Tsodyks model has its own built-in limitation on the peak magnitude of synaptic conductance, because the sum x+y+z is limited to 1. Now, the most "biological" interpretation of y might be that y represents the amount of transmitter that is released from the presynaptic terminal by any particular presynaptic spike, but there's no way that y can be > 1 so that puts a limit on g. Also, the amount by which y increases is only a fraction of x, and the size of x is governed by how quickly z (the "inactive" or "consumed resource") decays. So there's already an upper bound on y, and the more rapidly the presynaptic terminal fires, the less y will be able to increase on each spike. A rough emulation of receptor saturation could be achieved by slowing the recovery of z to x, or reducing u (the fraction of x that is converted to y when a spike occurs).

[maria diamantaki]

Hi Ted,

Thanks for the info. I used Tsodyks model and it seems to work ok with my model. However, I want to run some simulations for different synaptic conductances so I changed the BREAKPOINT block, as below

Code: Select all

BREAKPOINT {
	SOLVE state METHOD cnexp
	g0 = g*gmax
	i = g0*(v - e)
}

in order to be able to evaluate the output (EPSPs) of my model for different synaptic strengths (gmax). This is actually working for low values of gmax ~ 30 pS (i get the depression), but not for higher, like 300 or above (it is like voltage reaches a peak value and doesn't go above it).

So, first of all, is this change correct? If not, can you suggest me any other way of simulating different conductances? I guess this problem I have is due to the fact you explain in your last post, the limitation of x+y+z = 1. How can I overcome that?

maria

[ted]

I want to run some simulations for different synaptic conductances so I changed the BREAKPOINT block, as below

Code: Select all

BREAKPOINT {
	SOLVE state METHOD cnexp
	g0 = g*gmax
	i = g0*(v - e)
}

Absolutely unnecessary. Get rid of g0 and use the NetCon's weight parameter to control synaptic strength.

This is actually working for low values of gmax ~ 30 pS (i get the depression), but not for higher, like 300 or above (it is like voltage reaches a peak value and doesn't go above it).

First, what is the "peak value." Provide a number please.