Problems with randomizing point process noise

[bscoventry]

Hi all,

I have a question regarding the initialization of the normrand function. I believe that the random number sequence is based only on the initial state of the random number generator seed. I am using a piece of code to generate synaptic noise in a neuron(modified from Destexhe's 2005 model in ModelDB).

Code: Select all

NEURON {
	POINT_PROCESS Gfluct
	RANGE g_e, g_i, E_e, E_i, g_e0, g_i0, g_e1, g_i1
	RANGE std_e, std_i, tau_e, tau_i, D_e, D_i
	RANGE new_seed
	NONSPECIFIC_CURRENT i
}

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

PARAMETER {
	dt		(ms)

	E_e	= 0 	(mV)	: reversal potential of excitatory conductance
	E_i	= -75 	(mV)	: reversal potential of inhibitory conductance

	g_e0	= 0.0121 (umho)	: average excitatory conductance
	g_i0	= 0.0573 (umho)	: average inhibitory conductance

	std_e	= 0.0030 (umho)	: standard dev of excitatory conductance
	std_i	= 0.0066 (umho)	: standard dev of inhibitory conductance

	tau_e	= 2.728	(ms)	: time constant of excitatory conductance
	tau_i	= 10.49	(ms)	: time constant of inhibitory conductance
}

ASSIGNED {
	v	(mV)		: membrane voltage
	i 	(nA)		: fluctuating current
	g_e	(umho)		: total excitatory conductance
	g_i	(umho)		: total inhibitory conductance
	g_e1	(umho)		: fluctuating excitatory conductance
	g_i1	(umho)		: fluctuating inhibitory conductance
	D_e	(umho umho /ms) : excitatory diffusion coefficient
	D_i	(umho umho /ms) : inhibitory diffusion coefficient
	exp_e
	exp_i
	amp_e	(umho)
	amp_i	(umho)
}

INITIAL {
	g_e1 = 0
	g_i1 = 0
	if(tau_e != 0) {
		D_e = 2 * std_e * std_e / tau_e
		exp_e = exp(-dt/tau_e)
		amp_e = std_e * sqrt( (1-exp(-2*dt/tau_e)) )
	}
	if(tau_i != 0) {
		D_i = 2 * std_i * std_i / tau_i
		exp_i = exp(-dt/tau_i)
		amp_i = std_i * sqrt( (1-exp(-2*dt/tau_i)) )
	}
}

BREAKPOINT {
	SOLVE oup
	if(tau_e==0) {
	   g_e = std_e * normrand(0,1)
	}
	if(tau_i==0) {
	   g_i = std_i * normrand(0,1)
	}
	g_e = g_e0 + g_e1
	g_i = g_i0 + g_i1
	i = g_e * (v - E_e) + g_i * (v - E_i)
}


PROCEDURE oup() {		: use Scop function normrand(mean, std_dev)
   if(tau_e!=0) {
	g_e1 =  exp_e * g_e1 + amp_e * normrand(0,1)
   }
   if(tau_i!=0) {
	g_i1 =  exp_i * g_i1 + amp_i * normrand(0,1)
   }
}


PROCEDURE new_seed(seed) {		: procedure to set the seed
	set_seed(seed)
	VERBATIM
	  printf("Setting random generator with seed = %g\n", _lseed);
	ENDVERBATIM
}

Now, I am using the Matlab interface, where I call the governing hoc file (sustainedfiring.hoc) for each run.

Code: Select all

 [status,result] = system('C:\nrn72\bin\nrniv.exe -nobanner sustainedfiring.hoc -c quit()');

So every iteration produces the exact same noise waveform. Correct me if I am wrong, but this is because the random number seed initializes to the exact same state, then produces the same waveform. If this is the case, is there a good way to create a random initialization?

Thanks all!

[ted]

That _is_ "the good way" to initialize a model that involves pseudorandom number generators. Reproducibility is essential for development and debugging. Failure to produce the same result on each run is a bug (usually a symptom of incorrect initialization).

"Fine, but assuming I have a program that seems to work properly, and produces the same result on each run, how can I introduce deliberate variation from run to run?"

Use a precalculated list of seeds. These could be embedded in your source code, or (better) saved to a text file. Example:

Code: Select all

 . . . all the code that sets up your model . . .
 . . . statements that read the seed file into a Vector . . .
for i=0, seedvec.size()-1 {
  . . . use seedvec.x[i] to perturb your random number generator . . .
  run()
  . . . do whatever needs to be done after a run . . .
}

[bscoventry]

Hi Ted,

Ok, that is what I was thinking I needed to do. Thank you for your input. I'll try this and let you know how it works out.

Best,
Brandon

[bscoventry]

Hi all,

Ok, I'm following the same idea listed here, albeit slightly differently to suit my needs. Unforntunately I can't quite understand why a particular procedure is not working. So, again here is the mod file code I'm using. In particular I'm trying to use the new_seed function:

Code: Select all

COMMENT
-----------------------------------------------------------------------------

	Fluctuating conductance model for synaptic bombardment
	======================================================
Gfluct2: conductance cannot be negative


REFERENCE

  Destexhe, A., Rudolph, M., Fellous, J-M. and Sejnowski, T.J.  
  Fluctuating synaptic conductances recreate in-vivo--like activity in
  neocortical neurons. Neuroscience 107: 13-24 (2001).

  (electronic copy available at http://cns.iaf.cnrs-gif.fr)


  A. Destexhe, 1999
-----------------------------------------------------------------------------
ENDCOMMENT



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

NEURON {
	POINT_PROCESS Gfluct2
	RANGE g_e, g_i, E_e, E_i, g_e0, g_i0, g_e1, g_i1
	RANGE std_e, std_i, tau_e, tau_i, D_e, D_i
	RANGE new_seed
	NONSPECIFIC_CURRENT i
}

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

PARAMETER {
	dt		(ms)

	E_e	= 0 	(mV)	: reversal potential of excitatory conductance
	E_i	= -75 	(mV)	: reversal potential of inhibitory conductance

	g_e0	= 0.0121 (umho)	: average excitatory conductance
	g_i0	= 0.0573 (umho)	: average inhibitory conductance

	std_e	= 0.0030 (umho)	: standard dev of excitatory conductance
	std_i	= 0.0066 (umho)	: standard dev of inhibitory conductance

	tau_e	= 2.728	(ms)	: time constant of excitatory conductance
	tau_i	= 10.49	(ms)	: time constant of inhibitory conductance
}

ASSIGNED {
	v	(mV)		: membrane voltage
	i 	(nA)		: fluctuating current
	g_e	(umho)		: total excitatory conductance
	g_i	(umho)		: total inhibitory conductance
	g_e1	(umho)		: fluctuating excitatory conductance
	g_i1	(umho)		: fluctuating inhibitory conductance
	D_e	(umho umho /ms) : excitatory diffusion coefficient
	D_i	(umho umho /ms) : inhibitory diffusion coefficient
	exp_e
	exp_i
	amp_e	(umho)
	amp_i	(umho)
}

INITIAL {
	g_e1 = 0
	g_i1 = 0
	if(tau_e != 0) {
		D_e = 2 * std_e * std_e / tau_e
		exp_e = exp(-dt/tau_e)
		amp_e = std_e * sqrt( (1-exp(-2*dt/tau_e)) )
	}
	if(tau_i != 0) {
		D_i = 2 * std_i * std_i / tau_i
		exp_i = exp(-dt/tau_i)
		amp_i = std_i * sqrt( (1-exp(-2*dt/tau_i)) )
	}
}

BREAKPOINT {
	SOLVE oup
	if(tau_e==0) {
	   g_e = std_e * normrand(0,1)
	}
	if(tau_i==0) {
	   g_i = std_i * normrand(0,1)
	}
	g_e = g_e0 + g_e1
	if(g_e < 0) { g_e = 0 }
	g_i = g_i0 + g_i1
	if(g_i < 0) { g_i = 0 }
	i = g_e * (v - E_e) + g_i * (v - E_i)
}


PROCEDURE oup() {		: use Scop function normrand(mean, std_dev)
   if(tau_e!=0) {
	g_e1 =  exp_e * g_e1 + amp_e * normrand(0,1)
   }
   if(tau_i!=0) {
	g_i1 =  exp_i * g_i1 + amp_i * normrand(0,1)
   }
}


PROCEDURE new_seed(seed) {		: procedure to set the seed
	set_seed(seed)
	VERBATIM
	  printf("Setting random generator with seed = %g\n", _lseed);
	ENDVERBATIM
}

Now, I'm trying to call new_seed in the hoc file, as follows:

Code: Select all

//soma setdata_Gfluct2(0.5)
new_seed_Gfluct2(.1)

Note that I'm not using a range variable in the seed, so the first line shouldn't be needed. I have however tried it just in case. Here are the two error messages I receive when I call the code, first with the first line commented then with it uncommented.

rniv: new_seed_Gfluct2 undefined function
in SustainedFiring.hoc near line 162
new_seed_Gfluct2(.1)
^
new_seed_Gfluct2(0.1 )
xopen("SustainedFiring.hoc" )
NEURONMainMenu[0].execute1("{xopen("SustainedFiring.hoc")}" )
NEURONMainMenu[0].load_file(0"./SustainedFiring.hoc" , )
and others

nrniv: setdata_Gfluct2 undefined function
in SustainedFiring.hoc near line 161
soma setdata_Gfluct2(0.5)
^
setdata_Gfluct2(0.5 )
xopen("SustainedFiring.hoc" )
NEURONMainMenu[0].execute1("{xopen("SustainedFiring.hoc")}" )
NEURONMainMenu[0].load_file(0"./SustainedFiring.hoc" , )
and others

I know this isn't quite on topic, but its completely baffling me.

Brandon

[ted]

This is one time when it's simpler than you anticipated. So simple, in fact, that it is all too easy to forget exactly what to do (see if that doesn't happen some time in the future!).

Individual point processes are instances of a class. hoc uses standard dot syntax for accessing the public members of a class. Example:

Code: Select all

objref gf
soma gf = new Gfluct2(0.5)
gf.new_seed(0.1)

[bscoventry]

Hi Ted!

You're right, that was pretty simple. I'm not exactly new to neuron, but there are always a few things I miss, haha. I just wanted to let you know that the code is now up and running, and I've also now able to keep a log of the variation. Thank you for all your help!

Brandon