## basic questions on RxD

Extending NEURON to handle reaction-diffusion problems.

Moderators: hines, wwlytton, ramcdougal

Tuoma
Posts: 9
Joined: Mon Apr 23, 2018 6:59 am

### basic questions on RxD

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 https://neuron.yale.edu/neuron/static/d ... 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.
ramcdougal
Posts: 215
Joined: Fri Nov 28, 2008 3:38 pm
Location: Yale School of Public Health

### Re: basic questions on RxD

(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()

t.record(h._ref_t)
ca_val.record(ca.nodes(dend)(0.5)._ref_concentration)
buffer_val.record(buffer.nodes(dend)(0.5)._ref_concentration)
cabuf_val.record(cabuf.nodes(dend)(0.5)._ref_concentration)

h.finitialize(-65)
h.continuerun(10)

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

``````

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)

}

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

{
// define the graph
g = new Graph()
g.size(0, 10, 0, 1)
graphList.append(g)
}

{
// run the simulation
tstop = 20
run()
}
``````
(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 https://github.com/neuronsimulator/nrn/ ... 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 psection.py 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
psection(soma)
``````
And here is the psection.py 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())):
mt.select(i)
mypps = set()
pp = mt.pp_begin(sec=sec)
while pp is not None:
pp = mt.pp_next()
if mypps:
mt.selected(mname)
pps[mname] = 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())):
mt.select(i)
mt.selected(mname)
name = mname
if name in center_seg_dir:
mechs_present.append(name)

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(ms.name(mname, j))
name = mname
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]
else:
mechname = name#+ '_' + mech
for seg in sec:
if n > 1:
pvals.append([getattr(seg, mechname)[i] for i in range(n)])
else:
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
else:
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'] = [seg.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):
my_species.append(sp)
results['species'] = set(my_species)
results['name'] = sec.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.insert('hh')
soma.insert('pas')
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)
h.finitialize(-65)
pprint(psection(soma))
print('\n')
pprint(psection(dend))
``````
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
h.cvode.re_init()

t, ip3s = h.Vector(), h.Vector()
t.record(h._ref_t)
ip3s.record(ip3.nodes._ref_concentration)

h.finitialize(-65)

# 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))

h.continuerun(10)

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

``````
(4) The units should be mmol/litre/ms.
SurG
Posts: 11
Joined: Fri Jun 19, 2015 4:39 pm

### Re: basic questions on RxD

Hi,
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
Posts: 42
Joined: Tue Apr 18, 2017 10:05 am

### Re: basic questions on RxD

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;

https://neuron.yale.edu/neuron/docs/ini ... 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)
``````
Last edited by adamjhn on Sat Jul 06, 2019 10:59 am, edited 1 time in total.
SurG
Posts: 11
Joined: Fri Jun 19, 2015 4:39 pm

### Re: basic questions on RxD

This one still gives the same error "no such region"
I'm using Python 3.7 on a Windows 64bit system
Posts: 42
Joined: Tue Apr 18, 2017 10:05 am

### Re: basic questions on RxD

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)
``````
SurG
Posts: 11
Joined: Fri Jun 19, 2015 4:39 pm

### Re: basic questions on RxD

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.