Strange NEURON behavior with multithreading

When Python is the interpreter, what is a good
design for the interface to the basic NEURON
concepts.

Moderator: hines

rth
Posts: 41
Joined: Thu Jun 21, 2012 4:47 pm

Strange NEURON behavior with multithreading

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

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
}``````
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.

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()

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()
``````
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

Code: Select all

``````hpc = h.ParallelContext()