Restricting reactions to a part of a region

A version of this notebook may be run online via Google Colab at https://tinyurl.com/neuron-restricted-reaction (make a copy or open in playground mode).

Overview

Reactions may be restricted to only occur in a given part of an rxd.Region by multiplying the rate by an indicator variable, that is an rxd.Parameter that takes the value 1 in the subregion where the reaction is to occur and 0 elsewhere.

Let’s consider an example.

Setup NEURON library and imports

Let’s import our usual NEURON libraries and definitions. Remember you can use either um or µm for micron.

from neuron import h, rxd
from neuron.units import mV, ms, um, mM

## needed for standard run system
h.load_file("stdrun.hoc")

1.0

Now import plotly, a graphics library. (You could easily modify this code to use other graphics libraries like matplotlib, plotnine, or bokeh.)

import plotly.graph_objects as go

Setup the model

Morphology and discretization

left = h.Section(name="left")
right = h.Section(name="right")

left.nseg = right.nseg = 101
left.L = right.L = 101

right.connect(left)

right

The rxd.Region

Our region will be present on both the left and the right sections. Since they define our entire “cell”, we will pick one and put the rxd.Region on the .wholetree() that is connected to it.

cytosol = rxd.Region(left.wholetree())

The rxd.Species

Let’s begin by defining our initialization rule. The concentration will initially be 0 everywhere except for nodes between 90 and 110 microns from the left edge of left (that is, the last 10 microns of left and the first 10 microns of right).

def initial_protein(node):
    if 90 * um < h.distance(node, left(0)) < 110 * um:
        return 1 * mM
    else:
        return 0

Now the species definition itself.

protein = rxd.Species(cytosol, d=1, initial=initial_protein)

The rxd.Reaction and its localization

We need an rxd.Parameter indicator variable to indicate the sections where the reaction is to take place. In parallel to how we defined the initial concentration of the species above, it is often simplest to define the values using a function of the node (here 1 in the right section; 0 otherwise):

def active_region_value(node):
    if node in right:
        return 1
    else:
        return 0

Now the indicator variable rxd.Parameter itself:

in_region = rxd.Parameter(cytosol, value=active_region_value)

For the reaction (here an rxd.Rate but could be an rxd.Reaction or rxd.MultiCompartmentReaction as well), note that we multiply by the value of in_region, so when that value is 1 the reaction happens and when it is 0, the reaction does not.

production_rate = 0.002 * mM / ms

reaction = rxd.Rate(protein, production_rate * in_region)

Run and visualize the simulation

Always initialize your simulations.

h.finitialize(-65 * mV)

1.0

We’ll plot the protein concentrations against their distance from the left edge of the left section (i.e. from left(0)).

def plot(fig):
    y = protein.nodes.concentration
    x = [h.distance(node, left(0)) for node in protein.nodes]
    fig.add_trace(go.Scatter(x=x, y=y, name=f"t = {h.t:g} ms"))

Plot the initial conditions, then advance for 25 ms, plot, advance again, etc.

fig = go.Figure()
plot(fig)

for advance_count in range(1, 5):
    h.continuerun(advance_count * 25 * ms)
    plot(fig)

fig.update_layout(xaxis_title="Position (µm)", yaxis_title="Concentration (mM)")

fig.show(renderer="notebook_connected")