thresholds

Overview

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

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

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

In [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()