COMMENT

THIS IS THE MORE EFFICIENT VERSION of EPSPPlas.
Since short-term plasticity, and the synaptic dynamics is the same
for all the presynaptic terminals of a neuron, we can calculate these
parameters presynaptically, and then pass them to each synapse model.

For computational efficiency it is obviously better to process the release
parameters in the presynaptic mechanism, rather than do the same calculations
in each synaptic mechanism.  The problem is the synaptic delay.  
To deal with this problem I have created a history (histR, histG) of the 
synaptic conductances, that are accessed by the synaptic mechanisms. The vector 
index accessed corresponds to the delay.
In the vector the first index (0) is always the current time step.
Thus if dt = 0.1 to implement a delay of 1ms you should setpointer:
setpointer syn.R_1, IN[0].soma.histR_ExIAF[10]
for a zero ms delay
setpointer syn.R_1, IN[0].soma.histR_ExIAF[0]



ENDCOMMENT


INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}

NEURON { 
:    NONSPECIFIC_CURRENT i
   :POINT_PROCESS EPlasSom
   SUFFIX EPlasSom
	RANGE C, lastrelease, lastspike, releaseat, Delay
	GLOBAL Cdur, Deadtime, terror			 
	GLOBAL Alpha_1, Beta_1                  
	RANGE  ampa, R0_1, R1_1, Rinf_1, Rtau_1 
	GLOBAL Alpha_2, Beta_2 
   GLOBAL DurNMDA, tauNMDA
   :SHORT-TERM PLASTICITY
   GLOBAL U, trec, tfac
   RANGE  R, u, RG 

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

STATE {
	nmda				: fraction of open NMDA channels
	R_2
}

PARAMETER {

   terror
	Cdur	= 1	  (ms)		: transmitter duration (rising phase) 
	Deadtime = 1  (ms)		: mimimum time between release events 
 	Delay = 1.5     (ms)		: T - Synaptic Delay

	Alpha_1 = 1.5	(/ms mM)	: AMPA forward (binding) rate 
	Beta_1 = 0.75	(/ms)		: AMPA backward (unbinding) rate 
	Alpha_2 = 0.25 (/ms mM) 	: NMDA forward (binding) rate 
	Beta_2 = 0.025 (/ms)		: NMDA backward (unbinding) rate
   DurNMDA = 0.4               : used in sigmoid R_G, small values long time-peak
   tauNMDA	=50                 : in ms
	tfac = 0.01     (ms)		: this value should be close to 0 for no facilitation 
	trec = 500		 (ms) 	: recovery from depression time constant 
	U = 0.5 						: percent of transmitter released on first pulse
}


ASSIGNED { 
 
	dt		 (ms)
	v		 (mV)		: postsynaptic voltage
	i 		 (nA)		: current = g*(v - Erev)
	C		 (mM)		: transmitter concentration 
 
	ampa
	R0_1				: open channels at start of release
	R1_1				: open channels at end of release
	Rinf_1			: steady state channels open
	Rtau_1	(ms)		: time constant of channel binding
	RG				: reflects binding of Transm -> G protein
		 
	lastrelease	(ms)		: time of last spike 
	lastspike	(ms) 
	releaseat

	R				: Releasable pool 
	u				: for running value of U 
    
	
}

INITIAL {
   terror = dt/10
	C = 0 
	ampa = 0
	lastrelease = -9e4 
	lastspike   = -9e4 
	releaseat   = -9e4 
 
	R = 1 
	u = U 
 
}

BREAKPOINT {
    SOLVE release
}

PROCEDURE release() { LOCAL q
    :will crash if user hasn't set pre with the connect statement
	:FIND OUT THERE WAS A SPIKE 
	q = (t - lastspike)		: time since last release ended 
						 
	if (q > Deadtime) {		: ready for another release? 
		if (v > 0) {		: spike occured? 
			lastspike = t 
			releaseat = t + Delay
		} 
	} 
	
	: CALCULATE RELEASE PARAMETERS 
	q = (t - lastrelease -Cdur)			: time since last spike with delay 
	if (q > Deadtime) {          : start release 
 
		if (t > releaseat - terror && t < releaseat + terror) { 
			lastrelease = t 
			u = u*(exptable(-q/tfac)) + U*(1-u*exptable(-q/tfac))
			R = R*(1 - u)*exptable(-q/trec) + 1 - exptable(-q/trec)		 
			C = R*u				: start new release, turn on 
			Rinf_1 = C*Alpha_1 / (C*Alpha_1 + Beta_1) 
			Rtau_1 = 1 / ((Alpha_1 * C) + Beta_1) 
			R0_1 = ampa
			R_2=(1-R_2)*0.5+R_2
		}	 
						 
	} else if (q < 0) {			: still releasing? 
		: do nothing 
	} else if (C > 0) {			: in dead time after release, turn off 
		C = 0. 
		R1_1 = ampa 
	} 
 	if (C > 0) {				: transmitter being released? 
	   ampa = Rinf_1 + (R0_1 - Rinf_1) * exptable (- (t - lastrelease) / Rtau_1) 
	} else {					: no release occuring 
  	   ampa = R1_1 * exptable (- Beta_1 * (t - (lastrelease + Cdur))) 
	} 

    SOLVE G_protein METHOD cnexp : see http://www.neuron.yale.edu/phpBB/viewtopic.php?f=28&t=592

	VERBATIM
	return 0;
	ENDVERBATIM
} 

DERIVATIVE G_protein { 							: ready for anotherelease?

    R_2'=-(R_2/tauNMDA)

	if (R_2<0.01) {
		RG = 0.0
	} else {
		RG = 1/( 1 + exptable(-((R_2)-DurNMDA)/0.05) ) 		: binding of T -> G	
	}
	nmda' = Alpha_2 * RG * (1-nmda) - Beta_2 * nmda

    :printf("D-------->%f\n",t)

}

FUNCTION exptable(x) { 
	TABLE  FROM -10 TO 10 WITH 2000
 
	if ((x > -10) && (x < 10)) {
		exptable = exp(x)
	} else {
		exptable = 0.
	}
}