WATCH and "Init and Run"

duboismathieu · Post by **duboismathieu** » Thu Jan 10, 2008 11:08 am

Hello everybody using NEURON

I have a long lasting problem that I would like to solve.

I want to implement a conductance-based LIF. The model works this way:
-when the neuron receive an event, I increase instantaneously the value of conductance which then decays exponentially
-the membrane potential m obeys the classical conductance-based LIF

To solve the equation I use a BREAKPOINT block (lines 105-114).

The NET_RECEIVE block deals with refractory period and conductance updates etc...

I use WATCH to detect the threshold.

Here is the code:

Code: Select all


TITLE gIF4cd m version exc/inh
COMMENT
	Clock-based implementation of the gIF4 modeL with variable m.
	Use 2 types of synapses (excitatory and inhibitory).
	Default values of parameters are the same than in 'Simulation of networks of Spiking Neurons' Brette et al. 2006.
	During the refractory period, we don't do increase the conductance (simply ignore the incoming events).
ENDCOMMENT

COMMENT
	General constants
ENDCOMMENT
DEFINE TRUE 1
DEFINE FALSE 0

COMMENT
	Constants to control output flow
ENDCOMMENT
DEFINE NO_INFO 0
DEFINE EMIT_SPIKE 1
DEFINE SPIKE_SYNAPSE 2
DEFINE INTEGRATION 3
DEFINE DEBUG 4

COMMENT
	Types of events
ENDCOMMENT
DEFINE SYNAPTIC_EVENT 0
DEFINE SPIKE 1
DEFINE END_REFRAC 2
DEFINE INIT 3

NEURON {
	POINT_PROCESS gIF4cdPR
	RANGE m, refractory, refrac, thresh
	RANGE mL, mSe, mSi, itauSe, itauSi
	RANGE itaumL, iDtaumSe, iDtaumSi, itaumSe, itaumSi, itaum, C
	RANGE ID, verbose
}

STATE {
	m	(1)	://The membrane variable
}

PARAMETER {
	:Leak
	mL = 0			(1)	://Leak reversal potential
	itaumL = 0.05	  	(1/ms)  ://(inverse) Time constant at rest taumL

	:Synaptic component of time-constant
	mSe = 6			(1)	://Excitatory Synaptic reversal potential
	mSi = -2		(1)	://Inhibitory Synaptic reversal potential

	iDtaumSe = 0.03		(1/ms)	://Increase in itaumSe
	iDtaumSi = 0.335	(1/ms)	://Increase in itaumSi

	:Synaptic input
	itauSe = 0.2   		(1/ms)	://Decrease of the excitatory synaptic event
	itauSi = 0.1		(1/ms)	://Decrease of the inhibitory synaptic event

	:Neuron
	C = 200e-6		(microfarads)

	refrac = 5 		(ms)	://Duration of the refractory period
	thresh = 1		(1)	://Threshold of spike

	ID = 0			(1)	://ID
	verbose = 0		(1)	://Level of verbosity
}

ASSIGNED{
	://Time constants
	itaum		(1/ms)	://Total membrane time constant at time t
	itaumSe		(1/ms)	://Excitatory Synaptic component of membrane time constant at time t
	itaumSi		(1/ms)	://Inhibitory Synaptic component of membrane time constant at time t

	refractory	(1)	://If TRUE the neuron is in refractory period

	itaumSe0	(1/ms)	://Value of itaumSe at the arrival of a spike
	itaumSi0	(1/ms)	://Value of itaumSi at the arrival of a spike

	ts		(ms)	://Time of arrival of the spike
}

INITIAL {
	IF (verbose>=DEBUG) {printf("gIF4cd %g: Initialization\n", ID)}
	:Neuron
	refractory = FALSE
	t0 = 0
	m = mL

	:Time-constants
	itaumSe = 0
	itaumSe0 = itaumSe

	itaumSi = 0
	itaumSi0 = itaumSi

	itaum = itaumL + itaumSe + itaumSi

: 	printf("t=%g\n", t)
	net_send(0, INIT)
	IF (verbose>=DEBUG) {printf("gIF4cd %g: End of initialization\n", ID)}
}

BREAKPOINT {
	update_conductance()
	IF (refractory == FALSE) {
		SOLVE states METHOD cnexp
	} ELSE IF (fabs(t - ts) < .5) {
			m = 2
		} else {
			m = mL
		}
}

DERIVATIVE states{
	m' = -itaumL*(m-mL) -itaumSe*(m-mSe) -itaumSi*(m-mSi)
}

NET_RECEIVE(w) {

	IF (verbose >= SPIKE_SYNAPSE) {printf("*** gIF4cd %g: Net_Receive t=%g, m=%g itaum=%g, itaumSe=%g, itaumSi=%g, t0=%g, refractory=%g***\n", ID, t, m, itaum, itaumSe, itaumSi, t0, refractory)}

	IF (flag == SYNAPTIC_EVENT) { ://We have received a spike
		IF (verbose >= INTEGRATION) {printf("\tReceive a spike at t=%g, w=%g\n", t, w)}
		IF (verbose >= INTEGRATION) {printf("\tIntegration of the spike\n")}
		IF (verbose >= DEBUG) {
			printf("\tdriving forces: rest=%g, excitatory=%g, inhibitory=%g\n", (m-mL), (m-mSe), (m-mSi))
			printf("\ttime constants: rest=%g, excitatory=%g, inhibitory=%g\n", itaumL, itaumSe, itaumSi)
			printf("\tCurrent=%g\n", (itaumL*(m-mL) + itaumSe*(m-mSe) + itaumSi*(m-mSi))*C)
		}
		IF (refractory == TRUE) {
			IF (verbose >= INTEGRATION) {printf("\tRefactory period\n")}
		} ELSE {
			://Update the value of itau_m^Se and itau_m^Si just before the spike (eq. 2.8) and the value of itaum
			update_conductance()
			t0 = t
			process_event(w)
		}
	}
	if (flag == SPIKE) {
		IF (verbose >= EMIT_SPIKE) {printf("*** gIF4cd %g: Emit a spike at t=%g***\n", ID, t)}
		ts = t
		refractory = TRUE
		m = 2
		net_event(t)
		net_send(refrac, END_REFRAC)
	}
	if (flag == END_REFRAC) {
		m = mL
		refractory = FALSE
		IF (verbose >= SPIKE) {printf("\tEnd of refractory period at t=%g\n", t)}
	}
	if (flag == INIT) {
		IF (verbose >= DEBUG) {printf("\tInit\n")}
		WATCH (m > thresh) SPIKE
	}
	IF (verbose >= DEBUG) {printf("*** gIF4cd: End of Net_Receive ***\n")}
	IF (verbose >= SPIKE) {printf("\n")}
}

://
:// proc. that compute itaumSe and itaum at any time t
://
PROCEDURE update_conductance() {
	itaumSe = itaumSe0*exp(-(t-t0)*itauSe)
	itaumSi = itaumSi0*exp(-(t-t0)*itauSi)
	itaum = itaumSe + itaumSi + itaumL
}

://
:// proc. that update the synaptic time constants when the spike arrive
://
PROCEDURE process_event(w) {
	IF (w>0) {
		IF (verbose >= INTEGRATION) {printf("\tThis is an excitatory spike\n")}
		itaumSe = itaumSe + iDtaumSe
		itaum = itaumSe + iDtaumSi + itaumL
		IF (verbose >= INTEGRATION) {printf("\titaumSe updated by iDtaumSe=%g: itaumSe=%g\n", iDtaumSe, itaumSe)}
	}
	IF (w<0) {
		IF (verbose >= INTEGRATION) {printf("\tThis is an inhibitory spike\n")}
		itaumSi = itaumSi + iDtaumSi
		itaum = itaumSe + iDtaumSi + itaumL
		IF (verbose >= INTEGRATION) {printf("\titaumSi updated by iDtaumSi=%g: itaumSi=%g\n", iDtaumSi,itaumSi)}
	}

	://This is needed to ensure the correct computation of both components of synaptic conductance with only one variable to remember the previous synaptic event
	itaumSe0 = itaumSe
	itaumSi0 = itaumSi
}

://Removed not important code

The INITIAL block sends an event with flag==INIT to set the WATCH as explained in "The NEURON book".
I have also some hoc files to send events and so on.

The model works fine but sometime when I click "Init and Run", it doesn't detect the threshold. Debug messages show that the INIT-event is received far after t=0....

For instance with 50 excitatory spikes (weight=0.03 and one event is send each
ms) it gives a spike around t=24.8ms for the first run.
After clicking a second time "Init and Run", no spike is detected and one can see that the INIT event is sent at t=49ms...

This is not really a bug since it works most of the time but it's very strange to me.

Did I make a mistake somewhere? Is that a subtle initialization procedure?
I can send my hoc code (it's straightforward).

I use NEURON 6.1.2.

Post by **ted** » Thu Jan 10, 2008 3:33 pm

This is not really a bug since it works most of the time

"works most of the time" means there's a bug. The question is where.

I can send my hoc code

The simplest implementation that reproduces the problem would be useful.

Post by **hines** » Mon Jan 14, 2008 7:00 pm

Change every occurrence of t0 to t00
or any other name and declare it ASSIGNED.

Since the beginning of the modl language
for SCoP in the 80's, t0 has been a reserved word and defines the intitial value of t. This makes absolutely no sense in NEURON and is an embarassment I forgot about.
Another example that has caused trouble over the years is j0 which is the bessel
function. Anyway, I'm going to eliminate
t0 as a reserved word in the next
distribution. It is ok to use t0 even now in
ARTIFICIAL_CELL models since when
there is no BREAKPOINT, t0 is not
automatically declared static and there
is no t=t0 statement placed at the beginning of the INITIAL block.

duboismathieu · Post by **duboismathieu** » Tue Jan 15, 2008 4:24 am

Thanks a lot, it works perfectly now!

Effectively t0 was not declared in my mod file. This should have worried me.

Thanks again.
Mathieu