Strange NEURON behavior with multithreading
Posted: Tue Dec 11, 2018 11:53 am
Dear NEURON developers,
I have run into a very strange NEURON behavior, which I couldn't understand. It relates to the multithreading, but perhaps I did something stupid.
I have a simple mod mechanisms
So the code runs two lists of neurons with different parameter sets to check the difference in F/I curves. Each neuron in a list has a particular Iapp.
If I run it as a single thread the result is very reasonable.
(images are keeping disappear, please follow the link to see it https://drive.google.com/file/d/18o8OtB ... rk6s-/view)
however, if I uncomment lines 139 and 140
and run it on 8 cores, it produces a very strange result.
https://drive.google.com/file/d/1K5rRVw ... sp=sharing
I have tried to turn on/off TABLES, and GLOBALS for the tabulated variables, set/remove THREADSAFE in the mod file with no luck.
What did I do wrong?
-rth
I have run into a very strange NEURON behavior, which I couldn't understand. It relates to the multithreading, but perhaps I did something stupid.
I have a simple mod mechanisms
Code: Select all
UNITS {
(mA) = (milliamp)
(mV) = (millivolt)
(mS) = (millisiemens)
}
NEURON {
SUFFIX nak3dv02
NONSPECIFIC_CURRENT i
RANGE ninit, hinit : initial conditions for n and h
: if negative, it uses stady-state for given
: voltage
RANGE n0,nv12,nsl,nt0,nst,nv0,nsg
RANGE h0,hv12,hsl,ht0,hst,hv0,hsg
RANGE gk,ek,gna,ena,gl,el,a,b
:GLOBAL minf, ninf, ntau
}
PARAMETER {
v (mV)
gna = 120. (mS/cm2)
ena = 50. (mV)
gk = 36. (mS/cm2)
ek = -77. (mV)
gl (mS/cm2)
el (mV)
n0 = 0 (1) : foot of n
nsc = 1 (1) : scaler of n
nv12 (mV) : v12 for n
nsl (mV) : slope of n
nt0 (ms) : foot for tau_n
nst (ms) : scaler for tau_n
nv0 (mV) : center of max for tau_n
nsg (mV) : sigma for tau_n
h0 = 0 (1) : foot of h
hsc = 1 (1) : scaler of h
hv12 (mV) : v12 for h
hsl (mV) : slope of h
ht0 (ms) : foot for tau_h
hst (ms) : scaler for tau_h
hv0 (mV) : center of max for tau_h
hsg (mV) : sigma for tau_h
a (1) : intersect h/n
b (1) : slope h/n
ninit = -1.
hinit = -1.
}
STATE {
h
n
}
ASSIGNED {
i (mA/cm2)
minf
hinf
htau (ms)
ninf
ntau (ms)
}
BREAKPOINT {
SOLVE states METHOD cnexp
:----vvvv-- is needed to convert uA/cm2 to mA/cm2
i = (1e-3)*( gna*minf*minf*minf*(a+b*n)*h*(v-ena)+gk*n*n*n*n*(v-ek)+gl*(v-el) )
}
DERIVATIVE states {
rates(v)
h'= (hinf - h) / htau
n'= (ninf - n) / ntau
}
INITIAL {
rates(v)
if (ninit < 0 || ninit > 1){
n = ninf
} else {
n = ninit
}
if (hinit < 0 || hinit > 1){
h = hinf
} else {
h = hinit
}
nsc = 1 - n0
hsc = 1 - h0
}
PROCEDURE rates(v (mV)) {
UNITSOFF
:TABLE minf, hinf, htua, ninf, ntau FROM -100 TO 100 WITH 200
minf = 1./(1.+exp(-(v+40.)/9.5))
ninf = n0 + nsc/(1.+exp(-(v-nv12)/nsl ))
ntau = nt0 + nst*exp(-((v-nv0)/nsg)*((v-nv0)/nsg))
hinf = h0 + hsc/(1.+exp(-(v-hv12)/hsl ))
htau = ht0 + hst*exp(-((v-hv0)/hsg)*((v-hv0)/hsg))
UNITSON
}Code: Select all
import os,sys
from numpy import *
from matplotlib.pyplot import *
#== The main parameter set ==#
gna , ena = 120. , 50.0
gk , ek = 36.0, -77.0
gl, el = 0.1, -65.
hslp, hitr = -1.10692947808, 0.906483183915
n0 = 0.
nv12,nsl = -54.5, 19.
nt0, nsc = 0.5, 5.
nv0, nsg = -60., 30.
h0 = 1.
hv12,hsl = -55., -9.
ht0, hsc = 50., 75.
hv0, hsg = -60., 30.
#== An alternative parameter set ==#
Agna , Aena = 120., 50.0
Agk , Aek = 36.0, -77.
Agl, Ael = 0.1, -65.
Ahslp, Ahitr = -1.10692947808, 0.906483183915
An0 = 0.
Anv12,Ansl = -54.5, 19.
Ant0, Ansc = 0.5, 5.
Anv0, Ansg = -60., 30.
Ah0 = 0.
Ahv12,Ahsl = -55., -9.
Aht0, Ahsc = 50., 75.
Ahv0, Ahsg = -60., 30.
os.system("nrnivmodl")
from neuron import h
h.dt = 0.01
class neuron:
def __init__ (self, vinit=None, ninit=None, hinit=None, **params):
self.soma = h.Section()
self.soma.L = 100.
self.soma.diam = 10./pi
self.soma.insert('nak3dv02')
#=== neuron initialization ===#
if type(vinit) is float or type(vinit) is int:
self.soma(0.5).v = vinit
if type(ninit) is float or type(ninit) is int:
self.soma(0.5).nak3dv02.ninit = ninit
if type(hinit) is float or type(hinit) is int:
self.soma(0.5).nak3dv02.hinit = hinit
if len(params) != 0:
for n,v in params.items():
if n=="vinit" or n=="ninit" or n=="hinit": continue
exec "self.soma(0.5).nak3dv02.{}={}".format(n,v)
#=== set recordings ===#
self.v = h.Vector()
self.v.record(self.soma(0.5)._ref_v)
self.n = h.Vector()
self.n.record(self.soma(0.5).nak3dv02._ref_n)
self.h = h.Vector()
self.h.record(self.soma(0.5).nak3dv02._ref_h)
#=== record spike times ===#
self.spks = h.Vector()
self.sptr = h.APCount(.5, sec=self.soma)
self.sptr.thresh = 0.#25
self.sptr.record(self.spks)
#=== synapse ===#
self.isyn = h.Exp2Syn(0.5,sec=self.soma)
self.isyn.tau1 = 1.
self.isyn.tau2 = 3.
self.isyn.e = -75.
#=== Input current ===#
self.I = h.IClamp(0.5, sec=self.soma)
self.I.delay = 0
self.I.amp = 0.
self.I.dur = 1e9
def create_T2():
return neuron( gna=gna,ena=ena,gk=gk,ek=ek,gl=gl,el=el,a=hitr,b=hslp,
n0=n0,nv12=nv12,nsl=nsl,nt0=nt0,nst=nsc,nv0=nv0,nsg=nsg,
h0=h0,hv12=hv12,hsl=hsl,ht0=ht0,hst=hsc,hv0=hv0,hsg=hsg,
vinit = -67.91262149643644, ninit = 0.32971471805616)
def create_AT2():
return neuron( gna=Agna,ena=Aena,gk=Agk,ek=Aek,gl=Agl,el=Ael,a=Ahitr,b=Ahslp,
n0=An0,nv12=Anv12,nsl=Ansl,nt0=Ant0,nst=Ansc,nv0=Anv0,nsg=Ansg,
h0=Ah0,hv12=Ahv12,hsl=Ahsl,ht0=Aht0,hst=Ahsc,hv0=Ahv0,hsg=Ahsg,
vinit = -67.91262149643644, ninit = 0.32971471805616)
Fast = True
N,dur,Imin,Imax = 50,2000., 0., 0.5
#hpc = h.ParallelContext()
#hpc.nthread(8)
T1n,T2u = [], []
for i,Iu in enumerate( linspace(Imin,Imax ,N+1) ):
n1,up = create_T2(), create_AT2()
n1.Ia, n1.It = h.Vector(),h.Vector()
up.Ia, up.It = h.Vector(),h.Vector()
n1.It, up.It = \
n1.It.from_python( [0.,dur,dur*2.1] ),\
up.It.from_python( [0.,dur,dur*2.1] )
n1.Ia, up.Ia = \
n1.Ia.from_python( [Imin,Iu,Iu] ),\
up.Ia.from_python( [Imin,Iu,Iu] )
n1.Ia.play(n1.I._ref_amp,n1.It,0 )
up.Ia.play(up.I._ref_amp,up.It,0 )
n1._i_app = Iu
up._i_app = Iu
T1n.append( n1 )
T2u.append( up )
tmax = dur*2.1
t = h.Vector()
t.record(h._ref_t)
h.finitialize()
h.fcurrent()
h.frecord_init()
while h.t < tmax:h.fadvance()
T1FI, T2FI = [], []
for n1,up in zip(T1n,T2u):
n1sp, upsp = array(n1.spks),array(up.spks)
n1sp, upsp = n1sp[where(n1sp > dur*1.1)], upsp[where(upsp > dur*1.1)]
T1FI.append( (n1._i_app, float(n1sp.shape[0])*1000./dur))
T2FI.append( (up._i_app, float(upsp.shape[0])*1000./dur))
del n1,up
figure(2,figsize=(11,11))
plot(T1FI[:,0],T1FI[:,1],"o",label=r"The main parameter set")
plot(T2FI[:,0],T2FI[:,1],"x",label=r"An alternative parameter set")
legend(loc=0,fontsize=16)
ylabel("Frequency (Hz)")
xlabel(r"Current (uA/cm$^2$)")
savefig("FIcurve.png")
show()
(images are keeping disappear, please follow the link to see it https://drive.google.com/file/d/18o8OtB ... rk6s-/view)
however, if I uncomment lines 139 and 140
Code: Select all
hpc = h.ParallelContext()
hpc.nthread(8)
https://drive.google.com/file/d/1K5rRVw ... sp=sharing
I have tried to turn on/off TABLES, and GLOBALS for the tabulated variables, set/remove THREADSAFE in the mod file with no luck.
What did I do wrong?
-rth