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 mayalso be specified using HOC.

We can access the rxd module via:

Python

from neuron import rxd

HOC

objref pyobj, n, rxd

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

The above additionally provides access to an object called n which is traditionally how Python accesses core NEURON functionality (e.g. in Python one would use n. Vector instead of Vector). You might not need to use n since when working in HOC, but it does provide certain convenient functions like n.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, n, rxd, cyt, ca, buf, cabuf, buffering, g

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

{
    // define the domain and the dynamics
    create soma

    cyt = rxd.Region(n.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).

Note: In older Python code, you may find from neuron import crxd as rxd but this is equivalent to the above as crxd has been an alias for rxd for several years.

In general, a reaction-diffusion model specification involves answering three conceptual questions:

Where the dynamics are occurring (specified using an rxd.Region or rxd.Extracellular)
Who is involved (specified using an rxd.Species or rxd.State)
What the reactions are (specified using rxd.Reaction, rxd.Rate, or rxd.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.

Extracellular regions

class rxd.Extracellular

Python

Declare a extracellular region to be simulated in 3D; unlike rxd.Region, this always describes the extracellular region.

Syntax:

r = rxd.Extracellular(xlo, ylo, zlo, xhi, yhi, zhi, dx,
                    volume_fraction=1, tortuosity=None, permeability=None)
Here:

xlo, …, zhi – define the bounding box of the domain; it should typically contain the entire morphology… by default NEURON assumes reflective (Neumann) boundary conditions, so you may want the box to extend well beyond the cell morphology depending on your use case

dx – voxel edge size in µm

tortuosity – increase factor in path length due to obstacles, effective diffusion coefficient d/tortuosity^2; either a single value for the whole region or a Vector giving a value for each voxel. Default is 1 (no change).

volume_fraction – the free fraction of extracellular space; a volume_fraction of 1 assumes no cells; lower values are probably warranted for most simulations

Example:

A tutorial demonstrating extracellular diffusion is available here.

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

Python

Declare an ion/protein/etc that can react and diffuse.

Syntax:

s = rxd.Species(regions=None,
                d=0,
                name=None,
                charge=0,
                initial=None,
                atolscale=1,
                ecs_boundary_conditions=None,
                represents=None
)
Parameters:

regions – a Region or list of Region objects containing the species

d – the diffusion constant of the species (optional; default is 0, i.e. non-diffusing)

name – the name of the Species; used for syncing with NMODL and HOC (optional; default is none)

charge – the charge of the Species (optional; default is 0)

initial – the initial concentration or None (if None, then imports from HOC if the species is defined at finitialize, else 0); can be specified as a constant or as a function of an rxd.Node

atolscale – scale factor for absolute tolerance in variable step integrations

ecs_boundary_conditions – if Extracellular rxd is used ecs_boundary_conditions=None for zero flux boundaries or if ecs_boundary_conditions=the concentration at the boundary.

represents – optionally provide CURIE (Compact URI) to annotate what the species represents e.g. CHEBI:29101 for sodium(1+)

Note

Charge must match the charges specified in NMODL files for the same ion, if any. Common species charges include: sodium (+1), potassium (+1), calcium (+2), magnesium (+2), chloride (-1).

You probably want to adjust atolscale for species present at low concentrations (e.g. calcium).

NEURON does not require any specific ontology for identifiers, however CHEBI contains identifiers for many substances of interest in reaction-diffusion modeling. A number of ontology search providers are available on the internet, including BioPortal.

To refer to a given Species restricted to a certain region, specify the region in square brackets. e.g., er_calcium = ca[er].

property nodes

An rxd.nodelist.NodeList of all the nodes corresponding to the species.

This can then be further restricted using the callable property of NodeList objects.

Example: All nodes from sp on the Section dend:

nodes_on_dend = sp.nodes(dend)

Example: All nodes from sp on the Segment dend(0.5):

nodes_on_dend_center = sp.nodes(dend(0.5))
For 1D simulation with a species defined on 1 region, this will be a rxd.nodelist.NodeList of length 1 (or 0 if the species is not defined on the segment); for 3D simulation or with multiple regions, the list may be longer and further filtering may be required.

Note

The code

node_list = ca[cyt].nodes

is more efficient than the otherwise equivalent

node_list = ca.nodes(cyt)

because the former only creates the rxd.Node objects belonging to the restriction ca[cyt] whereas the second option constructs all rxd.Node objects belonging to the rxd.Species ca and then culls the list to only include those also belonging to the rxd.Region cyt.

property states: A list of a copy of the state values corresponding to this species; modifying the values in the list will not change the values used in the simulation.

re_init()

Syntax:

sp.re_init()
where sp is an instance of an rxd.Species.

Reinitialize the rxd concentration of this species to match the NEURON grid. Used when e.g. cai is modified directly instead of through a corresponding rxd.Species.

concentrations(): Deprecated. Do not use.

property charge

Get or set the charge of the Species.

Setting was added in NEURON 7.4 and is allowed only before the reaction-diffusion model is instantiated.

property regions

Get or set the regions where the Species is present.

Setting was added in NEURON 7.4 and is allowed only before the reaction-diffusion model is instantiated. Getting was added in NEURON 7.5.

indices()

Return the indices corresponding to this species in the given region.

This is occasionally useful, but it’s generally best to avoid explicit indexing when developing models.

Syntax:

sp.indices(r=None, secs=None)

where sp is an instance of an rxd.Species.

If r is None, then returns all species indices. If r is a list of regions return indices for only those sections that are on all the regions. If secs is a set return all indices on the regions for those sections.

defined_on_region()

Returns True if the rxd.Species sp is present on r, else False.

Syntax:

result = sp.defined_on_region(r)

property name

Get or set the name of the Species.

This is used only for syncing with the non-reaction-diffusion parts of NEURON (NMODL, HOC).

Setting requires NEURON 7.4+, and then only before the reaction-diffusion model is instantiated.

class rxd.State: Python

Like an rxd.Species but indicates the semantics of something that is not intended to diffuse.

For example, one might use a State to represent the gating variables on an IP3 receptor on the endoplasmic reticulum (ER).

class rxd.Parameter

python

Declares a parameter, that can be used in place of a rxd.Species, but unlike rxd.Species a parameter will not change.

Syntax:

s = rxd.Parameter(regions, name=None, charge=0,
                value=None, represents=None)

Parameters:

regions – a rxd.Region or list of rxd.Region objects containing the parameter

name – the name of the parameter; used for syncing with NMODL and HOC (optional; default is none)

charge – the charge of the Parameter (optional, probably only rarely needed; default is 0)

value – the value or None (if None, then imports from HOC if the parameter is defined at finitialize(), else 0)

represents – optionally provide CURIE (Compact URI) to annotate what the parameter represents e.g. CHEBI:29101 for sodium(1+)

Note

Charge must match the charges specified in NMODL files for the same ion, if any.

Attempting to specify a non-zero diffusion rate for an rxd.Parameter will raise an rxd.RxDException.

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

Python

Syntax:

r1 = rxd.Reaction(reactant_sum, product_sum, forward_rate,
                backward_rate=0, regions=region_list, custom_dynamics=False)
Here:

reactant_sum is a combination of rxd.Species, rxd.State,

rxd.Parameter, e.g. ca + 2 * cl representing the left-hand-side of the reaction. * product_sum is like reactant_sum but represdenting the right-hand-side. * forward_rate and backward_rate represent the reaction rates; reactions are assumed to be governed by mass-action kinetics with these as the rate constants unless custom_dynamics is true, in which case these are expressions fully defining the rate of change. In particular, these can be constants or expressions combining rxd.Species etc with constants, arithmetic, and rxd.v * region_list is a list of regions on which this reaction occurs. If ommitted or None, the reaction occurs on all regions where all involved species are defined.

Examples:

For the following, we assume that H, O, H2O, dimer, and decomposed are instances of rxd.Species.

For 2 * H + O > H2O in a mass action reaction at rate k:
r = rxd.Reaction(2 * H + O, H2O, k)
To constrain the reaction to a specified list of regions, say to just the extracellular space (ext) and the cytosol (cyt), use the regions keyword, e.g.
r = rxd.Reaction(2 * H + O, H2O, k, regions=[ext, cyt])
For a bi-directional reaction, specify a backward reaction rate. e.g. if kf is the forward rate and kb is the backward rate, then:
r = rxd.Reaction(2 * H + O, H2O, kf, kb)
To use dynamics other than mass-action, add that mass_action=False flag and put the full formula instead of a mass-action rate for kf (and kb). E.g. Michaelis-Menten degradation
r = rxd.Reaction(
    dimer, decomposed, dimer / (k + diamer), mass_action=False
)

property f_rate: Get or set the forward reaction rate.

property b_rate: Get or set the backward reaction rate.

HOC

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, and custom_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 a backward_rate or restricting to specific regions, pass 0 for backward_rate and pyobj.None for regions.

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 via nrnpython(). That is, we could write

// here: ca + buf <> cabuf, kf = 1, kb = 0.1
nrnpython("from neuron import n")
nrnpython("ca_plus_buf = n.ca + n.buf")
buffering = rxd.Reaction(pyobj.ca_plus_buf, cabuf, 1, 0.1)

This is admittedly longer than the previous example, but it simplifies the creation of and increases the readability of relatively complicated expressions for rate constants:

nrnpython("from neuron import n")
nrnpython("kf = n.ca ** 2 / (n.ca ** 2 + (1e-3) ** 2)")
// and then work with pyobj.kf

property f_rate: Get or set the forward reaction rate.

property b_rate: Get or set the backward reaction rate.

class rxd.Rate

Python

Declare a contribution to the rate of change of a species or other state variable.

Syntax:

r = rxd.Rate(species, rate, regions=None, membrane_flux=False)
If regions is None, then the rate applies on all regions.

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 an rxd.Reaction, etc) acting on the same species, then their effects are summed.

HOC

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 and membrane_flux are optional, but if membrane_flux is specified, then regions (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 when regions is ommitted or pyobj.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 an rxd.Reaction, etc) acting on the same species, then their effects are summed.

class rxd.MultiCompartmentReaction: Python

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 keyword argument scale_by_area=False. Pass in membrane_flux=True if the reaction produces a current across the plasma membrane that should affect the membrane potential. Unlike rxd.Reaction objects, the base units for the rates are in terms of molecules per square micron per ms.

HOC

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 the scale_by_area field. Pass the argument 1 for membrane_flux if the reaction produces a current across the plasma membrane that should affect the membrane potential. Unlike rxd.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.

Python

To use any of these, first do:

from neuron.rxd import rxdmath

Example:

degradation_switch = (1 + rxdmath.tanh((ip3 - threshold) * 2 * m)) / 2
degradation = rxd.Rate(ip3, -k * ip3 * degradation_switch)

HOC

objref pyobj, rxdmath

{
    // load rxdmath
    nrnpython("from neuron.rxd import rxdmath")
    pyobj = new PythonObject()
    rxdmath = pyobj.rxdmath
}

For a full runnable example, see this tutorial which as here uses the hyperbolic tangent as an approximate on/off switch for the reaction.

Full documentation on this submodule 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()

Membrane potential

property rxd.v

Python

A special object representing the local membrane potential in a reaction-rate expression. This can be used with rxd.Rate and rxd.MultiCompartmentReaction to build ion channel models as an alternative to using NMODL, NeuroML (and converting to NMODL via jneuroml), the ChannelBuilder, or KSChan.

(If you want a numeric value for the current membrane potential at a segment seg use seg.v instead.)

Example (adapted from the Hodgkin Huxley via rxd tutorial)

from neuron.rxd.rxdmath import vtrap, exp, log
from neuron import rxd

alpha = 0.01 * vtrap(-(rxd.v + 55.0), 10.0)
beta = 0.125 * exp(-(rxd.v + 65.0)/80.0)
ntau = 1.0/(alpha + beta)
ninf = alpha/(alpha + beta)

# ... define cyt, mem, sections ...

ngate = rxd.State([cyt, mem], name='ngate', initial=0.24458654944007166)
n_gate = rxd.Rate(ngate, (ninf - ngate)/ntau)

HOC

(If you want a numeric value for the current membrane potential at a segment sec(x) use sec.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()

Python

Reinitialize all rxd.Species, rxd.State, and rxd.Parameter from changes made to NEURON segment-level concentrations. This calls the corresponding rxd.Species.re_init() methods. Note that reaction-diffusion models may contain concentration data at a finer-resolution than that of a nrn.Segment (e.g., for models being simulated in 3D).

Syntax:

rxd.re_init()

HOC

Syntax:

rxd.re_init()

Numerical options

rxd.nthread()

Python

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.

HOC

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

Python

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=None, dimension=None, dx=None, nsubseg=None, method=None)
where:

domain – a Section or Python iterable of sections. If the domain is None or omitted, the specification will apply to the entire model.

dimension – 1 or 3

dx – not implemented; specify dx for 3D models when creating the rxd.Region

nsubseg – not implemented

method – not implemented

This function may be called multiple times; the last setting for any given field will be used. Different sections may be simulated in different dimensions (a so-called hybrid model).

Warning

For 3D reaction-diffusion simulations, we recommend upgrading to at least NEURON 8.1.

(Calculation of 3D concentration changes from MOD file activity in NEURON 7.8.x and 8.0.x was underestimated due to an inconsistency in surface voxel partial volume calculations.)

HOC

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 – a SectionList or other iterable of sections. Pass pyobj.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).

SBML Export

rxd.export.sbml()

Python

Export dynamics at a segment to an SBML file.

Syntax:

rxd.export.sbml(segment, filename=None, model_name=None, pretty=True)

This does not currently support rxd.MultiCompartmentReaction; attempting to export dynamics that involve such reactions will raise an rxd.RxDException.

Note

Prior to NEURON 9, export functionality required also calling import neuron.rxd.export. This is no longer necessary.

Saving and restoring state

Some simulations require a lengthy initialization before exploring various possible stimuli. In these situations, it is often convenient to run the initialization once, save the state, do an experiment, and revert back to the saved state.

Beginning in NEURON 8.1, reaction-diffusion states are included automatically when using SaveState which additionally saves many other model states.

If one wants to save and restore only reaction-diffusion states, this can be done via the following functions:

The reaction-diffusion state data returned by rxd.save_state() and expected by rxd.restore_state() consists of 16 bytes of metadata (8 bytes for a version identifier and 8 bytes for the length of the remaining portion) followed by gzip-compresssed state values. In particular, not every binary string of a given length is a valid state vector, nor is every state vector for a given model necessarily the same length (as the compressability may be different).

Error handling

Most errors in the usage of the neuron.rxd module should raise a rxd.RxDException.

class rxd.RxDException: Python

An exception originating from the neuron.rxd module due to invalid usage. This allows distinguishing such exceptions from other errors.

The text message of an rxd.RxDException e may be read as str(e).

HOC

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 as e.__str__().

Basic Reaction-Diffusion

Overview

Specifying the domain

Intracellular regions and regions in Frankenhauser-Hodgkin space

For specifying the geometry

Extracellular regions

Defining proteins, ions, etc

Defining reactions

Mathematical functions for rate expressions

Manipulating nodes

Membrane potential

Synchronization with segments

Numerical options

SBML Export

Saving and restoring state

Error handling

Basic Reaction-Diffusion

Overview

Related resources

Specifying the domain

Intracellular regions and regions in Frankenhauser-Hodgkin space

For specifying the geometry

Extracellular regions

Defining proteins, ions, etc

Defining reactions

Mathematical functions for rate expressions

Manipulating nodes

Membrane potential

Synchronization with segments

Numerical options

SBML Export

Saving and restoring state

Error handling