(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
h.load_file('stdrun.hoc')
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)[0]._ref_concentration)
buffer_val.record(buffer.nodes(dend)(0.5)[0]._ref_concentration)
cabuf_val.record(cabuf.nodes(dend)(0.5)[0]._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)
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()
load_file("stdrun.hoc")
}
{
// 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)
graphList[0].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:
mypps.add(pp)
pp = mt.pp_next()
if mypps:
mt.selected(mname)
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())):
mt.select(i)
mt.selected(mname)
name = mname[0]
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[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]
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.
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
h.load_file('stdrun.hoc')
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[0]._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.