Reaction-Diffusion: Thresholds





Suppose we have an rxd.Reaction or rxd.Rate that should only occur when the concentration is above (or below) a certain threshold. These functions, however, only support continuous rate functions. What can we do?


One approach is to use a sigmoid function such as $\tanh(x)$:


Consider the following transformation of $\tanh(x)$:


$$f(x) = \frac{1 + \tanh(2m(x-a))}{2}$$


One can show that $\displaystyle \lim_{x \to \infty} f(x) = 1$, $\displaystyle \lim_{x \to -\infty} f(x) = 0$, $\displaystyle f(a) = 0.5,$ and $\displaystyle f'(a) = m$. Furthermore $f$ is a sigmoid function that shifts between $0$ and $1$ arbitrarily quickly (parameterized by $m$) around $x=a$.


Here, for example, is the graph of $\displaystyle g(x) = \frac{1 + \tanh(2\cdot 10(x-2))}{2}$:


Using this logic, we can scale reaction rates by a function of the form $f(x)$ for suitably chosen $a$ and $m$ to approximately threshold them by a concentration.


For example, suppose we wish to model a substance (we'll arbitrarily call it IP3) that degrades exponentially (i.e. $y'=-k y$) but only when the concentration is above $0.25$:

In [2]:
from neuron import h, rxd
from matplotlib import pyplot as plt

soma = h.Section(name='soma')
cyt = rxd.Region([soma], name='cyt', nrn_region='i')
ip3 = rxd.Species(cyt, name='ip3', charge=0)
k = 2 # degradation rate
threshold = 0.25 # mM... called 'a' in f(x)
m = 100 # steepness of switch
degradation_switch = (1 + rxd.rxdmath.tanh((ip3 - threshold) * 2 * m)) / 2
degradation = rxd.Rate(ip3, -k * ip3 * degradation_switch)

# prior to NEURON 7.7, this first finitialize is necessary for the pointers to exist below

t = h.Vector()
ip3_conc = h.Vector()


plt.plot(t, ip3_conc)