Restricting reactions to certain sections¶
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()