MultiCompartmentReaction across sections

[Tuoma]

Is it possible to make a MultiCompartmentReaction that involves species in non-connected sections?

I get an error "RxDException: Multicompartment reactions the membrane and sources and destinations must share common sections." which suggests this is not possible. Here's a MWE where I would like the neurotransmitter diffuse from a presynaptic spine to postsynaptic spine:

Code: Select all

from neuron import h,rxd

h("""
load_file("stdlib.hoc")
load_file("stdrun.hoc")
create presyn, postsyn
presyn {nseg=1 L=1 diam=0.5}
postsyn {nseg=1 L=1 diam=0.5}
tstop = 100
""")

cytpre = rxd.Region([h.presyn], name='presyn', nrn_region='i')
cytpost = rxd.Region([h.postsyn], name='postsyn', nrn_region='i')
NT = rxd.Species([cytpre,cytpost], name='NT', charge=0, initial=0.1)

conceptualMembrane = rxd.Region([h.presyn,h.postsyn], name='mem', geometry=rxd.membrane())

NTdiffusion = rxd.MultiCompartmentReaction(NT[cytpre],NT[cytpost], 1.0, membrane=conceptualMembrane, scale_by_area=False)

h.init()
h.continuerun(h.tstop)

The nrn_region was 'i' in these species but the message doesn't change even if it's 'o'.

[adamjhn]

You are right about multi-compartment reactions, they are used for reactions between regions of a section or with the extracellular space.
A possible solution to the problem of diffusion between a pre- and post-synaptic section is to use the node.include_flux. For example; Diffusion across a synapse colab.

[Tuoma]

Thanks, that seems to work.

[Tuoma]

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?

[hines]

Code: Select all

python3 mycode.py 2 10

Gives me a segfault for NEURON version 8.1.0 but not for 8.2.2 or current master.

I also built the master using the address sanitizer and there were no substantive complaints (only some memory leaks), so there is a possibility that the bug was fixed after 8.1.0.

[Tuoma]

Thanks, with 8.2.2 I always get past the init() without segfault.

However, when I stop the simulation with CTR+C I always get a segfault. That wouldn't necessarily be a problem if my simulations always finished.

It seems the simulation is very sensitive to the NTdiffusion value, which is expected I guess - with large enough NTdiffusion values, which would lead to neurotransmitter concentration to be very similar between pre and post sides of the synaptic cleft, the simulation becomes super slow. This suggests that it's probably not advisable to model the neurotransmitter diffusion from presynaptic side to the postsynaptic side, at least this way, since it would imply that one has to use super small time steps to ensure correct integration. If you have alternative ideas to force the concentration to be the same or almost the same at pre and post sides I'd be happy to hear. I'm guessing the wisest thing to do is to use pointers in the .mod mechanisms at the postsynaptic side to the NT values at the presynaptic side.

[hines]

I've run the address sanitizer several times and upon ^C have always ended up at

Code: Select all

nrn-enable-sanitizer --preload python3 -i test.py
...
Adding flux to exc. synapse 500/500 in spine 2/2: prend seg = spine[1](0.5), postnd seg = dend[65](0.5)
Starting simulation
^CAddressSanitizer:DEADLYSIGNAL
=================================================================
==429635==ERROR: AddressSanitizer: UNKNOWN SIGNAL on unknown address 0x634000000844 (pc 0x7fd8a5ed7c18 bp 0x7ffc3ab68a20 sp 0x7ffc3ab68a10 T0)
    #0 0x7fd8a5ed7c18 in _Py_IS_TYPE /home/hines/.pyenv/versions/3.10.5/include/python3.10/object.h:146
    #1 0x7fd8a5ed7c4a in _PyObject_TypeCheck /home/hines/.pyenv/versions/3.10.5/include/python3.10/object.h:247
    #2 0x7fd8a5ed9abb in apply_node_flux(int, long*, double*, _object**, double, double*) ../src/nrnpython/rxd.cpp:327
 ...

I think, but am not certain, that an arbitrary break of the execution leaves the system in an unstable state. A possible workaround is to set the stoprun
variable. That is watched by cvode and stops integrating at safe points. At least that worked the one time I tried it. Try experimenting with

Code: Select all

diff --git a/test.py b/test.py
index 22fa87a..e490128 100644
--- a/test.py
+++ b/test.py
@@ -1,9 +1,14 @@
-from neuron import h, rxd
+from neuron import h, rxd, gui
+
+h.xpanel("stoprun")
+h.xbutton("Stop","stoprun=1")
+h.xpanel()
+
 import sys
 from pylab import *