rxd.nthread · rxd.re_init · rxd.set_solve_type
Basic Reaction-Diffusion
Overview
NEURON provides the rxd
submodule to simplify and standardize the specification of
models incorporating reaction-diffusion dynamics, including ion accumulation.
The interface is implemented using Python, however as long as Python is available to
NEURON, reaction-diffusion dynamics may be specified using HOC.
We can access the rxd
module in HOC via:
objref pyobj, h, rxd
{
// load reaction-diffusion support and get convenient handles
nrnpython("from neuron import h, rxd")
pyobj = new PythonObject()
rxd = pyobj.rxd
h = pyobj.h
}
The above additionally provides access to an object called h
which is traditionally
how Python accesses core NEURON functionality (e.g. in Python one would use h. Vector
instead of Vector
). You might not need to use h since when working in HOC,
but it does provide certain convenient functions like h.allsec()
, which returns
an iterable of all sections usable with rxd
without having to explicitly construct
a SectionList
.
The main gotchas of using rxd in HOC is that (1) rxd
in Python uses operator overloading to
specify reactants and products; in HOC, one must use __add__
, etc instead.
(2) rxd in Python is usually used with keyword arguments; in HOC, everything must be
specified using positional notation.
Here’s a full working example that simulates a calcium buffering reaction:
Ca + Buf <> CaBuf
:
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__(buf), cabuf, 1, 0.1)
}
{
// 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()
}
In particular, note that instead of ca + buf
one must write
ca.__add__(buf)
.
In general, a reaction-diffusion model specification involves answering three conceptual questions:
Where the dynamics are occurring (specified using an
rxd.Region
orrxd.Extracellular
)Who is involved (specified using an
rxd.Species
orrxd.State
)What the reactions are (specified using
rxd.Reaction
,rxd.Rate
, orrxd.MultiCompartmentReaction
)
Another key class is rxd.Parameter
for defining spatially varying parameters.
Integration options may be specified using rxd.set_solve_type()
.
Specifying the domain
NEURON provides two main classes for defining the domain where a given set of reaction-diffusion rules
applies: rxd.Region
and rxd.Extracellular
for intra- and extracellular domains,
respectively. Once defined, they are generally interchangeable in the specification of the species involved,
the reactions, etc. The exact shape of intracellular regions may be specified using any of a number of
geometries, but the default is to include the entire intracellular space.
Intracellular regions and regions in Frankenhauser-Hodgkin space
- class rxd.Region
Declares a conceptual intracellular region of the neuron.
Syntax:
r = rxd.Region(secs, nrn_region, geometry, dimension, dx, name)
In NEURON 7.4+,
secs
is optional at initial region declaration, but it must be specified before the reaction-diffusion model is instantiated.All arguments are optional, but all prior arguments must be specified. To use the default values for the prior arguments, specify their values as
pyobj.None
.Here:
secs
is aSectionList
or any Python iterable of sections (e.g.h.allsec()
)nrn_region
specifies the classic NEURON region associated with this object and must be either"i"
for the region just inside the plasma membrane,"o"
for the region just outside the plasma membrane orpyobj.None
for none of the above.name
is the name of the region (e.g.cyt
orer
); this has no effect on the simulation results but it is helpful for debuggingdx
deprecated; when specifyingname
pass inpyobj.None
heredimension
deprecated; when specifyingname
pass inpyobj.None
here
- property rxd.Region.nrn_region
Get or set the classic NEURON region associated with this object.
There are three possible values:
'i'
– just inside the plasma membrane'o'
– just outside the plasma membranepyobj.None
– none of the above
Setting requires NEURON 7.4+, and then only before the reaction-diffusion model is instantiated.
- property rxd.Region.secs
Get or set the sections associated with this region.
The sections may be expressed as a
SectionList
or as any Python iterable of sections.- Note: The return value is a copy of the internal section list; modifying
it will not change the Region.
Setting requires NEURON 7.4+ and allowed only before the reaction-diffusion model is instantiated.
- property rxd.Region.geometry
Get or set the geometry associated with this region.
Setting the geometry to
None
will cause it to default torxd.geometry.inside
.Added in NEURON 7.4. Setting allowed only before the reaction-diffusion model is instantiated.
- property rxd.Region.name
Get or set the Region’s name.
Added in NEURON 7.4.
For specifying the geometry
NEURON provides several built-in geometries for intracellular simulation that may be specified
by passing a geometry
argument to the rxd.Region
constructor. New region shapes may
be defined by deriving from neuron.rxd.geometry.RxDGeometry
.
See the Python programmer’s reference for more on rxd.inside
. rxd.membrane
,
rxd.DistributedBoundary
, rxd.FractionalVolume
, rxd.Shell
,
rxd.FixedCrossSection
, rxd.FixedPerimeter
, and
rxd.ScalableBorder
.
Extracellular regions
- class rxd.Extracellular
Declare a extracellular region to be simulated in 3D; unlike
rxd.Region
, this always describes the extracellular region.See the entry for
rxd.Extracellular
in the Python programmer’s reference for more information.
Defining proteins, ions, etc
Values that are distributed spatially on an rxd.Region
or rxd.Extracellular
may be defined using
an rxd.Species
if they represent things that change and diffuse, an rxd.State
if they’re in fixed locations but changeable
(e.g. gates in an IP3R), or an rxd.Parameter
if
they are just fixed values.
- class rxd.Species
Declare an ion/protein/etc that can react and diffuse.
See the entry for
rxd.Species
in the Python programmer’s reference for more information.
- class rxd.State
Like an
rxd.Species
but indicates the semantics of something that is not intended to diffuse.See the entry for
rxd.State
in the Python programmer’s reference for more information.
- class rxd.Parameter
Declares a parameter, that can be used in place of a
rxd.Species
, but unlikerxd.Species
a parameter will not change.See the entry for
rxd.Parameter
in the Python programmer’s reference for more information.
Defining reactions
NEURON provides three classes for specifying reaction dynamics: rxd.Reaction
for single-compartment (local)
reactions; rxd.MultiCompartmentReaction
for reactions spanning multiple compartments (e.g. a pump that
moves calcium from the cytosol into the ER changes concentration in two regions), and rxd.Rate
for
specifying changes to a state variable directly by an expression to be added to a differential equation.
Developers may introduce new forms of reaction specification by subclassing
neuron.rxd.generalizedReaction.GeneralizedReaction
, but this is not necessary for typical modeling usage.
It is sometimes necessary to build rate expressions including mathematical functions. To do so, use the
functions defined in neuron.rxd.rxdmath
as listed below.
- class rxd.Reaction
Reaction at a point. May be mass-action or defined via custom dynamics.
Syntax:
r1 = rxd.Reaction(reactant_sum, product_sum, forward_rate, backward_rate, regions, custom_dynamics)
backward_rate
,regions
, andcustom_dynamics
are optional, but when used from HOC, all previous parameters must be specified. To specify that the dynamics should be custom (i.e. fully defined by the rates) without abackward_rate
or restricting to specific regions, pass0
forbackward_rate
andpyobj.None
forregions
.Example:
// here: ca + buf <> cabuf, kf = 1, kb = 0.1 buffering = rxd.Reaction(ca.__add__(buf), cabuf, 1, 0.1)
Note the need to use
__add__
instead of+
. To avoid this cumbersome notation, consider defining the rate expression in Python vianrnpython()
. That is, we could write// here: ca + buf <> cabuf, kf = 1, kb = 0.1 nrnpython("from neuron import h") nrnpython("ca_plus_buf = h.ca + h.buf") buffering = rxd.Reaction(pyobj.ca_plus_buf, cabuf, 1, 0.1)
This is admittedly longer than the previous example, but it allows the creation of relatively complicated expressions for rate constants:
nrnpython("from neuron import h") nrnpython("kf = h.ca ** 2 / (h.ca ** 2 + (1e-3) ** 2)") // and then work with pyobj.kf
For more, see the
rxd.Reaction
entry in the Python Programmer’s reference.
- class rxd.Rate
Declare a contribution to the rate of change of a species or other state variable.
Syntax:
r = rxd.Rate(species, rate, regions, membrane_flux)
regions
andmembrane_flux
are optional, but ifmembrane_flux
is specified, thenregions
(the set of regions where the rate occurs) must also be specified. The default behavior is that the rate applies on all regions where all involved species are present; this region rule applies whenregions
is ommitted orpyobj.None
.Example:
constant_production = rxd.Rate(protein, k)
If this was the only contribution to protein dynamics and there was no diffusion, the above would be equivalent to:
dprotein/dt = k
If there are multiple
rxd.Rate
objects (or anrxd.Reaction
, etc) acting on the same species, then their effects are summed.
- class rxd.MultiCompartmentReaction
Specify a reaction spanning multiple regions to be added to the system. Use this for, for example, pumps and channels, or interactions between species living in a volume (e.g. the cytosol) and species on a membrane (e.g. the plasma membrane). For each species/state/parameter, you must specify what region you are referring to, as it could be present in multiple regions. You must also specify a membrane or a border (these are treated as synonyms) that separates the regions involved in your reaction. This is necessary because the default behavior is to scale the reaction rate by the border area, as would be expected if one of the species involved is a pump that is binding to a species in the volume. If this is not the desired behavior, pass the argument
0
for thescale_by_area
field. Pass the argument1
formembrane_flux
if the reaction produces a current across the plasma membrane that should affect the membrane potential. Unlikerxd.Reaction
objects, the base units for the rates are in terms of molecules per square micron per ms.
Mathematical functions for rate expressions
NEURON’s neuron.rxd.rxdmath
module provides a number of mathematical functions that
can be used to define reaction rates. These generally mirror the functions available
through Python’s math
module but support rxd.Species
objects.
To use any of these, first do:
objref pyobj, rxdmath { // load rxdmath nrnpython("from neuron.rxd import rxdmath") pyobj = new PythonObject() rxdmath = pyobj.rxdmath }For a full runnable example, see this Python tutorial which as here uses the hyperbolic tangent as an approximate on/off switch for the reaction.
Full documentation on this submodule (under the Python programmer’s reference, but the HOC interface is identical) is available here or you may go
directly to the documentation for any of its specific functions:
rxdmath.acos()
, rxdmath.acosh()
, rxdmath.asin()
,
rxdmath.asinh()
, rxdmath.atan()
, rxdmath.atan2()
,
rxdmath.ceil()
, rxdmath.copysign()
,
rxdmath.cos()
, rxdmath.cosh()
,
rxdmath.degrees()
, rxdmath.erf()
,
rxdmath.erfc()
, rxdmath.exp()
,
rxdmath.expm1()
, rxdmath.fabs()
,
rxdmath.factorial()
, rxdmath.floor()
,
rxdmath.fmod()
, rxdmath.gamma()
,
rxdmath.lgamma()
, rxdmath.log()
,
rxdmath.log10()
, rxdmath.log1p()
,
rxdmath.pow()
, rxdmath.pow()
,
rxdmath.sin()
, rxdmath.sinh()
,
rxdmath.sqrt()
, rxdmath.tan()
,
rxdmath.tanh()
, rxdmath.trunc()
,
rxdmath.vtrap()
Manipulating nodes
A rxd.node.Node
represents a particular state value or rxd.Parameter
in a particular location. Individual rxd.node.Node
objects are typically obtained either from being passed to an initialization function or by filtering or selecting from an rxd.nodelist.NodeList
returned by rxd.Species.nodes
. Node objects are often used for recording concentration using rxd.node.Node._ref_concentration
.
- class rxd.nodelist.NodeList
An
rxd.nodelist.NodeList
is a subclass of a Python list containingrxd.node.Node
objects. It is not intended to be created directly in a model, but rather is returned byrxd.Species.nodes
.Standard Python list methods are supported, including
.append(node)
,.extend(node_list)
,.insert(i, node)
,.index(node)
. To access the item in the ith position (0-indexed) of a NodeListnl
in HOC, usenl.__get__item(i)
. (In Python, one could saynl[i]
, but this notation is not supported by HOC.)A key added functionality is the ability to filter the nodes by rxd property (returning a new
rxd.nodelist.NodeList
). Any filter object supported by the.satifies
method of the node types present in therxd.nodelist.NodeList
may be passed in parentheses; e.g.To filter a
new_node_list
to only contain nodes present in therxd.Region
er
:just_er = new_node_list(er)
See the entry for
rxd.nodelist.NodeList
in the Python programmer’s reference for more information.
Membrane potential
- property rxd.v
A special object representing the local membrane potential in a reaction-rate expression. This can be used with
rxd.Rate
andrxd.MultiCompartmentReaction
to build ion channel models as an alternative to using NMODL, NeuroML (and converting to NMODL via jneuroml), the ChannelBuilder, orKSChan
.(If you want a numeric value for the current membrane potential at a segment
sec(x)
usesec.v(x)
instead; this syntax is slightly different from the Python convention for accessing membrane potential.)
Synchronization with segments
Changes to rxd.Species
node concentrations are propagated to segment-level concentrations automatically no later
than the next time step. This is generally the right direction for information to flow, however NEURON also provides
a rxd.re_init()
function to transfer data from segments to rxd.Species
.
- rxd.re_init()
Reinitialize all
rxd.Species
,rxd.State
, andrxd.Parameter
from changes made to NEURON segment-level concentrations. This calls the correspondingrxd.Species.re_init()
methods. Note that reaction-diffusion models may contain concentration data at a finer-resolution than that of a segment (e.g. for models being simulated in 3D).Syntax:
rxd.re_init()
Numerical options
- rxd.nthread()
Specify a number of threads to use for extracellular and 3D intracellular simulation. Currently has no effect on 1D reaction-diffusion models.
Syntax:
rxd.nthread(num_threads)
Example:
To simulate using 4 threads:
rxd.nthread(4)
Thread scaling performance is discussed in the NEURON extracellular and 3D intracellular methods papers.
- rxd.set_solve_type()
Specify numerical discretization and solver options. Currently the main use is to indicate Sections where reaction-diffusion should be simulated in 3D.
Syntax:
rxd.set_solve_type(domain, dimension)
where:
domain
– aSectionList
or other iterable of sections. Passpyobj.None
to apply the specification to the entire model.dimension
– 1 or 3
This function may be called multiple times; the last setting for dimension for a given section will apply. Different sections may be simulated in different dimensions (a so-called hybrid model).
Error handling
Most errors in the usage of the rxd
module should
raise a rxd.RxDException
.
- class rxd.RxDException
An exception originating from the
rxd
module due to invalid usage. This allows distinguishing such exceptions from other errors. HOC’s support for Error Handling is limited, so it is generally difficult to get access to these objects inside HOC, but they might be passed to HOC via a function called in Python.The text message of an
rxd.RxDException
e
may be read in HOC ase.__str__()
.