docs/rxd-tutorials/thresholds.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reaction-diffusion thresholds\n",
    "<p>Suppose we have an <code>rxd.Reaction</code> or <code>rxd.Rate</code> 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?</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One approach is to use a <a href=\"https://en.wikipedia.org/wiki/Sigmoid_function\">sigmoid function</a> such as $\\tanh(x)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import numpy\n",
    "x = numpy.linspace(-5, 5)\n",
    "y = numpy.tanh(x)\n",
    "plt.grid()\n",
    "plt.plot(x, y)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('tanh(x)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the following transformation of $\\tanh(x)$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$f(x) = \\frac{1 + \\tanh(2m(x-a))}{2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can show that\n",
    "$\\displaystyle \\lim_{x \\to \\infty} f(x) = 1$,\n",
    "$\\displaystyle \\lim_{x \\to -\\infty} f(x) = 0$,\n",
    "$\\displaystyle f(a) = 0.5,$ and \n",
    "$\\displaystyle f'(a) = m$. Furthermore $f$ is a sigmoid function that shifts between $0$ and $1$ arbitrarily quickly (parameterized by $m$) around $x=a$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, for example, is the graph of $\\displaystyle g(x) = \\frac{1 + \\tanh(2\\cdot 10(x-2))}{2}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = numpy.linspace(0, 4, 1000)\n",
    "y = (1+numpy.tanh(2*10*(x-2)))/2\n",
    "plt.grid()\n",
    "plt.plot(x, y)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('g(x)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from neuron import h, rxd\n",
    "from neuron.units import mV, ms, mM\n",
    "from matplotlib import pyplot as plt\n",
    "h.load_file('stdrun.hoc')\n",
    "\n",
    "soma = h.Section(name='soma')\n",
    "cyt = rxd.Region([soma], name='cyt', nrn_region='i')\n",
    "ip3 = rxd.Species(cyt, name='ip3', charge=0, initial=1 * mM)\n",
    "k = 2            # degradation rate\n",
    "threshold = 0.25 # mM... called 'a' in f(x)\n",
    "m = 100          # steepness of switch\n",
    "degradation_switch = (1 + rxd.rxdmath.tanh((ip3 - threshold) * 2 * m)) / 2\n",
    "degradation = rxd.Rate(ip3, -k * ip3 * degradation_switch)\n",
    "\n",
    "t = h.Vector().record(h._ref_t)\n",
    "ip3_conc = h.Vector().record(soma(0.5)._ref_ip3i)\n",
    "\n",
    "h.finitialize(-65 * mV)\n",
    "h.continuerun(2 * ms)\n",
    "\n",
    "plt.plot(t, ip3_conc)\n",
    "plt.xlabel('t (ms)')\n",
    "plt.ylabel('[IP3] (mM)')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}