basic questions on RxD

Extending NEURON to handle reaction-diffusion problems.

Moderators: hines, wwlytton, ramcdougal

Post Reply
Posts: 8
Joined: Mon Apr 23, 2018 6:59 am

basic questions on RxD

Post by Tuoma » Tue Apr 24, 2018 10:10 am

I have some questions to which I didn't find answers in the tutorial.

1) How can I record concentrations? Say, I would like to record the hydrogen concentration at location dend(0.5) into a vector in the example of ... ction-rate. Is there both a Python and a HOC way of doing that?

2) Is there a way to see which species are inserted into which section in a similar fashion as with 'psection()' in HOC? And is there a command to list all the reactions?

3) How can I introduce molecular flux into the intracellular compartments, e.g. a constant flux x mM/ms of hydrogen during certain time window? I noticed I can make a Ca2+ flux e.g. by
"reaction_Ca_flux = rxd.Reaction(Ca > 2*Ca, kf, mass_action = False)" and changing the kf in the middle of the simulation, but it looks a bit messy - I wonder if there's a better way of doing it?

4) When I set mass_action=False, are the units really molecules/μm3/ms as stated in the Frontiers article or mmol/litre/ms? The latter seems to fit better my test simulations.

Posts: 144
Joined: Fri Nov 28, 2008 3:38 pm
Location: Yale School of Medicine

Re: basic questions on RxD

Post by ramcdougal » Tue Apr 24, 2018 10:53 am

(1) To get a concentration pointer in Python, use ._concentration_ref; you can pass that to Vector.record. An example follows:

Code: Select all

from neuron import h, rxd
from matplotlib import pyplot


dend = h.Section(name='dend')
cyt = rxd.Region([dend], name='cyt', nrn_region='i')
ca = rxd.Species(cyt, name='ca', charge=2, initial=1)
buffer = rxd.Species(cyt, name='buffer', initial=1)
cabuf = rxd.Species(cyt, name='cabuf', initial=0)

buffering = rxd.Reaction(2 * ca + buffer > cabuf, 1)

t = h.Vector()
ca_val = h.Vector()
buffer_val = h.Vector()
cabuf_val = h.Vector()



pyplot.plot(t, ca_val, label='ca')
pyplot.plot(t, buffer_val, label='buffer')
pyplot.plot(t, cabuf_val, label='cabuf')

For HOC: If you have a named species that lives on a place where nrn_region='i', then you can use the regular & syntax; see HOC code below. Note that this still requires HOC to be able to connect to Python, which with recent versions of NEURON depends on an environment variable to specify which Python to use. To test if the connection works, run nrniv and type nrnpython("print('hi')") If the connection was successful, it will display hi and return 1; otherwise it returns 0.

Code: Select all

objref pyobj, h, rxd, cyt, ca, buf, cabuf, buffering, g

    // load reaction-diffusion support and get convenient handles
    nrnpython("from neuron import h, rxd")
    pyobj = new PythonObject()
    rxd = pyobj.rxd
    h = pyobj.h

    // define the domain and the dynamics
    create soma
    cyt = rxd.Region(h.allsec(), "i")
    ca = rxd.Species(cyt, 0, "ca", 2, 1)
    buf = rxd.Species(cyt, 0, "buf", 0, 1)
    cabuf = rxd.Species(cyt, 0, "cabuf", 0, 0)

    buffering = rxd.Reaction(ca.__add__(ca).__add__(buf), cabuf, 1, 0)

    // if launched with nrniv, we need this to get graph to update automatically
    // and to use run()

    // define the graph
    g = new Graph()
    g.addvar("ca", &soma.cai(0.5), 1, 1)
    g.addvar("cabuf", &soma.cabufi(0.5), 2, 1)
    g.size(0, 10, 0, 1)

    // run the simulation
    tstop = 20
(2) In the released versions, there is no psection equivalent that displays the reaction-diffusion information, but we have added it to the development version; see ... b659eccab4

Using the development version, you can just do a sec.psection() and get a dictionary that you can mine to identify reaction-diffusion information and more. You can get similar functionality in the released versions and older version using the code below. You can run it directly to get an example of what it can do (look at the if __name__ == '__main__' bit). From your own code or the terminal you can:

Code: Select all

import pprint
from psection import psection

# the following assumes you have a section defined called soma
And here is the file:

Code: Select all

def psection(sec):
    from neuron import h, rxd
    from neuron.rxd import region, species
    from neuron.rxd import rxd as rxd_module

    results = {}

    mname = h.ref('')

    # point processes
    pps = {}
    mt = h.MechanismType(1)
    for i in range(int(mt.count())):
        mypps = set()
        pp = mt.pp_begin(sec=sec)
        while pp is not None:
            pp = mt.pp_next()
        if mypps:
            pps[mname[0]] = mypps
    results['point_processes'] = pps

    center_seg_dir = dir(sec(0.5))

    mechs_present = []

    # membrane mechanisms
    mt = h.MechanismType(0)

    for i in range(int(mt.count())):
        name = mname[0]
        if name in center_seg_dir:

    results['density_mechs'] = {}
    results['ions'] = {}

    for mech in mechs_present:
        my_results = {}
        ms = h.MechanismStandard(mech, 0)
        for j in range(int(ms.count())):
            n = int(, j))
            name = mname[0]
            pvals = []
            # TODO: technically this is assuming everything that ends with _ion
            #       is an ion. Check this.
            if mech.endswith('_ion'):
                pvals = [getattr(seg, name) for seg in sec]
                mechname = name#+ '_' + mech
                for seg in sec:
                    if n > 1:
                        pvals.append([getattr(seg, mechname)[i] for i in range(n)])
                        pvals.append(getattr(seg, mechname))
            my_results[name[:-(len(mech) + 1)] if name.endswith('_' + mech) else name] = pvals
        # TODO: should be a better way of testing if an ion
        if mech.endswith('_ion'):
            results['ions'][mech[:-4]] = my_results
            results['density_mechs'][mech] = my_results

    morphology = {'L': sec.L,
                  'diam': [seg.diam for seg in sec],
                  'pts3d': [(sec.x3d(i), sec.y3d(i), sec.z3d(i), sec.diam3d(i)) for i in range(sec.n3d())]}

    results['morphology'] = morphology
    results['nseg'] = sec.nseg
    results['Ra'] = sec.Ra
    results['cm'] = [ for seg in sec]

    regions = {r() for r in region._all_regions if r() is not None and sec in r().secs}
    results['regions'] = regions

    my_species = []
    for sp in species._all_species:
        sp = sp()
        if sp is not None:
            sp_regions = sp._regions
            if not hasattr(sp_regions, '__len__'):
                sp_regions = [sp_regions]
            if any(r in sp_regions for r in regions):
    results['species'] = set(my_species)
    results['name'] =
    results['hoc_internal_name'] = sec.hoc_internal_name()
    results['cell'] = sec.cell()

    #all_active_reactions = [r() for r in rxd_module._all_reactions if r() is not None]

    return results

if __name__ == '__main__':
    from pprint import pprint
    from neuron import h, rxd

    soma = h.Section(name='soma')
    soma.nseg = 3
    iclamp = h.IClamp(soma(0.3))
    gaba = h.ExpSyn(soma(0.7))
    dend = h.Section(name='dend')
    cyt = rxd.Region([dend], name='cyt', nrn_region='i')
    er = rxd.Region([dend], name='er')
    na = rxd.Species(cyt, name='na', charge=1)
    ca = rxd.Species([cyt, er], name='ca', charge=2)
One minor comment: the psection definition above and the version in the repository at the immediate moment does not return active reactions, but if you look at the commented out line about all_active_regions you'll see how to get everything. Why doesn't it do everything? There's a minor bug that interferes with printing the reactions if they have constant rates. In the released versions this affects all such reactions; in the development version, it only affects MultiCompartmentReactions, but if you look at the pull requests, you'll see the patch for it, and we'll merge this in soon.

(3) Use rxd.Rate for constant rates, not rxd.Reaction… rxd.Rate basically just adds whatever argument you give it to the derivative. e.g.

Code: Select all

ca_flux = rxd.Rate(ca, 1)
causes ca to increase at a rate of 1 mM per ms in all compartments. One way of doing this is demonstrated below:

Code: Select all

from neuron import h, rxd
from matplotlib import pyplot


dend = h.Section(name='dend')
cyt = rxd.Region([dend])
ip3 = rxd.Species(cyt)

k = rxd.Parameter(cyt, initial=1)
rxd_prod = rxd.Rate(ip3, k)

def set_param(param, val):
    param.nodes.value = val
    # the below is needed to work with variable step
    # future versions of NEURON should not need this

t, ip3s = h.Vector(), h.Vector()


# define parameter change points at t=4 and t=8
h.cvode.event(4, lambda: set_param(k, 3))
h.cvode.event(8, lambda: set_param(k, -2))


pyplot.plot(t, ip3s)
pyplot.ylabel('[IP3] (mM)')
pyplot.xlabel('t (ms)')

(4) The units should be mmol/litre/ms.

Post Reply