MultiCompartmentReaction Section List Mismatch

Extending NEURON to handle reaction-diffusion problems.

Moderators: hines, wwlytton, ramcdougal

Post Reply
Posts: 5
Joined: Wed Jul 31, 2019 5:14 pm

MultiCompartmentReaction Section List Mismatch

Post by amerg »

I've been running into an issue and I'm curious what sort of work arounds may be available. My model previously had all sections with a region that was the cytosol and a region that was the endoplasmic reticulum (er). Now I'm making synaptic spines which don't always have ER, so I tried making regions where the cytosol region had more sections than the er region. But as soon as I tried to set up a vector to record a value I got this error

Code: Select all

raise RxDException("Rates and Reactions with species defined on different regions are not currently supported.") 
I made a simple test script and was able to replicate this error when the two regions don't completely overlap in terms of sections and there is a MultiCompartmentReaction between the two regions.

Is there a work around for this? I could make a separate region for the synaptic spine sections, but how would I implement axial diffusion between two regions? Is there a way to connect regions that I'm overlooking?
Posts: 43
Joined: Tue Apr 18, 2017 10:05 am

Re: MultiCompartmentReaction Section List Mismatch

Post by adamjhn »

The current work around is to create the regions everywhere as you had before, then alter the diffusion coefficient and use parameters to exclude the unwanted region. Unfortunately, this approach does not allow the volume fraction to vary; that is, if there is no ER in the spine, you might expect cytosol to occupy a different fraction of the volume, but this approach would still have it being the same percent as in the dendrite. For example;

Code: Select all

from neuron import h, rxd
from neuron.units import μM
from matplotlib import pyplot as plt

# parameters (from
fraction_cyt = 0.83
fraction_er = 0.17
caD = 0.233
gleak = 3 * 1000 # increased leak rate for this example

# morphology
dend = h.Section(name='dend')
dend.nseg = 101
spine = h.Section(name='spine')

# regions
cyt = rxd.Region(h.allsec(), name='cyt', nrn_region='i',
er = rxd.Region(h.allsec(), name='er',
cyt_er_membrane = rxd.Region(h.allsec(), geometry=rxd.ScalableBorder(1))

# helper function
def excluded_er(node):
    if node.region == er and node.sec == spine:
        return 0
    return 1

# species
ca = rxd.Species([cyt,er], name='ca', charge=2, 
                 initial=lambda nd: 100 * μM if nd.region == cyt and 
                                    nd.sec == spine else 60*μM)

dendonly = rxd.Parameter([er], name='dendonly',
                          initial=lambda nd: excluded_er(nd))

# prevent spine ER diffusion
for nd in ca.nodes:
    nd.d = caD * excluded_er(nd) 

# reactions
leak = rxd.MultiCompartmentReaction(ca[er] + dendonly[er],
                                    ca[cyt] + dendonly[er],
                                    gleak, gleak,

# record some concentrations
t_vec = h.Vector().record(h._ref_t)
caer_vec = h.Vector().record(ca[er].nodes(dend(0.5))._ref_concentration)
cacyt_vec = h.Vector().record(spine(0.5)._ref_cai)

# run the simulation

# plot the results
plt.plot(t_vec, cacyt_vec/μM, label='Ca$^{2+}_{Cyt}$ spine')
plt.plot(t_vec, caer_vec/μM, label='Ca$^{2+}_{ER}$ dendrite')
plt.xlabel('time (ms)')
plt.ylabel('concentration (μM)')

plt.plot([seg.x * dend.L for seg in dend], [seg.cai/μM for seg in dend],
plt.plot([nd.x * dend.L for nd in ca[er].nodes(dend)], 
         [nd.value/μM for nd in ca[er].nodes(dend)], label='ER')
plt.xlabel('dendrite (μm)')
plt.ylabel('concentration (μM)')

If you wish, you may experiment with this code example online: ... aDXIu0TcCp
Post Reply