I'm getting some Segmentation fault errors after using the node.include_flux in a bit longer script that simulates a number of presynaptic spines and postsynaptic sections and random connections between them. This may be a technical issue somewhere else in my code but I haven't found the cause despite many attempts. Sorry for a long MWE:
Code: Select all
from neuron import h, rxd
import sys
from pylab import *
tstop = 4000
Nspines = 20 #Number of presynaptic spines
NsynTot = 1000 #Number of connections between spines and postsynaptic dendritic sections
NTdiffusion = 1e-17 #Diffusion parameter for neurotransmitter between pre and postsynaptic sections
if len(sys.argv) > 1:
Nspines = int(float(sys.argv[1]))
if len(sys.argv) > 2:
NsynTot= int(float(sys.argv[2]))
Nsyn_exc_per_spine = int(NsynTot/Nspines) #Number of Exc. synaptic inputs to be placed for each spine
seed(1)
#Initialize the post-synaptic dendritic sections and the region containing 84 post-synaptic dendritic sections:
h("""
load_file("stdlib.hoc")
load_file("stdrun.hoc")
objref cvode
cvode = new CVode()
cvode.active(1)
cvode.atol(1e-8)
create dend[84]
for (i=0;i<84;i+=1) { dend[i] {L=1 diam=1 nseg=1 } }
tstop = """+str(tstop)+"""
""")
seclist = []
for sec in h.dend:
seclist.append(sec)
cyt_postsyn_exc = rxd.Region(seclist, name="postsyn", nrn_region='i')
#Pick locations where synapses will placed in the post-synaptic sections based on their diameter (here all diameters equal):
isecs = [] #This will have the indices (to secnames) where the synapses will be located
isegs = [] #This will have the indices (0 to nseg-1) where the synapses will be located
xsecs = [] #This will have the x along the section where the synapses will be located
for isyn in range(0,Nsyn_exc_per_spine*Nspines):
myi = int(84*rand())
secname = "dend["+str(int(myi))+"]"
isecname = [i for i in range(0,len(seclist)) if secname == seclist[i].name()][0]
h(secname+" thisnseg = nseg")
isecs.append(isecname)
#Initialize presynaptic spines and the region containing the presynaptic spines
h("create spine["+str(Nspines)+"]")
h("for (ispine=0;ispine<"+str(Nspines)+";ispine+=1) { spine[ispine] { L = 0.8 diam = 0.4 nseg = 1 } }")
spinecyt = rxd.Region(h.spine, name='cyt', nrn_region='i')
#Define our only species, the neurotransmitter (NT) that will be present in all pre- and postsynaptic sections, and the only reaction (decay of NT)
NT = rxd.Species([cyt_postsyn_exc,spinecyt], name='NT', charge=0, initial=0.0, atolscale=1.0)
myreaction = rxd.Reaction(NT, NT*0, 10.0)
#Save the NT nodes into a list and save the index array for each synapse so that iNTnodes[isyn] will tell which index of NTnodes is the section to which the synapse projects
NTnodes = NT.nodes(cyt_postsyn_exc)
iNTnodes = []
for isyn in range(0,len(isecs)):
iNTnodes.append(isecs[isyn])
#Define NT fluxes between pre- and post-synaptic sections. Each of the spines (Nspines) is connected to a fixed number (Nsyn_exc_per_spine) of randomly picked post-synaptic sections/segments
for ispine in range(0,Nspines):
prend = NT.nodes(spinecyt)[ispine]
for isyn in range(0,Nsyn_exc_per_spine):
postnd = NTnodes[iNTnodes[ispine*Nsyn_exc_per_spine+isyn]]
print('Adding flux to exc. synapse '+str(isyn+1)+"/"+str(Nsyn_exc_per_spine)+" in spine "+str(ispine+1)+"/"+str(Nspines)+": prend seg = "+str(prend.segment)+", postnd seg = "+str(postnd.segment))
postnd.include_flux(lambda postnd = postnd, prend = prend:(prend.value-postnd.value)*postnd.volume*NTdiffusion, units='mol/ms')
prend.include_flux(lambda postnd = postnd, prend = prend:(postnd.value-prend.value)*prend.volume*NTdiffusion, units='mol/ms')
#Introduce flux parameter and a NT rate to each presynaptic spine whose value depends on the parameter
NTflux_rates = []
NTreaction_fluxes = []
for ispine in range(0,Nspines):
NTflux_rates.append(rxd.Parameter(spinecyt, initial=0))
NTreaction_fluxes.append(rxd.Rate(NT, NTflux_rates[-1]))
NTrec_pre = h.Vector()
NTrec_post = h.Vector()
vec_t = h.Vector()
NTrec_pre.record(NT.nodes(h.spine[0])(0.5)[0]._ref_concentration)
NTrec_post.record(NT.nodes(h.dend[0])(0.5)[0]._ref_concentration)
vec_t.record(h._ref_t)
h("""v_init = -80""")
h.init()
def set_param(param, val):
param.nodes.value = val
h.cvode.re_init()
for ispine in range(0,Nspines):
for ipst in range(0,int(tstop/300)):
h.cvode.event(300*ipst+(ispine+1)*300/(1+Nspines), lambda ispine=ispine: set_param(NTflux_rates[ispine], 2.0)) #Increase the presynaptic NT concentration by 2.0 mM/ms for 5 ms every 300 ms
h.cvode.event(300*ipst+(ispine+1)*300/(1+Nspines)+5, lambda ispine=ispine: set_param(NTflux_rates[ispine], 0.0))
print("Starting simulation")
h.continuerun(tstop)
f,axarr = subplots(2,1)
axarr[0].plot(array(vec_t),array(NTrec_pre))
axarr[1].plot(array(vec_t),array(NTrec_post))
f.savefig("test.eps")
The code runs if I use a small number of synaptic connections, e.g.
python3 mycode.py 2 10
However, when I use a large number of synaptic connections I get a segmentation fault at h.init(), for example
python3 mycode.py 2 1000
The reason I thought there's some error related to my include_flux formalism is that if I comment out the two lines with include_flux() then the code is running fine.
Any ideas what could be wrong?