Recording variables from ARTIFICIAL_CELL

[rth]

I need to use the Wilson-Cowan Firing rate model as a counterpart to my multicompartment model. Borrowing some code from NetStim and ignoring that the e and i variables are inside the sigmoid function, I got something like this:

Code: Select all

NEURON { 
    THREADSAFE
    ARTIFICIAL_CELL WilsonCowanStim   
    RANGE internal_dt
    RANGE etau, itau
    RANGE eI12, eIscl
    RANGE iI12, iIscl
    RANGE wee, wei, wie, wii
    RANGE e, i
    RANGE e_inf, i_inf, e_in_tau, i_in_tau
    RANGE e_in, i_in
    :Generator
    BBCOREPOINTER donotuse
    RANGE ratescale
}

PARAMETER {
    internal_dt =   0.2
    etau        =   4.0
    itau        =  12.0
    eI12        =   1.3
    eIscl       =   4.
    iI12        =   2.
    iIscl       =   3.7
    wee         =  16. 
    wei         =  15.
    wie         = -12.
    wii         =  -3.
    e_inf       =   0.
    i_inf       =   0.
    e_in_tau    =   2.
    i_in_tau    =  10.
    ratescale   = 200.
}
ASSIGNED {
    e
    i
    eS0
    iS0
    e_in
    i_in
: Generator
	cspike
    event (ms)
	donotuse
}

INITIAL {
    e    = 0
    i    = 0
    e_in = e_inf
    i_in = i_inf
    eS0  = S(0,eI12,eIscl)
    iS0  = S(0,iI12,iIscl)
    
    event = invl(e*ratescale)
    cspike = 10
    net_send(internal_dt, 2)
    net_send(event      , cspike)
}

NET_RECEIVE (w) {
    LOCAL evt

	if (flag == 0) { : external event
        if (w > 0) {
            e_in = e_in + w
        } else {
            i_in = i_in - w
        }
    }
    if (flag == 2) { : compute next time stape
        update()
        net_send(internal_dt, 2)
        evt = t + invl(e*ratescale)
        if (evt < event){
            event = evt
            cspike = cspike + 1
            if (cspike > 10000){
                cspike = 10
            }
            net_send(event-t, cspike)
            :DB>>
            :printf("Send  t=%g, e=%g, x=%g @ %g\n",t,e, cspike,event)
            :<<DB
        }
    } else {
        :DB>>
        :printf("Event t=%g, f=%g:c=%g\n",t,flag,cspike)
        :<<DB
        :net_event(t)
        if (flag == cspike){
            net_event(t)
            event = t + invl(e*ratescale)
            cspike = cspike + 1
            if (cspike > 10000){
                cspike = 10
            }
            net_send(event-t, cspike)
            :DB>>
            :printf("Send  t=%g, e=%g, x=%g @ %g\n",t,e, cspike,event)
            :<<DB
        }
    }
}

PROCEDURE update(){
    LOCAL Fe, Fi, ze, zi
    Fe = Se(wee*e+wie*i+e_in)
    Fi = Si(wei*e+wii*i+i_in)
    ze = 1+Fe
    zi = 1+Fi
    e = Fe/ze + (e-Fe/ze)*exp(-ze*internal_dt/etau)
    i = Fi/zi + (i-Fi/zi)*exp(-zi*internal_dt/itau)
    e_in = e_inf + (e_in - e_inf)*exp(-internal_dt/e_in_tau)
    i_in = i_inf + (i_in - i_inf)*exp(-internal_dt/i_in_tau)
}

FUNCTION S(x,a,b){
    S=1./(1.+exp(-a*(x-b)))
}

FUNCTION Se(x){
    Se = S(x,eI12,eIscl) - eS0
}

FUNCTION Si(x){
    Si = S(x,iI12,iIscl) - iS0
}

VERBATIM
#if NRNBBCORE /* running in CoreNEURON */

#define IFNEWSTYLE(arg) arg

#else /* running in NEURON */

/*
   1 means noiseFromRandom was called when _ran_compat was previously 0 .
   2 means noiseFromRandom123 was called when _ran_compat was previously 0.
*/
static int _ran_compat; /* specifies the noise style for all instances */
#define IFNEWSTYLE(arg) if(_ran_compat == 2) { arg }

#endif /* running in NEURON */
ENDVERBATIM

:backward compatibility
PROCEDURE seed(x) {
VERBATIM
#if !NRNBBCORE
ENDVERBATIM
	set_seed(x)
VERBATIM
#endif
ENDVERBATIM
}


FUNCTION invl(rate) (ms) {
    LOCAL mean
    if (rate <= 0.){
        mean = 1e19
    } else {
        mean = 1000/rate
    }
	if (mean <= 0.) {
		mean = .1 (ms)
	}
    invl = mean*erand()
}

FUNCTION erand() {
VERBATIM
	if (_p_donotuse) {
		/*
		:Supports separate independent but reproducible streams for
		: each instance. However, the corresponding hoc Random
		: distribution MUST be set to Random.negexp(1)
		*/
#if !NRNBBCORE
		if (_ran_compat == 2) {
			_lerand = nrnran123_negexp((nrnran123_State*)_p_donotuse);
		}else{
			_lerand = nrn_random_pick(_p_donotuse);
		}
#else
		_lerand = nrnran123_negexp((nrnran123_State*)_p_donotuse);
#endif
		return _lerand;
	}else{
#if NRNBBCORE
		assert(0);
#else
		/*
		: the old standby. Cannot use if reproducible parallel sim
		: independent of nhost or which host this instance is on
		: is desired, since each instance on this cpu draws from
		: the same stream
		*/
#endif
	}
#if !NRNBBCORE
ENDVERBATIM
	erand = exprand(1)
VERBATIM
#endif
ENDVERBATIM
}

PROCEDURE noiseFromRandom() {
VERBATIM
#if !NRNBBCORE
 {
	void** pv = (void**)(&_p_donotuse);
	if (_ran_compat == 2) {
		fprintf(stderr, "NetStim.noiseFromRandom123 was previously called\n");
		assert(0);
	}
	_ran_compat = 1;
	if (ifarg(1)) {
		*pv = nrn_random_arg(1);
	}else{
		*pv = (void*)0;
	}
 }
#endif
ENDVERBATIM
}


PROCEDURE noiseFromRandom123() {
VERBATIM
#if !NRNBBCORE
 {
	nrnran123_State** pv = (nrnran123_State**)(&_p_donotuse);
	if (_ran_compat == 1) {
		fprintf(stderr, "NetStim.noiseFromRandom was previously called\n");
		assert(0);
	}
	_ran_compat = 2;
	if (*pv) {
		nrnran123_deletestream(*pv);
		*pv = (nrnran123_State*)0;
	}
	if (ifarg(3)) {
		*pv = nrnran123_newstream3((uint32_t)*getarg(1), (uint32_t)*getarg(2), (uint32_t)*getarg(3));
	}else if (ifarg(2)) {
		*pv = nrnran123_newstream((uint32_t)*getarg(1), (uint32_t)*getarg(2));
	}
 }
#endif
ENDVERBATIM
}

DESTRUCTOR {
VERBATIM
	if (_p_donotuse) {
#if NRNBBCORE
		{ /* but note that mod2c does not translate DESTRUCTOR */
#else
		if (_ran_compat == 2) {
#endif
			nrnran123_State** pv = (nrnran123_State**)(&_p_donotuse);
			nrnran123_deletestream(*pv);
			*pv = (nrnran123_State*)0;
		}
	}
ENDVERBATIM
}

VERBATIM
static void bbcore_write(double* x, int* d, int* xx, int *offset, _threadargsproto_) {
	/* error if using the legacy scop_exprand */
	if (!_p_donotuse) {
		fprintf(stderr, "NetStim: cannot use the legacy scop_negexp generator for the random stream.\n");
		assert(0);
	}
	if (d) {
		char which;
		uint32_t* di = ((uint32_t*)d) + *offset;
#if !NRNBBCORE
		if (_ran_compat == 1) {
			void** pv = (void**)(&_p_donotuse);
			/* error if not using Random123 generator */
			if (!nrn_random_isran123(*pv, di, di+1, di+2)) {
				fprintf(stderr, "NetStim: Random123 generator is required\n");
				assert(0);
			}
			nrn_random123_getseq(*pv, di+3, &which);
			di[4] = (int)which;
		}else{
#else
    {
#endif
			nrnran123_State** pv = (nrnran123_State**)(&_p_donotuse);
			nrnran123_getids3(*pv, di, di+1, di+2);
			nrnran123_getseq(*pv, di+3, &which);
			di[4] = (int)which;
#if NRNBBCORE
			/* CORENeuron does not call DESTRUCTOR so... */
			nrnran123_deletestream(*pv);
                        *pv = (nrnran123_State*)0;
#endif
		}
		/*printf("Netstim bbcore_write %d %d %d\n", di[0], di[1], di[3]);*/
	}
	*offset += 5;
}

static void bbcore_read(double* x, int* d, int* xx, int* offset, _threadargsproto_) {
	/* Generally, CoreNEURON, in the context of psolve, begins with
           an empty model so this call takes place in the context of a freshly
           created instance and _p_donotuse is not NULL.
	   However, this function
           is also now called from NEURON at the end of coreneuron psolve
           in order to transfer back the nrnran123 sequence state. That
           allows continuation with a subsequent psolve within NEURON or
           properly transfer back to CoreNEURON if we continue the psolve
           there. So now, extra logic is needed for this call to work in
           a NEURON context.
        */

	uint32_t* di = ((uint32_t*)d) + *offset;
#if NRNBBCORE
	nrnran123_State** pv = (nrnran123_State**)(&_p_donotuse);
	assert(!_p_donotuse);
	*pv = nrnran123_newstream3(di[0], di[1], di[2]);
	nrnran123_setseq(*pv, di[3], (char)di[4]);
#else
	uint32_t id1, id2, id3;
	assert(_p_donotuse);
	if (_ran_compat == 1) { /* Hoc Random.Random123 */
		void** pv = (void**)(&_p_donotuse);
		int b = nrn_random_isran123(*pv, &id1, &id2, &id3);
		assert(b);
		nrn_random123_setseq(*pv, di[3], (char)di[4]);
	}else{
		assert(_ran_compat == 2);
		nrnran123_State** pv = (nrnran123_State**)(&_p_donotuse);
		nrnran123_getids3(*pv, &id1, &id2, &id3);
		nrnran123_setseq(*pv, di[3], (char)di[4]);
	}
        /* Random123 on NEURON side has same ids as on CoreNEURON side */
	assert(di[0] == id1 && di[1] == id2 && di[2] == id3);
#endif
	*offset += 5;
}
ENDVERBATIM

Now I need to record population firing rates, i.e. e and i variables.

Code: Select all

from numpy import *
from neuron import h
from matplotlib.pyplot import *



wcs = h.WilsonCowanStim()
wcs.ratescale = 200
wcs.e_inf     = 1.5
rds = h.Random()
rds.negexp(1)             # set random # generator using negexp(1) - avg interval in NetStim
sead = random.randint(0,32562)
rds.MCellRan4(sead,sead)  # seeds are in order, shouldn't matter                       
wcs.noiseFromRandom(rds)  # use random # generator for this NetStim 


trec = h.Vector()
trec.record(h._ref_t)
erec = h.Vector()
erec.record(wcs._ref_e)
irec = h.Vector()
irec.record(wcs._ref_i)

srec = h.Vector()
ncout = h.NetCon(wcs, None)
ncout.record(srec)

hpc  = h.ParallelContext()
hpc.nthread(2)

h.finitialize()
h.fcurrent()
h.frecord_init()
h.dt = 0.1

while h.t < 1000 :
    h.fadvance()

plot(array(trec),array(erec),'-')
plot(array(trec),array(irec),'-')
srec = array(srec)
plot(srec,ones(srec.shape[0])*0.5,'k|')

show()

This code returns an error
NEURON: Section access unspecified
near line 0
objref hoc_obj_[2]
^
Vector[0].record(...)
Traceback (most recent call last):
File "/media/rth/rth-core/simulations/exampls/neuron-controledP/WilsonCowanStim_test.py", line 18, in <module>
trec.record(h._ref_t)
RuntimeError: hocobj_call error


The error message changes if I create a "fake" section and associate everything with it.

Code: Select all

from numpy import *
from neuron import h
from matplotlib.pyplot import *


sec = h.Section()
wcs = h.WilsonCowanStim(sec=sec)
wcs.ratescale = 200
wcs.e_inf     = 1.5
rds = h.Random(sec=sec)
rds.negexp(1)             # set random # generator using negexp(1) - avg interval in NetStim
sead = random.randint(0,32562)
rds.MCellRan4(sead,sead)  # seeds are in order, shouldn't matter                       
wcs.noiseFromRandom(rds)  # use random # generator for this NetStim 


trec = h.Vector()
trec.record(h._ref_t,sec=sec)
erec = h.Vector()
erec.record(wcs._ref_e,sec=sec)
irec = h.Vector()
irec.record(wcs._ref_i,sec=sec)

srec = h.Vector()
ncout = h.NetCon(wcs, None)
ncout.record(srec)

hpc  = h.ParallelContext()
hpc.nthread(2)

h.finitialize()
h.fcurrent()
h.frecord_init()
h.dt = 0.1

while h.t < 1000 :
    h.fadvance()

plot(array(trec),array(erec),'-')
plot(array(trec),array(irec),'-')
srec = array(srec)
plot(srec,ones(srec.shape[0])*0.5,'k|')

show()

Which returns
NEURON: We were unable to associate a PlayRecord item with a thread
near line 0
objref hoc_obj_[2]
^
finitialize()
Traceback (most recent call last):
File "/media/rth/rth-core/simulations/exampls/neuron-controledP/WilsonCowanStim_test.py", line 31, in <module>
h.finitialize()
RuntimeError: hocobj_call error
The error above does not appear if ParallelContext and multithreading are commented out.

Question:
How to record variables in ARTIFICIAL_CELL, which is technically not associated with any section of the compartment-based model?

[ted]

There are two problems. First is accessing the values of the variables that you want. Second is accessing them at the times that you want. Typically you want to capture these variables at fixed intervals, not just when the ARTIFICIAL_CELL receives an event.

The easiest case is when the variables are analytic functions of time, e.g. IntFire1's "membrane activation variable" m. All you have to do is add those functions to your mod file, and force them to be evaluated at every time step. For IntFire1, the function would be

Code: Select all

FUNCTION M() {
  M = m*exp(-(t - t0)/tau)
}

which takes advantage of the fact that IntFire1 updates m only when it receives an event--which means that m always has whatever value it had immediately after the most recent event arrived.

And I would have to give my modified IntFire1 its own class name to avoid a conflict with the built-in IntFire1. Maybe IntFire1x.

Great. Now how to call IntFire1x.M after each fadvance(), and how to get the value it returns into a Vector? You could use an FinitializeHandler and CVode.event() to force evaluation of M at regular intervals, or you could write your own procedure that forces M to be evaluated at regular intervals. In hoc or Python I would do something like this bit of pseudocode (which assumes that ifx is the objref or Python variable that references the IntFirex instance) :

Code: Select all

create a Vector called mvec

define a procedure called xrun which does this:
  resize mvec to 0 elements
  finitialize() (or h.finitialize() if you're using Python)
  while (t<tstop) {
    append ifx.M() to mvec
    fadvance() (or h.fadvance() if you're using Python)
  }
  append ifx.M() to mvec // to get that last value of m

Then you simply call xrun() and wait for the results.

[rth]

Thank you, Ted, for your thoughts.

The overall bulky construction is because the e-variable of the WC model controls the firing rate of the Poisson generator. Therefore, the rate needs to be updated, and a new time to the next spike needs to be recalculated. If the new time is 'sooner' than the previous one, a new event is issued, and the old event is ignored. So, there are many net_send calls with different flags.

WC is difficult to convert into anything analytical, but by ignoring the e and i variables in the sigmoid function, equations can be solved by the exp Euler method, which I call every 0.2 ms. The procedure update in the mod file actually does that. Net events with flag=2 update dynamic variables, recompute new event time evt, and send new "time to spike" event with flag > 10. So if a received event has a flag that is not equal to the current event counter, that is an old event, and we should ignore it. Finally, if the flag is zero, it is an external input spike, and we should update the input firing rate.

Unfortunately, performance here is a key element, the simulations run from a few tens of minutes to hours in model time (we model a very peculiar developmental processes https://elifesciences.org/articles/84333), so it is not a good idea to put anything in the main loop. Vector.record mechanisms could be the perfect choice for recording variables.

The NEURON documentation says

For the local variable timestep method, CVode.use_local_dt() and/or multiple threads, ParallelContext.nthread() , it is often helpful to provide specific information about which cell the var pointer is associated with by inserting as the first arg some POINT_PROCESS object which is located on the cell. This is necessary if the pointer is not a RANGE variable and is much more efficient if it is. The fixed step and global variable time step method do not need or use this information for the local step method but will use it for multiple threads. It is therefore a good idea to supply it if possible.

However, what to do with ARTIFICIAL_CELL is not clear. I tried to use POINT_PROCESS approach for ARTIFICIAL_CELL and insert it into a compartment, but it returns an error.

Code: Select all

from numpy import *
from neuron import h
from matplotlib.pyplot import *

wcs_sec = h.Section()
wcs = h.WilsonCowanStim(0.5, sec=wcs_sec)
wcs.ratescale = 200
wcs.e_inf     = 1.5
rds = h.Random(0.5, sec=wcs_sec)
rds.negexp(1)             # set random # generator using negexp(1) - avg interval in NetStim
sead = random.randint(0,32562)
rds.MCellRan4(sead,sead)  # seeds are in order, shouldn't matter                       
wcs.noiseFromRandom(rds)  # use random # generator for this NetStim 


trec = h.Vector(0.5, sec=wcs_sec)
trec.record(h._ref_t)
erec = h.Vector(0.5, sec=wcs_sec)
erec.record(wcs._ref_e)
irec = h.Vector(0.5, sec=wcs_sec)
irec.record(wcs._ref_i)

srec = h.Vector(0.5, sec=wcs_sec)
ncout = h.NetCon(wcs, None)
ncout.record(srec)

hpc  = h.ParallelContext()
hpc.nthread(8)

h.finitialize()
h.fcurrent()
h.frecord_init()
h.dt = 0.1

while h.t < 1000 :
    h.fadvance()

plot(array(trec),array(erec),'-')
plot(array(trec),array(irec),'-')
srec = array(srec)
plot(srec,ones(srec.shape[0])*0.5,'k|')

show()

NEURON: We were unable to associate a PlayRecord item with a thread
near line 0
objref hoc_obj_[2]
^
finitialize()
Traceback (most recent call last):
File "/media/rth/rth-core/simulations/exampls/neuron-controledP/WilsonCowanStim_test.py", line 30, in <module>
h.finitialize()
RuntimeError: hocobj_call error


The only option I have at this point is to grab code from the Vector class and staff my mod file with it, but I still hope there is a workaround.

[rth]

It behaves very strange, for example the code below works

Code: Select all

from numpy import *
from neuron import h
from matplotlib.pyplot import *

wcs = h.WilsonCowanStim()
wcs.ratescale = 200
wcs.e_inf     = 1.5
rds = h.Random()
rds.negexp(1)             # set random # generator using negexp(1) - avg interval in NetStim
sead = random.randint(0,32562)
rds.MCellRan4(sead,sead)  # seeds are in order, shouldn't matter                       
wcs.noiseFromRandom(rds)  # use random # generator for this NetStim 


trec = h.Vector()
trec.record(h._ref_t, 1)
erec = h.Vector()
erec.record(wcs._ref_e,1)
irec = h.Vector()
irec.record(wcs._ref_i,1)

srec = h.Vector()
ncout = h.NetCon(wcs, None)
ncout.record(srec)

hpc  = h.ParallelContext()
hpc.nthread(8)

h.finitialize()
h.fcurrent()
h.frecord_init()
h.dt = 0.1

while h.t < 1000 :
    h.fadvance()

plot(array(trec),array(erec),'-')
plot(array(trec),array(irec),'-')
srec = array(srec)
plot(srec,ones(srec.shape[0])*0.5,'k|')

show()

But as soon as I add any section, it returns the same error

Code: Select all

from numpy import *
from neuron import h
from matplotlib.pyplot import *

# I just added this section!
#vvvvvvvvvvvvvvvvvvv
wcs_sec = h.Section()
#^^^^^^^^^^^^^^^^^^^^

wcs = h.WilsonCowanStim()
wcs.ratescale = 200
wcs.e_inf     = 1.5
rds = h.Random()
rds.negexp(1)             # set random # generator using negexp(1) - avg interval in NetStim
sead = random.randint(0,32562)
rds.MCellRan4(sead,sead)  # seeds are in order, shouldn't matter                       
wcs.noiseFromRandom(rds)  # use random # generator for this NetStim 


trec = h.Vector()
trec.record(h._ref_t, 1)
erec = h.Vector()
erec.record(wcs._ref_e,1)
irec = h.Vector()
irec.record(wcs._ref_i,1)

srec = h.Vector()
ncout = h.NetCon(wcs, None)
ncout.record(srec)

hpc  = h.ParallelContext()
hpc.nthread(8)

h.finitialize()
h.fcurrent()
h.frecord_init()
h.dt = 0.1

while h.t < 1000 :
    h.fadvance()

plot(array(trec),array(erec),'-')
plot(array(trec),array(irec),'-')
srec = array(srec)
plot(srec,ones(srec.shape[0])*0.5,'k|')

show()

NEURON: We were unable to associate a PlayRecord item with a thread
near line 0
objref hoc_obj_[2]
^
finitialize()
Traceback (most recent call last):
File "/media/rth/rth-core/simulations/exampls/neuron-controledP/WilsonCowanStim_test.py", line 30, in <module>
h.finitialize()
RuntimeError: hocobj_call error

[hines]

First, it seems to me to be a mistake in the internal NEURON implementation that the following generates an error

Code: Select all

trec.record(wcs, h._ref_t)
erec.record(wcs, wcs._ref_e)

That is, I believe the thread that wcs (an ARTIFICIAL CELL) is in is definitely known just as any POINT_PROCESS. An extra section should not be needed.
I'll look into this in more detail, and if I'm correct, will make the needed code changes. This will happen in the master branch of the repository. What version of NEURON are you using? Indeed, I see that the if statement that generates the error is logically incorrect in its intent to allow ARTIFICIAL_CELL.

Second, I think you can reduce the use of the event queue by making use of net_move to move the existing SelfEvent to an earlier time. I.e. no need to use your cspike count strategy for the flag. Instead , in the INITIAL block use

Code: Select all

net_send(event, 1) : SelfEvents that can be moved must have flag 1.

and in the NET_RECEIVE block get rid of all mention of cspike and use

Code: Select all

    if (flag == 2) { : compute next time step
        update()
        net_send(internal_dt, 2)
        evt = t + invl(e*ratescale)
        if (evt < event){
            event = evt
            net_move(event
        }
    } else if (flag == 1) {
            net_event(t)   
            event = t + invl(e*ratescale)
            net_send(event-t, 1)
    }

Third, I suggest you use the Random123 generator instead of MCellRan4 as the former is easier to use for parallel simulations as each stream is statistically independent and reproducible.

[rth]

What version of NEURON are you using

The most fresh from pip.

Code: Select all

 python
Python 3.10.12 (main, Nov 20 2023, 15:14:05) [GCC 11.4.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> from neuron import h
>>> h.nrnversion()
'NEURON -- VERSION 8.2.3 HEAD (f0ed37010) 2023-09-14'
>>> 

The second and third are perfect suggestions! Thank you very much!

[rth]

This will happen in the master branch of the repository

Should I compile locally? Will it be pushed to pip soon?

[hines]

I've started a pull request https://github.com/neuronsimulator/nrn/pull/2858 . For now, you can work around the
"Optional first arg is not a POINT_PROCESS" error by using
hpc.nthread(1)
or else creating a section for each thread and putting a POINT_PROCESS in the middle of each Section. One would use the latter as the first arg for the record statements. Unfortunately, ARTIFICIAL_CELLs are distributed in round-robin fashion among threads and it is not easy for the user to determine which thread the ARTIFICIAL_CELL is in. So it is not clear which thread time the ARTIFICIAL_CELL should be associated with when plotting. I suppose
it does not matter if you are using the fixed step method.

I see that threads do matter for performance, I created 10000 more instances of your WilsonCowanStim() and ran three times each for nthread=1,2,4,8
The time (seconds) for a 1000ms run is

Code: Select all

nthread  1   time  5.21 6.31 5.12
nthread  2   time  2.89 3.33 3.12 
nthread  4   time  1.79 1.67 1.96
nthread  8   time  1.04 1.12 1.08

[hines]

I applied the pull request change above to the 8.2 branch (for a future 8.2.5 release) and observe better performance than the master (probably due to ARTIFICIAL_CELLS having better performance when the memory layout is "Array of Structures" rather than the newer "Structure of Array" we are experimenting with for a future 9.0 release). That is

Code: Select all

nthread  1   time  4.43 4.78 4.63
nthread  2   time  2.59 2.55 2.73
nthread  4   time  1.54 1.51 1.64
nthread  8   time  0.78 0.83 0.79

I also copied wcfr.mod (ARTIFICIAL_CELL WilsonCowanStim) to wcfrmove.mod (ARTIFICIAL_CELL WilsonCowanStimMove)
and modified to use net_move in place of one of the net_send statements to get times of

Code: Select all

nthread  1   time  3.60 3.50 3.73
nthread  2   time  1.95 1.96 2.08
nthread  4   time  1.37 1.25 1.37
nthread  8   time  0.65 0.64 0.67

[rth]

I have a bit of trouble with net_move, specifically when firing rates are high. It can miss events and then stops complitely. I would probably return to cspike counter. It generates lots of events, but somewhat more robust.

I'll try to make an example to show the problem.

[hines]

I'll try to make an example to show the problem.

I'll be interested to see the example.

[rth]

Here's how it looks. I use spikes from my multicompartment neurons to activate the WCS firing rate model.
The plots below are: (1) on the very top are input spikes, (2) spikes under it are WCS spikes, (3) blue and orange curves are e and i firing rates of WCS model, and (4) the wcs.e_in variable is plotted at the bottom. Again, I can record continuous variables from WCS only if I set it for recording, and there are no sections in the model (see above).

If WCS uses net_move (left plot), it fails to generate spikes, while if I do my spike counter (on the right), WCS works.


I'm trying to figure out when something like that is happening and to cook a simple example.

[rth]

Here's an example

My mod file

Code: Select all

NEURON { 
    THREADSAFE
    ARTIFICIAL_CELL WilsonCowanStim_cmove
    RANGE internal_dt
    RANGE etau, itau
    RANGE eI12, eIscl
    RANGE iI12, iIscl
    RANGE wee, wei, wie, wii
    RANGE e, i
    RANGE e_inf, i_inf, e_in_tau, i_in_tau
    RANGE e_in, i_in
    :Generator
    BBCOREPOINTER donotuse
    RANGE ratescale
}

PARAMETER {
    internal_dt =   0.2
    etau        =   4.0
    itau        =  12.0
    eI12        =   1.3
    eIscl       =   4.
    iI12        =   2.
    iIscl       =   3.7
    wee         =  16. 
    wei         =  15.
    wie         = -12.
    wii         =  -3.
    e_inf       =   0.
    i_inf       =   0.
    e_in_tau    =   2.
    i_in_tau    =  10.
    ratescale   = 200.
}
ASSIGNED {
    e
    i
    eS0
    iS0
    e_in
    i_in
: Generator
    event (ms)
    donotuse
}

INITIAL {
    e    = 0
    i    = 0
    e_in = e_inf
    i_in = i_inf
    eS0  = S(0,eI12,eIscl)
    iS0  = S(0,iI12,iIscl)
    
    event = invl(e*ratescale)
    net_send(internal_dt, 2)
    net_send(event      , 1)
}

NET_RECEIVE (w) {
    LOCAL evt

    if (flag == 0) { : external event
        if (w > 0) {
            e_in = e_in + w
        } else {
            i_in = i_in - w
        }
    }
    if (flag == 2) { : compute next time stape
        update()
        net_send(internal_dt, 2)
        evt = t + invl(e*ratescale)
        if (evt < event){
            event = evt
	    net_move(event)
        }
    }
    if (flag == 1) {
        net_event(t)
	event = t + invl(e*ratescale)
	net_send(event, 1)        
    }
}

PROCEDURE update(){
    LOCAL Fe, Fi, ze, zi
    Fe = Se(wee*e+wie*i+e_in)
    Fi = Si(wei*e+wii*i+i_in)
    ze = 1+Fe
    zi = 1+Fi
    e = Fe/ze + (e-Fe/ze)*exp(-ze*internal_dt/etau)
    i = Fi/zi + (i-Fi/zi)*exp(-zi*internal_dt/itau)
    e_in = e_inf + (e_in - e_inf)*exp(-internal_dt/e_in_tau)
    i_in = i_inf + (i_in - i_inf)*exp(-internal_dt/i_in_tau)
}

FUNCTION S(x,a,b){
    S=1./(1.+exp(-a*(x-b)))
}

FUNCTION Se(x){
    Se = S(x,eI12,eIscl) - eS0
}

FUNCTION Si(x){
    Si = S(x,iI12,iIscl) - iS0
}

VERBATIM
#if NRNBBCORE /* running in CoreNEURON */

#define IFNEWSTYLE(arg) arg

#else /* running in NEURON */

/*
   1 means noiseFromRandom was called when _ran_compat was previously 0 .
   2 means noiseFromRandom123 was called when _ran_compat was previously 0.
*/
static int _ran_compat; /* specifies the noise style for all instances */
#define IFNEWSTYLE(arg) if(_ran_compat == 2) { arg }

#endif /* running in NEURON */
ENDVERBATIM

:backward compatibility
PROCEDURE seed(x) {
VERBATIM
#if !NRNBBCORE
ENDVERBATIM
	set_seed(x)
VERBATIM
#endif
ENDVERBATIM
}


FUNCTION invl(rate) (ms) {
    LOCAL mean
    if (rate <= 0.){
        mean = 1e19
    } else {
        mean = 1000/rate
    }
	if (mean <= 0.) {
		mean = .1 (ms)
	}
    invl = mean*erand()
}

FUNCTION erand() {
VERBATIM
	if (_p_donotuse) {
		/*
		:Supports separate independent but reproducible streams for
		: each instance. However, the corresponding hoc Random
		: distribution MUST be set to Random.negexp(1)
		*/
#if !NRNBBCORE
		if (_ran_compat == 2) {
			_lerand = nrnran123_negexp((nrnran123_State*)_p_donotuse);
		}else{
			_lerand = nrn_random_pick(_p_donotuse);
		}
#else
		_lerand = nrnran123_negexp((nrnran123_State*)_p_donotuse);
#endif
		return _lerand;
	}else{
#if NRNBBCORE
		assert(0);
#else
		/*
		: the old standby. Cannot use if reproducible parallel sim
		: independent of nhost or which host this instance is on
		: is desired, since each instance on this cpu draws from
		: the same stream
		*/
#endif
	}
#if !NRNBBCORE
ENDVERBATIM
	erand = exprand(1)
VERBATIM
#endif
ENDVERBATIM
}

PROCEDURE noiseFromRandom() {
VERBATIM
#if !NRNBBCORE
 {
	void** pv = (void**)(&_p_donotuse);
	if (_ran_compat == 2) {
		fprintf(stderr, "NetStim.noiseFromRandom123 was previously called\n");
		assert(0);
	}
	_ran_compat = 1;
	if (ifarg(1)) {
		*pv = nrn_random_arg(1);
	}else{
		*pv = (void*)0;
	}
 }
#endif
ENDVERBATIM
}


PROCEDURE noiseFromRandom123() {
VERBATIM
#if !NRNBBCORE
 {
	nrnran123_State** pv = (nrnran123_State**)(&_p_donotuse);
	if (_ran_compat == 1) {
		fprintf(stderr, "NetStim.noiseFromRandom was previously called\n");
		assert(0);
	}
	_ran_compat = 2;
	if (*pv) {
		nrnran123_deletestream(*pv);
		*pv = (nrnran123_State*)0;
	}
	if (ifarg(3)) {
		*pv = nrnran123_newstream3((uint32_t)*getarg(1), (uint32_t)*getarg(2), (uint32_t)*getarg(3));
	}else if (ifarg(2)) {
		*pv = nrnran123_newstream((uint32_t)*getarg(1), (uint32_t)*getarg(2));
	}
 }
#endif
ENDVERBATIM
}

DESTRUCTOR {
VERBATIM
	if (_p_donotuse) {
#if NRNBBCORE
		{ /* but note that mod2c does not translate DESTRUCTOR */
#else
		if (_ran_compat == 2) {
#endif
			nrnran123_State** pv = (nrnran123_State**)(&_p_donotuse);
			nrnran123_deletestream(*pv);
			*pv = (nrnran123_State*)0;
		}
	}
ENDVERBATIM
}

VERBATIM
static void bbcore_write(double* x, int* d, int* xx, int *offset, _threadargsproto_) {
	/* error if using the legacy scop_exprand */
	if (!_p_donotuse) {
		fprintf(stderr, "NetStim: cannot use the legacy scop_negexp generator for the random stream.\n");
		assert(0);
	}
	if (d) {
		char which;
		uint32_t* di = ((uint32_t*)d) + *offset;
#if !NRNBBCORE
		if (_ran_compat == 1) {
			void** pv = (void**)(&_p_donotuse);
			/* error if not using Random123 generator */
			if (!nrn_random_isran123(*pv, di, di+1, di+2)) {
				fprintf(stderr, "NetStim: Random123 generator is required\n");
				assert(0);
			}
			nrn_random123_getseq(*pv, di+3, &which);
			di[4] = (int)which;
		}else{
#else
    {
#endif
			nrnran123_State** pv = (nrnran123_State**)(&_p_donotuse);
			nrnran123_getids3(*pv, di, di+1, di+2);
			nrnran123_getseq(*pv, di+3, &which);
			di[4] = (int)which;
#if NRNBBCORE
			/* CORENeuron does not call DESTRUCTOR so... */
			nrnran123_deletestream(*pv);
                        *pv = (nrnran123_State*)0;
#endif
		}
		/*printf("Netstim bbcore_write %d %d %d\n", di[0], di[1], di[3]);*/
	}
	*offset += 5;
}

static void bbcore_read(double* x, int* d, int* xx, int* offset, _threadargsproto_) {
	/* Generally, CoreNEURON, in the context of psolve, begins with
           an empty model so this call takes place in the context of a freshly
           created instance and _p_donotuse is not NULL.
	   However, this function
           is also now called from NEURON at the end of coreneuron psolve
           in order to transfer back the nrnran123 sequence state. That
           allows continuation with a subsequent psolve within NEURON or
           properly transfer back to CoreNEURON if we continue the psolve
           there. So now, extra logic is needed for this call to work in
           a NEURON context.
        */

	uint32_t* di = ((uint32_t*)d) + *offset;
#if NRNBBCORE
	nrnran123_State** pv = (nrnran123_State**)(&_p_donotuse);
	assert(!_p_donotuse);
	*pv = nrnran123_newstream3(di[0], di[1], di[2]);
	nrnran123_setseq(*pv, di[3], (char)di[4]);
#else
	uint32_t id1, id2, id3;
	assert(_p_donotuse);
	if (_ran_compat == 1) { /* Hoc Random.Random123 */
		void** pv = (void**)(&_p_donotuse);
		int b = nrn_random_isran123(*pv, &id1, &id2, &id3);
		assert(b);
		nrn_random123_setseq(*pv, di[3], (char)di[4]);
	}else{
		assert(_ran_compat == 2);
		nrnran123_State** pv = (nrnran123_State**)(&_p_donotuse);
		nrnran123_getids3(*pv, &id1, &id2, &id3);
		nrnran123_setseq(*pv, di[3], (char)di[4]);
	}
        /* Random123 on NEURON side has same ids as on CoreNEURON side */
	assert(di[0] == id1 && di[1] == id2 && di[2] == id3);
#endif
	*offset += 5;
}
ENDVERBATIM

The test example

Code: Select all

from numpy import *
from neuron import h
from matplotlib.pyplot import *

# wcs_sec = h.Section()
wcs = h.WilsonCowanStim_cmove()
wcs.ratescale = 5000
wcs.e_inf     = 0.5
rds = h.Random()
rds.negexp(1)             # set random # generator using negexp(1) - avg interval in NetStim
sead = random.randint(0,32562)
rds.Random123(sead,sead)  # seeds are in order, shouldn't matter                       
wcs.noiseFromRandom(rds)  # use random # generator for this NetStim 

inspk = h.NetStim()
inspk.start    = 200
inspk.interval = 0.2
inspk.number   = 10000

incon = h.NetCon(inspk,wcs)
incon.weight[0] = 0.1
isrec = h.Vector()
incon.record(isrec)

trec = h.Vector()
trec.record(h._ref_t,1)
erec = h.Vector()
erec.record(wcs._ref_e,1)
irec = h.Vector()
irec.record(wcs._ref_i,1)

srec = h.Vector()
ncout = h.NetCon(wcs, None)
ncout.record(srec)



hpc  = h.ParallelContext()
hpc.nthread(8)

h.finitialize()
h.fcurrent()
h.frecord_init()
h.dt = 0.1

while h.t < 1000 :
    h.fadvance()

plot(array(trec),array(erec),'-')
plot(array(trec),array(irec),'-')
srec = array(srec)
isrec = array(isrec)
plot(srec ,ones(srec.shape[0] )*0.5,'k|')
plot(isrec,ones(isrec.shape[0])*0.6,'k|')

show()

And result.

[rth]

I found a bug in my code.
Here's my NET_RECEIVE function, which works without problems

Code: Select all

NET_RECEIVE (w) {
    LOCAL evt

    if (flag == 0) { : external event
        if (w > 0) {
            e_in = e_in + w
        } else {
            i_in = i_in - w
        }
    }
    if (flag == 2) { : compute next time step
        update()
        net_send(internal_dt, 2)
        evt = t + invl(e*ratescale)
        if (evt < event){
            event = evt
	    net_move(event)
        }
    }
    if (flag == 1) { : my spike event
        net_event(t)
	evt = invl(e*ratescale)
	event = t + evt
	net_send(evt, 1)        
    }
}