mod files with assert in verbatim failed under mswin

NMODL and the Channel Builder.
Post Reply
tom_morse
Posts: 43
Joined: Wed May 18, 2005 10:23 pm
Location: Yale University School of Medicine
Contact:

mod files with assert in verbatim failed under mswin

Post by tom_morse » Tue Aug 28, 2007 10:31 am

We received a request for help from a modeler who cound not get a NEURON simulation with mod files with assert statements such as this nmda.mod to work under mswin :

Code: Select all

COMMENT
-----------------------------------------------------------------------------
Simple synaptic mechanism derived for first order kinetics of
binding of transmitter to postsynaptic receptors.

A. Destexhe & Z. Mainen, The Salk Institute, March 12, 1993.

$Log: rnmda.mod,v $
Revision 1.6  1997/11/28 20:11:34  karchie
Phase calculation was still incorrect.  Fixed it (hopefully) for good.

Revision 1.5  1997/11/28 19:54:51  karchie
The INITIAL block was using phase incorrectly: phase 0+e now means
that we just fired, and phase 1-e that we're imminently firing.

Revision 1.4  1996/06/23 18:26:12  karchie
Negative phase values now indicate that spike arrivals are a Poisson process.

Revision 1.3  1996/05/07  23:32:20  karchie
Set R1 in INITIAL block to prevent weirdness just after t=0 when phase
is near 1.0.

Revision 1.2  1996/05/07  22:22:20  karchie
Modified to fire regularly with specified period and phase.
Removed pointer trigger mechanism.

-----------------------------------------------------------------------------


During the arrival of the presynaptic spike (detected by threshold 
crossing), it is assumed that there is a brief pulse (duration=Cdur)
of neurotransmitter C in the synaptic cleft (the maximal concentration
of C is Cmax).  Then, C is assumed to bind to a receptor Rc according 
to the following first-order kinetic scheme:

Rc + C ---(Alpha)--> Ro							(1)
       <--(Beta)--- 

where Rc and Ro are respectively the closed and open form of the 
postsynaptic receptor, Alpha and Beta are the forward and backward
rate constants.  If R represents the fraction of open gates Ro, 
then one can write the following kinetic equation:

dR/dt = Alpha * C * (1-R) - Beta * R					(2)

and the postsynaptic current is given by:

Isyn = gmax * R * (V-Erev)						(3)

where V is the postsynaptic potential, gmax is the maximal conductance 
of the synapse and Erev is the reversal potential.

If C is assumed to occur as a pulse in the synaptic cleft, such as

C     _____ . . . . . . Cmax
      |   |
 _____|   |______ . . . 0 
     t0   t1

then one can solve the kinetic equation exactly, instead of solving
one differential equation for the state variable and for each synapse, 
which would be greatly time consuming...  

Equation (2) can be solved as follows:

1. during the pulse (from t=t0 to t=t1), C = Cmax, which gives:

   R(t-t0) = Rinf + [ R(t0) - Rinf ] * exp (- (t-t0) / Rtau )		(4)

where 
   Rinf = Alpha * Cmax / (Alpha * Cmax + Beta) 
and
   Rtau = 1 / (Alpha * Cmax + Beta)

2. after the pulse (t>t1), C = 0, and one can write:

   R(t-t1) = R(t1) * exp (- Beta * (t-t1) )				(5)

The variables "phase" and "mean_ia_time" control presynaptic spike arrival:
  * phase, if in [0, 1), determines the time of the first of the regular
    arrivals.  If phase < 0, arrivals are a Poisson process.
  * mean_ia_time > 0 is the mean interspike interval.  If arrivals are
    regular, mean_ia_time > Cdur + Deadtime.

ENDCOMMENT

VERBATIM
#include <assert.h>

static char rcsid[] = "$Id: rnmda.mod,v 1.6 1997/11/28 20:11:34 karchie Exp $";

ENDVERBATIM


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

NEURON {
	POINT_PROCESS rNMDA
	RANGE C, R, R0, R1, g, gmax, lastrelease, phase, mean_ia_time, Alpha, Beta
	NONSPECIFIC_CURRENT i
	GLOBAL Cmax, Cdur, Erev, Deadtime, Rinf, Rtaur
}
UNITS {
	(nA) = (nanoamp)
	(mV) = (millivolt)
	(umho) = (micromho)
	(mM) = (milli/liter)
}

PARAMETER {

	Cmax	= 1	(mM)		: max transmitter concentration
	Cdur	= 25	(ms)		: transmitter duration (rising phase)
	Alpha	= .20333(/ms mM)	: forward (binding) rate
	Beta	= .00746(/ms)		: backward (unbinding) rate
	Erev	= 0	(mV)		: reversal potential
	Deadtime = 0	(ms)		: mimimum time between release events
	gmax		(umho)		: maximum conductance
	eta     = 0.33  (/mM)
	mag     = 1     (mM)
	gamma   = 0.06  (/mV)
	steep   = 8
	vhalf   = -20

	phase	= 0			: event time offset / poisson flag
	mean_ia_time = 100 (ms)		: event interarrival time
}

ASSIGNED {
	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 release
	R1				: open channels at end of release
	Rinf				: steady state channels open
	Rtaur		(ms)		: time constant of channel binding
	lastrelease	(ms)		: time of last spike
	ia_time		(ms)		: time between last and next spike
}

INITIAL {
	R = 0
	R1 = R
	C = 0
	Rinf = Cmax*Alpha / (Cmax*Alpha + Beta)
	Rtaur = (Alpha * Cmax) + Beta

	: Make sure mean IA time and phase are valid.
	: Valid phase values for regular firing are 0.0 <= phase < 1.0
	: If (phase < 0.0), events are a Poisson process.
VERBATIM
	assert(phase < 1.0);
	if (phase < 0.0) {
		assert(mean_ia_time > 0.0);
	} else {
		assert(mean_ia_time > Cdur + Deadtime);
	}
ENDVERBATIM

	: Set the first spike time.
	if (phase < 0.0) {
		: Events are a Poisson process.  We begin, on average,
                : halfway between events.
		ia_time = next_ia_time()
		lastrelease = t - ia_time / 2
		ia_time = ia_time / 2
	} else {
		: Use phase to set first spike time.
		ia_time = mean_ia_time
		lastrelease = t - phase * mean_ia_time
	}	
}

BREAKPOINT {
	SOLVE release
	g = (gmax * R)/(1 + eta * mag * exp((vhalf - v) / steep))
	i = g*(v - Erev)
}

PROCEDURE release() { LOCAL q
	: First determine whether we should start the next release,
	: end the previous release, or neither.
	q = t - lastrelease
	if (q >= ia_time) {			: time for the next release?
		C = Cmax
		R0 = R
		lastrelease = t
		ia_time = next_ia_time()	: calculate next ia time.
	} else if (q < Cdur) {			: still releasing?
		: do nothing
	} else if (C == Cmax) {			: in dead time after release
		R1 = R
		C = 0.
	}		
	
	if (C > 0) {				: transmitter being released?
	   R = Rinf + (R0 - Rinf) * exptable (- (t - lastrelease) * Rtaur)
	} else {				: no release occuring
  	   R = R1 * exptable (- Beta * (t - (lastrelease + Cdur)))
	}

	VERBATIM
	return 0;
	ENDVERBATIM
}


FUNCTION next_ia_time() {
	if (phase < 0.0) {
	} else {
		: Arrivals are regular
		next_ia_time = mean_ia_time
	}
}
FUNCTION exptable(x) { 
	TABLE  FROM -10 TO 10 WITH 2000
	if ((x > -10) && (x < 10)) {
		exptable = exp(x)
	} else {
		exptable = 0.
	}
}
A fix is to comment out all the assert statements.
The windows version of NEURON does not have the assert.h file
available for VERBATIM include statements.

Commenting out is OK as long as
the code using the mod files provides the proper values of the variables because
the assert macro is always only used for error checking. It would be better
coding practice to provide statements that check the values (to be true, i.e.non-zero), replacing the
assert statements with ordinary C language statements rather than omitting them
entirely as I have done.

Here is the mod file with the assert statements commented out:

Code: Select all

COMMENT
-----------------------------------------------------------------------------
Simple synaptic mechanism derived for first order kinetics of
binding of transmitter to postsynaptic receptors.

A. Destexhe & Z. Mainen, The Salk Institute, March 12, 1993.

$Log: rnmda.mod,v $
Revision 1.6  1997/11/28 20:11:34  karchie
Phase calculation was still incorrect.  Fixed it (hopefully) for good.

Revision 1.5  1997/11/28 19:54:51  karchie
The INITIAL block was using phase incorrectly: phase 0+e now means
that we just fired, and phase 1-e that we're imminently firing.

Revision 1.4  1996/06/23 18:26:12  karchie
Negative phase values now indicate that spike arrivals are a Poisson process.

Revision 1.3  1996/05/07  23:32:20  karchie
Set R1 in INITIAL block to prevent weirdness just after t=0 when phase
is near 1.0.

Revision 1.2  1996/05/07  22:22:20  karchie
Modified to fire regularly with specified period and phase.
Removed pointer trigger mechanism.

-----------------------------------------------------------------------------

During the arrival of the presynaptic spike (detected by threshold 
crossing), it is assumed that there is a brief pulse (duration=Cdur)
of neurotransmitter C in the synaptic cleft (the maximal concentration
of C is Cmax).  Then, C is assumed to bind to a receptor Rc according 
to the following first-order kinetic scheme:

Rc + C ---(Alpha)--> Ro							(1)
       <--(Beta)--- 

where Rc and Ro are respectively the closed and open form of the 
postsynaptic receptor, Alpha and Beta are the forward and backward
rate constants.  If R represents the fraction of open gates Ro, 
then one can write the following kinetic equation:

dR/dt = Alpha * C * (1-R) - Beta * R					(2)

and the postsynaptic current is given by:

Isyn = gmax * R * (V-Erev)						(3)

where V is the postsynaptic potential, gmax is the maximal conductance 
of the synapse and Erev is the reversal potential.

If C is assumed to occur as a pulse in the synaptic cleft, such as

C     _____ . . . . . . Cmax
      |   |
 _____|   |______ . . . 0 
     t0   t1

then one can solve the kinetic equation exactly, instead of solving
one differential equation for the state variable and for each synapse, 
which would be greatly time consuming...  

Equation (2) can be solved as follows:

1. during the pulse (from t=t0 to t=t1), C = Cmax, which gives:

   R(t-t0) = Rinf + [ R(t0) - Rinf ] * exp (- (t-t0) / Rtau )		(4)

where 
   Rinf = Alpha * Cmax / (Alpha * Cmax + Beta) 
and
   Rtau = 1 / (Alpha * Cmax + Beta)

2. after the pulse (t>t1), C = 0, and one can write:

   R(t-t1) = R(t1) * exp (- Beta * (t-t1) )				(5)

The variables "phase" and "mean_ia_time" control presynaptic spike arrival:
  * phase, if in [0, 1), determines the time of the first of the regular
    arrivals.  If phase < 0, arrivals are a Poisson process.
  * mean_ia_time > 0 is the mean interspike interval.  If arrivals are
    regular, mean_ia_time > Cdur + Deadtime.

tmm 20070823 this assert include removed from verbatim for windows compatibility:

#include <assert.h>

ENDCOMMENT

VERBATIM

static char rcsid[] = "$Id: rnmda.mod,v 1.6 1997/11/28 20:11:34 karchie Exp $";

ENDVERBATIM


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

NEURON {
	POINT_PROCESS rNMDA
	RANGE C, R, R0, R1, g, gmax, lastrelease, phase, mean_ia_time, Alpha, Beta
	NONSPECIFIC_CURRENT i
	GLOBAL Cmax, Cdur, Erev, Deadtime, Rinf, Rtaur
}
UNITS {
	(nA) = (nanoamp)
	(mV) = (millivolt)
	(umho) = (micromho)
	(mM) = (milli/liter)
}

PARAMETER {

	Cmax	= 1	(mM)		: max transmitter concentration
	Cdur	= 25	(ms)		: transmitter duration (rising phase)
	Alpha	= .20333(/ms mM)	: forward (binding) rate
	Beta	= .00746(/ms)		: backward (unbinding) rate
	Erev	= 0	(mV)		: reversal potential
	Deadtime = 0	(ms)		: mimimum time between release events
	gmax		(umho)		: maximum conductance
	eta     = 0.33  (/mM)
	mag     = 1     (mM)
	gamma   = 0.06  (/mV)
	steep   = 8
	vhalf   = -20

	phase	= 0			: event time offset / poisson flag
	mean_ia_time = 100 (ms)		: event interarrival time
}

ASSIGNED {
	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 release
	R1				: open channels at end of release
	Rinf				: steady state channels open
	Rtaur		(ms)		: time constant of channel binding
	lastrelease	(ms)		: time of last spike
	ia_time		(ms)		: time between last and next spike
}

INITIAL {
	R = 0
	R1 = R
	C = 0
	Rinf = Cmax*Alpha / (Cmax*Alpha + Beta)
	Rtaur = (Alpha * Cmax) + Beta

	: Make sure mean IA time and phase are valid.
	: Valid phase values for regular firing are 0.0 <= phase < 1.0
	: If (phase < 0.0), events are a Poisson process.

: tmm 20070823 these statements removed for windows compatibility VERBATIM
:	assert(phase < 1.0);
:	if (phase < 0.0) {
:		assert(mean_ia_time > 0.0);
:	} else {
:		assert(mean_ia_time > Cdur + Deadtime);
:	}
: ENDVERBATIM

	: Set the first spike time.
	if (phase < 0.0) {
		: Events are a Poisson process.  We begin, on average,
                : halfway between events.
		ia_time = next_ia_time()
		lastrelease = t - ia_time / 2
		ia_time = ia_time / 2
	} else {
		: Use phase to set first spike time.
		ia_time = mean_ia_time
		lastrelease = t - phase * mean_ia_time
	}	
}

BREAKPOINT {
	SOLVE release
	g = (gmax * R)/(1 + eta * mag * exp((vhalf - v) / steep))
	i = g*(v - Erev)
}

PROCEDURE release() { LOCAL q
	: First determine whether we should start the next release,
	: end the previous release, or neither.
	q = t - lastrelease
	if (q >= ia_time) {			: time for the next release?
		C = Cmax
		R0 = R
		lastrelease = t
		ia_time = next_ia_time()	: calculate next ia time.
	} else if (q < Cdur) {			: still releasing?
		: do nothing
	} else if (C == Cmax) {			: in dead time after release
		R1 = R
		C = 0.
	}		
	
	if (C > 0) {				: transmitter being released?
	   R = Rinf + (R0 - Rinf) * exptable (- (t - lastrelease) * Rtaur)
	} else {				: no release occuring
  	   R = R1 * exptable (- Beta * (t - (lastrelease + Cdur)))
	}

	VERBATIM
	return 0;
	ENDVERBATIM
}


FUNCTION next_ia_time() {
	if (phase < 0.0) {
	} else {
		: Arrivals are regular
		next_ia_time = mean_ia_time
	}
}


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

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


Post Reply