rxd.nthread · rxd.re_init · rxd.set_solve_type

rxd.Extracellular

rxd.MultiCompartmentReaction

rxd.nodelist.NodeList

rxd.Parameter

rxd.Rate

rxd.Reaction

rxd.Region

rxd.RxDException

rxd.Species

rxd.State

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:

  1. Where the dynamics are occurring (specified using an rxd.Region or rxd.Extracellular)

  2. Who is involved (specified using an rxd.Species or rxd.State)

  3. 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.

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 a SectionList 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 or pyobj.None for none of the above.

  • name is the name of the region (e.g. cyt or er); this has no effect on the simulation results but it is helpful for debugging

  • dx deprecated; when specifying name pass in pyobj.None here

  • dimension deprecated; when specifying name pass in pyobj.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 membrane

  • pyobj.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 to rxd.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 unlike rxd.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, 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 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 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

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.

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 containing rxd.node.Node objects. It is not intended to be created directly in a model, but rather is returned by rxd.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 NodeList nl in HOC, use nl.__get__item(i). (In Python, one could say nl[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 the rxd.nodelist.NodeList may be passed in parentheses; e.g.

To filter a new_node_list to only contain nodes present in the rxd.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 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 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()

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 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 – 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).

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