# 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?

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

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

[1]:

from matplotlib import pyplot as plt
import numpy

x = numpy.linspace(-5, 5)
y = numpy.tanh(x)
plt.grid()
plt.plot(x, y)
plt.xlabel("x")
plt.ylabel("tanh(x)")
plt.show()


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}$$:

[2]:

x = numpy.linspace(0, 4, 1000)
y = (1 + numpy.tanh(2 * 10 * (x - 2))) / 2
plt.grid()
plt.plot(x, y)
plt.xlabel("x")
plt.ylabel("g(x)")
plt.show()


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$$:

[3]:

from neuron import h, rxd
from neuron.units import mV, ms, mM
from matplotlib import pyplot as plt

h.load_file("stdrun.hoc")

soma = h.Section(name="soma")
cyt = rxd.Region([soma], name="cyt", nrn_region="i")
ip3 = rxd.Species(cyt, name="ip3", charge=0, initial=1 * mM)
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)

t = h.Vector().record(h._ref_t)
ip3_conc = h.Vector().record(soma(0.5)._ref_ip3i)

h.finitialize(-65 * mV)
h.continuerun(2 * ms)

plt.plot(t, ip3_conc)
plt.xlabel("t (ms)")
plt.ylabel("[IP3] (mM)")
plt.show()