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

Post Reply
rth
Posts: 38
Joined: Thu Jun 21, 2012 4:47 pm

Strange NEURON behavior with multithreading

Post by rth » 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

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()
#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()
If I run it as a single thread the result is very reasonable.
Image
(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)
and run it on 8 cores, it produces a very strange result.
Image
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

Post Reply