Page 1 of 1

basic questions on RxD

Posted: Tue Apr 24, 2018 10:10 am
by Tuoma
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.

Re: basic questions on RxD

Posted: Tue Apr 24, 2018 10:53 am
by ramcdougal
(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.

Re: basic questions on RxD

Posted: Sat Jul 06, 2019 9:54 am
by SurG
Is it possible to initialise the calcium concentrations of both cytosol and er to separate values?
From what I understand, the following statement initializes the concentrations of ca[cyt] and ca[er] to the same value, cac_init
ca = rxd.Species([cyt, er], d=caDiff, name='ca' ,charge=2,initial=cac_init)
But I would like to initialise ca[cyt] to 100 nM and ca[er] to 70 uM. I tried splitting the above command but I got an error "RxDException: no such region"

Any help will be appreciated.
Thank you

Re: basic questions on RxD

Posted: Sat Jul 06, 2019 10:38 am
by adamjhn
The initial keyword can take a function with a rxd node as an argument. This allows initial concentrations to depended on the location and region. Some examples are here; ... strategies

In your case, I think you want to define calcium like this;

Code: Select all

from neuron.units import nM, uM
ca = rxd.Species(cyt, name='ca', charge=2, initial=lambda nd: 100 * nM if nd.region == cyt else 70 * uM)

Re: basic questions on RxD

Posted: Sat Jul 06, 2019 10:54 am
by SurG
This one still gives the same error "no such region"
I'm using Python 3.7 on a Windows 64bit system

Re: basic questions on RxD

Posted: Sat Jul 06, 2019 2:21 pm
by adamjhn
Sorry, of course you still need to define it on both the cyt and the er;

Code: Select all

from neuron.units import nM, uM
ca = rxd.Species([cyt, er], name='ca', charge=2, initial=lambda nd: 100 * nM if nd.region == cyt else 70 * uM)

Re: basic questions on RxD

Posted: Sat Jul 06, 2019 3:12 pm
by SurG
This worked. Thank you so much!

I have another query. According to the RxD Tutorial 1.0.0 documentation, the unit of concentration is mM. However, in the same tutorial, "Calcium Waves in RxD", even though cac_init is initialised to 0.0001 mM (~100 nM) and kserca is initialised to 0.1 mM, ca[cyt] is multiplied by 1000 in the expression for SERCA. Why has this been done? Are the default units of ca[cyt] (and ca[er]) in uM?
I'm actually building a model with .mod mechanisms inserted into RxD. So I'm following the units of NEURON (NMODL), i.e., the fluxes are in mM/ms and the concentrations are in mM. But I am not able to understand why ca[cyt] is multiplied by 1000 but kserca is not.

Re: basic questions on RxD

Posted: Mon Jul 08, 2019 10:35 am
by adamjhn
The parameters used in the tutorial are python floats, NEURON does not assign them any particular units. This makes the example somewhat opaque and is part of the reason why in NEURON 7.7 the units modules was introduced.

In this specific example you can infer that kserca must have units μM as ca[cyt] is multiple by 1000 giving its concentration in units μM. Then, as it is a multicompartment reaction gserca must have units of molecules per μm² per ms.