Updated rxd example on initial conditions and parameters (#1036)

docs/rxd-tutorials/changing intial conditions.ipynb → docs/rxd-tutorials/Changing initial conditions and parameters.ipynb

@@ -1,101 +1,135 @@
 {
  "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "id": "nVDaFx76diTC"
+   },
+   "source": [
+    "A version of this notebook may be run online via Google Colab at https://tinyurl.com/rxd-initial-and-parameters\n",
+    " (make a copy or open in playground mode)."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "id": "qE7zg9QAa7Iz"
+   },
    "outputs": [],
    "source": [
     "%matplotlib inline"
    ]
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "id": "IXXDXQDXa7I1"
+   },
    "source": [
     "Often we will want to see how the choice of initial conditions affects the dynamics. We can do this by setting the initial attribute of an `rxd.Species` and rerunning."
    ]
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "id": "10GRCniZa7I2"
+   },
    "source": [
     "For example, suppose at a single point we have the bistable dynamics introduced in the first part of this tutorial. That is, $u'=-u(1-u)(\\alpha - u)$. (Here we use $u$ instead of a specific molecule name to indicate that we are not describing any particular molecule's kinetics, but rather modeling a class of phenomena.)"
    ]
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "id": "iLp5nQJDa7I2"
+   },
    "source": [
     "This time, we'll use an `rxd.Parameter` for $\\alpha$ instead of a constant. This offers two advantages: (1) this allows $\\alpha$ to vary spatially, and (2) this allows us to change the parameter values and rerun without changing the reaction/rate specification."
    ]
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "id": "daqhJHCva7I2"
+   },
    "source": [
     "We start by defining the model and setting up the recordings:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "id": "SgcbOYX7a7I2"
+   },
    "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",
     "u = rxd.Species(cyt, name='u')\n",
-    "α = rxd.Parameter(cyt, initial=0.3)\n",
+    "α = rxd.Parameter(cyt, value=0.3 * mM)\n",
     "\n",
-    "rate = rxd.Rate(u, -u * (1 - u) * (α - u))\n",
+    "rate = rxd.Rate(u, -u * (1 * mM - u) * (α - u))\n",
     "\n",
-    "h.finitialize(-65)\n",
+    "h.finitialize(-65 * mV)\n",
     "\n",
-    "t = h.Vector()\n",
-    "t.record(h._ref_t)\n",
-    "y = h.Vector()\n",
-    "y.record(soma(0.5)._ref_ui)"
+    "t = h.Vector().record(h._ref_t)\n",
+    "y = h.Vector().record(soma(0.5)._ref_ui)"
    ]
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "id": "0wTdBqSYa7I7"
+   },
    "source": [
     "Now we define a function that sets the initial value for ca, runs the simulation, and plots the curves:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "id": "oXMtBwtYa7I7"
+   },
    "outputs": [],
    "source": [
     "def plot_it(ca_init):\n",
     "    u.initial = ca_init\n",
-    "    h.finitialize(-65)\n",
-    "    h.continuerun(10)\n",
+    "    h.finitialize(-65 * mV)\n",
+    "    h.continuerun(10 * ms)\n",
     "    plt.plot(t, y)"
    ]
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "id": "O1zJcBL1a7I7"
+   },
    "source": [
     "And now let's run it a few times and plot it:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "colab": {
+     "base_uri": "https://localhost:8080/",
+     "height": 279
+    },
+    "id": "XjRfcODga7I8",
+    "outputId": "5c4926c9-0090-4098-8fb0-bcbb4e790899"
+   },
    "outputs": [],
    "source": [
-    "for u0 in range(-4, 15, 2):\n",
+    "for u0 in range(-4 * mM, 15 * mM, 2 * mM):\n",
     "    plot_it(u0 * 0.1)\n",
     "\n",
     "plt.xlabel('t (ms)')\n",
@@ -105,42 +139,52 @@
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "id": "ogWzDHuUa7I8"
+   },
    "source": [
     "Here we see the stable fixed points at $0$ and $1$ and how all other trajectories (except for the one starting exactly at $\\alpha$, not shown) tend towards one of those values."
    ]
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "id": "g-RmMz_Fa7I8"
+   },
    "source": [
     "<hr/>"
    ]
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "id": "G-j1FRxsa7I8"
+   },
    "source": [
-    "Alternatively, we might want to know how the kinetics change if we start at a fixed value (say 0.5) and vary $\\alpha$ instead. We can do this in the same way as above by setting the `initial` property of $\\alpha$:"
+    "Alternatively, we might want to know how the kinetics change if we start at a fixed value (say 0.5) and vary $\\alpha$ instead. We can do this in the same way as above by setting the `value` property of $\\alpha$:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "id": "czCVkzqAa7I9"
+   },
    "outputs": [],
    "source": [
     "def plot_α(α_init):\n",
-    "    u.initial = 0.5\n",
-    "    α.initial = α_init\n",
-    "    h.finitialize(-65)\n",
-    "    h.continuerun(10)\n",
-    "    plt.gca().plot(t, y, label='α = %g' % α_init)"
+    "    u.initial = 0.5 * mM\n",
+    "    α.value = α_init\n",
+    "    h.finitialize(-65 * mV)\n",
+    "    h.continuerun(10 * ms)\n",
+    "    plt.gca().plot(t, y, label=f'α = {α_init} mM')"
    ]
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "id": "rDIjkWAma7I9"
+   },
    "source": [
     "Here we have added a `label` argument to the plot and explicitly specify the `gca()` axes to allow displaying a legend; see the <a href=\"https://matplotlib.org/api/_as_gen/matplotlib.pyplot.legend.html\">matplotlib.pyplot.legend</a> documentation for more.\n",
     "\n",
@@ -150,10 +194,17 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "colab": {
+     "base_uri": "https://localhost:8080/",
+     "height": 279
+    },
+    "id": "rCECS2wea7I9",
+    "outputId": "b02cabce-2ab1-4668-e0d2-e71d9c1d9175"
+   },
    "outputs": [],
    "source": [
-    "for my_α in [0, 0.2, 0.4, 0.6, 0.8, 1]:\n",
+    "for my_α in [0 * mM, 0.2 * mM, 0.4 * mM, 0.6 * mM, 0.8 * mM, 1 * mM]:\n",
     "    plot_α(my_α)\n",
     "\n",
     "plt.gca().legend()\n",
@@ -164,6 +215,10 @@
   }
  ],
  "metadata": {
+  "colab": {
+   "name": "changing intial conditions and parameters.ipynb",
+   "provenance": []
+  },
   "kernelspec": {
    "display_name": "Python 3",
    "language": "python",
@@ -179,9 +234,9 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.7"
+   "version": "3.8.3"
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 1
 }