.. _hoc_neuron_rxd: 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: .. code:: hoc 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. :class:`Vector` instead of :hoc:class:`Vector`). You might not need to use h since when working in HOC, but it does provide certain convenient functions like :func:`h.allsec`, which returns an iterable of all sections usable with ``rxd`` without having to explicitly construct a :hoc:class:`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``: .. code:: hoc 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. :ref:`Where ` the dynamics are occurring (specified using an :hoc:class:`rxd.Region` or :hoc:class:`rxd.Extracellular`) 2. :ref:`Who ` is involved (specified using an :hoc:class:`rxd.Species` or :hoc:class:`rxd.State`) 3. :ref:`What ` the reactions are (specified using :hoc:class:`rxd.Reaction`, :hoc:class:`rxd.Rate`, or :hoc:class:`rxd.MultiCompartmentReaction`) Another key class is :hoc:class:`rxd.Parameter` for defining spatially varying parameters. Integration options may be specified using :hoc:func:`rxd.set_solve_type`. .. _hoc_rxd_where: Specifying the domain --------------------- NEURON provides two main classes for defining the domain where a given set of reaction-diffusion rules applies: :hoc:class:`rxd.Region` and :hoc:class:`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 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. hoc:class:: rxd.Region Declares a conceptual intracellular region of the neuron. Syntax: .. code:: hoc 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 :hoc:class:`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 :hoc:class:`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 :class:`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 :attr:`rxd.inside`. :attr:`rxd.membrane`, :class:`rxd.DistributedBoundary`, :class:`rxd.FractionalVolume`, :class:`rxd.Shell`, :class:`rxd.FixedCrossSection`, :class:`rxd.FixedPerimeter`, and :class:`rxd.ScalableBorder`. Extracellular regions --------------------- .. hoc:class:: rxd.Extracellular Declare a extracellular region to be simulated in 3D; unlike :hoc:class:`rxd.Region`, this *always* describes the extracellular region. See the entry for :class:`rxd.Extracellular` in the Python programmer's reference for more information. .. _hoc_rxd_who: Defining proteins, ions, etc ---------------------------- Values that are distributed spatially on an :hoc:class:`rxd.Region` or :hoc:class:`rxd.Extracellular` may be defined using an :hoc:class:`rxd.Species` if they represent things that change and diffuse, an :hoc:class:`rxd.State` if they're in fixed locations but changeable (e.g. gates in an IP3R), or an :hoc:class:`rxd.Parameter` if they are just fixed values. .. hoc:class:: rxd.Species Declare an ion/protein/etc that can react and diffuse. See the entry for :class:`rxd.Species` in the Python programmer's reference for more information. .. hoc:class:: rxd.State Like an :hoc:class:`rxd.Species` but indicates the semantics of something that is not intended to diffuse. See the entry for :class:`rxd.State` in the Python programmer's reference for more information. .. hoc:class:: rxd.Parameter Declares a parameter, that can be used in place of a :hoc:class:`rxd.Species`, but unlike :hoc:class:`rxd.Species` a parameter will not change. See the entry for :class:`rxd.Parameter` in the Python programmer's reference for more information. .. _hoc_rxd_what: Defining reactions ------------------ NEURON provides three classes for specifying reaction dynamics: :hoc:class:`rxd.Reaction` for single-compartment (local) reactions; :hoc:class:`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 :hoc:class:`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 :ref:`below `. .. hoc:class:: rxd.Reaction Reaction at a point. May be mass-action or defined via custom dynamics. Syntax: .. code:: hoc 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: .. code:: hoc // 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 :hoc:func:`nrnpython`. That is, we could write .. code:: hoc // 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: .. code:: hoc 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 :class:`rxd.Reaction` entry in the Python Programmer's reference. .. hoc:class:: rxd.Rate Declare a contribution to the rate of change of a species or other state variable. Syntax: .. code:: hoc 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: .. code:: hoc 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: .. code:: dprotein/dt = k If there are multiple :hoc:class:`rxd.Rate` objects (or an :hoc:class:`rxd.Reaction`, etc) acting on the same species, then their effects are summed. .. hoc: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 :hoc:class:`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 :hoc:class:`rxd.Species` objects. To use any of these, first do: .. code:: 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 Python tutorial <../../../rxd-tutorials/thresholds.html>`_ 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 :ref:`here ` or you may go directly to the documentation for any of its specific functions: :func:`rxdmath.acos`, :func:`rxdmath.acosh`, :func:`rxdmath.asin`, :func:`rxdmath.asinh`, :func:`rxdmath.atan`, :func:`rxdmath.atan2`, :func:`rxdmath.ceil`, :func:`rxdmath.copysign`, :func:`rxdmath.cos`, :func:`rxdmath.cosh`, :func:`rxdmath.degrees`, :func:`rxdmath.erf`, :func:`rxdmath.erfc`, :func:`rxdmath.exp`, :func:`rxdmath.expm1`, :func:`rxdmath.fabs`, :func:`rxdmath.factorial`, :func:`rxdmath.floor`, :func:`rxdmath.fmod`, :func:`rxdmath.gamma`, :func:`rxdmath.lgamma`, :func:`rxdmath.log`, :func:`rxdmath.log10`, :func:`rxdmath.log1p`, :func:`rxdmath.pow`, :func:`rxdmath.pow`, :func:`rxdmath.sin`, :func:`rxdmath.sinh`, :func:`rxdmath.sqrt`, :func:`rxdmath.tan`, :func:`rxdmath.tanh`, :func:`rxdmath.trunc`, :func:`rxdmath.vtrap` Manipulating nodes ------------------ A :hoc:class:`rxd.node.Node` represents a particular state value or :class:`rxd.Parameter` in a particular location. Individual :hoc:class:`rxd.node.Node` objects are typically obtained either from being passed to an initialization function or by filtering or selecting from an :hoc:class:`rxd.nodelist.NodeList` returned by :attr:`rxd.Species.nodes`. Node objects are often used for recording concentration using :attr:`rxd.node.Node._ref_concentration`. .. hoc:class:: rxd.nodelist.NodeList An :hoc:class:`rxd.nodelist.NodeList` is a subclass of a Python `list `_ containing :hoc:class:`rxd.node.Node` objects. It is not intended to be created directly in a model, but rather is returned by :hoc:attr:`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 :hoc:class:`rxd.nodelist.NodeList`). Any filter object supported by the ``.satifies`` method of the node types present in the :hoc:class:`rxd.nodelist.NodeList` may be passed in parentheses; e.g. To filter a ``new_node_list`` to only contain nodes present in the :hoc:class:`rxd.Region` ``er``: .. code:: hoc just_er = new_node_list(er) See the entry for :class:`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 :hoc:class:`rxd.Rate` and :hoc:class:`rxd.MultiCompartmentReaction` to build ion channel models as an alternative to using NMODL, NeuroML (and converting to NMODL via `jneuroml `_), the ChannelBuilder, or :hoc:class:`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 :hoc:class:`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 :hoc:func:`rxd.re_init` function to transfer data from segments to :hoc:class:`rxd.Species`. .. hoc:function:: rxd.re_init Reinitialize all :hoc:class:`rxd.Species`, :hoc:class:`rxd.State`, and :hoc:class:`rxd.Parameter` from changes made to NEURON segment-level concentrations. This calls the corresponding :hoc:meth:`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: .. code:: hoc rxd.re_init() Numerical options ----------------- .. hoc:function:: 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: .. code:: hoc rxd.nthread(num_threads) Example: To simulate using 4 threads: .. code:: hoc rxd.nthread(4) Thread scaling performance is discussed in the NEURON `extracellular `_ and `3D intracellular `_ methods papers. .. hoc:function:: 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: .. code:: hoc rxd.set_solve_type(domain, dimension) where: - ``domain`` -- a :hoc:class:`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 :hoc:class:`rxd.RxDException`. .. hoc: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 :hoc:class:`rxd.RxDException` ``e`` may be read in HOC as ``e.__str__()``.