Time step issue with normrand() in the mod file

[Skyfire79]

Hi, Ted,

I used normrand() in the mod file to generate a noisy current injection:

Code: Select all

NEURON {
  POINT_PROCESS IClampNoise
  RANGE i,del,dur,f0,f1,r,torn,std,bias
  ELECTRODE_CURRENT i
}

UNITS {
  (nA) = (nanoamp)
}

PARAMETER {
  del = 50    (ms)
  dur = 200   (ms)
  torn= 500   (ms)
  std = 0.2   (nA)
  f0  = 0.2   (nA)
  f1  = 0.8   (nA)
  r   = 60
  bias = 0    (nA)
  pi = 3.1415926
}

ASSIGNED {
  ival (nA)
  i (nA)
  amp (nA)
  noise (nA)
  on (1)
}

INITIAL {
  i  = 0
  on = 0
  net_send(del, 1)
}

BEFORE BREAKPOINT {
  if  (on) {
    noise = normrand(0,std*1(/nA))*1(nA)
    amp = f0 + 0.5*(f1-f0)*(tanh((t-torn)/(r/3)/(1(ms))-3)+1)
    ival = amp + noise + bias
  } else {
    ival = 0
  }
}

BREAKPOINT {
  i = ival
}

NET_RECEIVE (w) {
  if (flag == 1) {
    if (on == 0) {
      : turn it on
      on = 1
      : prepare to turn it off
      net_send(dur, 1)
    } else {
      : turn it off
      on = 0
    }
  }
}

In my hoc file, I also used the NetStim object to generate random spike trains to my network neurons. Now to my surprise, I discovered that the random spike trains change when I change the time step dt, even if I used the same seed for the NetStim. When I removed the normrand() function in the mod file, this won't happen. Since the random process changes with dt, it is hard for to me know whether the dt I am using is small enough for my simulation. Any way I can fix the random process from NetStim if I am using normrand() in the mod file?

Thank you!

Best,
Guoshi

[ted]

The "random process" does not change with dt. The problem you encountered is a consequence of the fact that a software "random number generator" is really a "pseudorandom sequence generator." Think of it as a function that returns a new value every time you call it. The sequence that it returns depends on its intial condition. The value that it returns on the nth call has nothing to do with dt--it depends only on the pseudorandom sequence generator's initial condition, and the fact that you have already called it n-1 times.

So how do you determine whether dt is short enough for accurate solutions? Easy. Turn off your "noise" sources before you run your test simulations. If dt is short enough with noise off, it will be short enough with noise on.

If you're planning to run a series of computational experiments that use "noise" generated by pseudorandom sequence generators, you better plan to use the same dt in all simulations.

[Skyfire79]

Thank you very much for your prompt response, Ted.

Yes, you are right that the "random process" does not change with dt. What I meant in the earlier post is that the spike trains from NetStim are not identical (the timing and number change) when I change dt if I have the normrand() function in the mod file. That suggests that somehow normrand() interferes with the spike generation from the Netstim, making it vary with dt. I am just wondering if there is any way to remove that interference.

Your approach to test dt sounds good. But is it possible that introducing noise would change the dynamics of the system which requires smaller dt to converge? In my simulation, sometime even if I used dt = 0.005 ms (Crank-Nicholson method), I still could not get identical spike output with dt = 0.01 ms (similar though). But the simulation will run very slow for such small dt. I have graded synapses in my network. Will that require a smaller dt to converge compared to the NetCon synapses? What is a reasonable dt for the CN integration method?

Thanks again!

Best,
GL

[ted]

What I meant in the earlier post is that the spike trains from NetStim are not identical (the timing and number change) when I change dt if I have the normrand() function in the mod file. That suggests that somehow normrand() interferes with the spike generation from the Netstim, making it vary with dt. I am just wondering if there is any way to remove that interference.

Ordinarily everything that uses a pseudorandom sequence generator draws from the same generator. Consequently if you have a NetStim with noise parameter > 0, plus an instance of your "noise current source" mechanism, changing dt makes the noise current source call the generator a different number of times per millisecond of model time. So the NetStim will get a different number from the generator. The solution is to give each mehanism instance its own pseudorandom generator instance. It's easy to do for the NetStim class, and is discussed elsewhere in the Forum. However, the NMODL code for your "noise current source" mechanism would have to be modified so that each instance of it could be paired with its own pseudorandom sequence generator.

Your approach to test dt sounds good. But is it possible that introducing noise would change the dynamics of the system which requires smaller dt to converge?

No. A current source does not affect the system equations or the eigenvalues of the system matrix.

In my simulation, sometime even if I used dt = 0.005 ms (Crank-Nicholson method), I still could not get identical spike output with dt = 0.01 ms (similar though).

Different with results with different dt is to be expected from a numerical solution.

I have graded synapses in my network. Will that require a smaller dt to converge compared to the NetCon synapses?

Not necessarily.

What is a reasonable dt for the CN integration method?

Whatever you decide it takes to get sufficient accuracy.

[Skyfire79]

Ted: thank you so much for your detailed answers to my questions. Your response makes perfect sense to me. I will try to find from the forum how to give each NetStim instance its own pseudorandom generator. However, I am not sure how to do that for the NMODL code for the "noise current source" mechanism. Could you give me more information or guidance on that?

Thanks!
GL

[hines]

You should modify your noise current source mod file to use the same code pattern for randomness that is implemented in the nrn/src/nrnoc/netstim.mod file.
The relevant fragment, which you can put at the end of your mod file, is

Code: Select all

NEURON { POINTER donotuse   THREADSAFE} : to hold a reference to a Hoc Random instance.
ASSIGNED{ donotuse }
VERBATIM
double nrn_random_pick(void* r);
void* nrn_random_arg(int argpos);
ENDVERBATIM
FUNCTION myrand() {
VERBATIM  
    if (_p_donotuse) {
        _lmyrand = nrn_random_pick(_p_donotuse);
    }else{
        /* only can be used in main thread */
        if (_nt != nrn_threads) {
            hoc_execerror("multithread random in NetStim"," only via hoc Random");
        }
ENDVERBATIM
        :your old random stuff goes here which is used if you are not useing hoc Random
        erand = ....
VERBATIM  
    }
ENDVERBATIM
}
PROCEDURE noiseFromRandom() {
VERBATIM
 {
        void** pv = (void**)(&_p_donotuse);
        if (ifarg(1)) {
                *pv = nrn_random_arg(1); 
        }else{
                *pv = (void*)0;
        }  
 }
ENDVERBATIM
}

Then in hoc, if your POINT_PROCESS is called Foo, you would attach an instance of Random to an instance of foo as in
objref foo, ran
foo = new Foo(.5)
ran = new Random()
foo.noise_from_random(ran)

and then set the proper distribution for ran ( ran.norm(...)?)
and now for statistically independent, reproducible random streams use
ran.MCellRan4(...) or ran.Random123(id1, id2).
http://www.neuron.yale.edu/neuron/stati ... #MCellRan4
http://www.neuron.yale.edu/neuron/stati ... #Random123
The latter is new in 7.3 but I highly recommend trying it out. id1 and id2 give you tremendous flexibility ind
defining your independent reproducible, streams. Another integer (ran.Random123_globalindex(id3)) means that if you run
n simulations with
for id3=0, n { anyraninstance.Random123_globalindex(id3) run() } then all streams of all runs will be statistically independent.
Finally there is an internal id4 (all the id0,1,2,3 range from 0 to 2^32-1) which is the sequence index of the stream.
You can see that this is administratively simpler than MCellRan4 where all aspects including sequence have to be mapped into
just two integers.

By the way. Whatever you end up with, be sure to test with printf statements so there are no surprises on your part about what is
going on.
I generally don't like picking a random number per time step unless you have carefully considered the effect of the dt dependence of the
noise frequency spectrum. The extreme case is that for constant parameters of the distribution, noisiness eventually disappears as dt -> 0.
ie. I faintly recall that noise parameters have to be related to the square root of dt. Perhapse someone could weigh in on this.
It would seem that one would want the noise spectrum to be independent of dt. I wonder if one should only pick at fixed Dt intervals
as one decreases dt.

[Skyfire79]

Hines: Thank you very much for your detailed information on the time step issue and I really appreciate it. I will try to see if I can work it out. You are right that maybe I should pick the random number at a fixed interval DT instead of dt. If the random number is picked at each dt, changing dt seems to change the noisy input by adding more random numbers (due to more time steps). But I am not sure how I can implement it at this time.

Thanks again for your help!

Best,
GL

[ted]

maybe I should pick the random number at a fixed interval DT instead of dt.
. . .
I am not sure how I can implement it at this time.

It could be done by taking advantage of events, which can be used to force execution of particular code at specific times. For example, events can be used to implement a signal source, such as a current or voltage clamp, that applies a waveform repeatedly or at certain times--see
I just want a current clamp that will deliver a sequence of current pulses at regular intervals.
and
I want a current clamp that will generate a pulse when I send it an event, or that I can use to produce pulses at precalculated times.
in the FAQ list http://www.neuron.yale.edu/neuron/faq/general-questions

However, the NMODL source code for your noise current source must have a NET_RECEIVE block, and at this time only point processes can have NET_RECEIVE blocks. So instead of inserting a density mechanism noise source into each section, your model setup code would have to attach a separate point process noise current source to every segment in each section that is supposed to have noise. That's easy to do--if every section is supposed to have noise in all segments, and the noise mechanism point process is called NIClamp (for "noisy current clamp"), this

Code: Select all

objref nsrc // will be a List that holds all instances of the noise mechanism
nsrc = new List()
forall for (x,0) nsrc.append(new NIClamp(x))

will do the trick.

Also, since a point process delivers current in nA, not in mA/cm2, the amplitude of the point process's output would have to be scaled by segment area. This is easily done inside the NMODL source code, since segment area in um2 is automatically available as an ASSIGNED with the name "area".

NIClamp would have a PARAMETER Dt that specifies the interval at which it picks a new value from its associated random number generator ("RNG"). At initialization, each instance of the NIClamp class would launch a self event that returns at t==Dt. This causes execution of code in the NET_RECEIVE block that does two things:
(1) picks a value from the associated RNG, which is used to calculate the new noise current that is injected into the segment,
and
(2) sends another self event that will return at t==t+Dt (so that a new value will be picked from the RNG)

[hines]

If your noise is used in the BREAKPOINT block an alternative to events is

BEFORE BREAKPOINT {
if (t >= tpick) {
r = pick the random number
tpick = t + Dt
}
}

and initialize tpick to 0.
Make tpick an ASSIGNED and Dt a PARAMETER
I didn't execute this so there may be typos but the idea is sound.

[ted]

If your noise is used in the BREAKPOINT block an alternative to events is

This would allow sticking with a density mechanism. Sounds good to me.

[Skyfire79]

Hines and Ted: Thanks a lot for taking time to address my question and they are all very helpful. I will try the approachs you suggested and hopefully can work it out.

Best,
GL