Add astrocyte_lr_1994 to NEST for the modeling of astrocytes and neuron-astrocyte interactions

doc/htmldoc/examples/index.rst

@@ -54,6 +54,13 @@ PyNEST examples
            * :doc:`../auto_examples/pong/run_simulations`
            * :doc:`../auto_examples/pong/generate_gif`
 
+    .. grid-item-card:: Astrocytes
+           :img-top: ../static/img/astrocyte_tripartite.png
+
+           * :doc:`../auto_examples/astrocyte_single`
+           * :doc:`../auto_examples/astrocyte_tripartite`
+
+
 .. grid:: 1 1 2 3
 
     .. grid-item-card:: Random balanced networks (Brunel)
@@ -316,6 +323,8 @@ PyNEST examples
    ../auto_examples/csa_example
    ../auto_examples/csa_spatial_example
    ../auto_examples/hpc_benchmark
+   ../auto_examples/astrocyte_single
+   ../auto_examples/astrocyte_tripartite
 
 .. toctree::
    :hidden:

doc/htmldoc/model_details/astrocyte_model_implementation.ipynb

@@ -0,0 +1,326 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "17ef5228",
+   "metadata": {},
+   "source": [
+    "# NEST implementation of the `astrocyte_lr_1994` model\n",
+    "\n",
+    "The purpose of this notebook is to provide a reference solution for the the _astrocyte_lr_1994_ model in NEST. The model is based on Li, Y. X., & Rinzel, J. (1994), De Young, G. W., & Keizer, J. (1992), and Nadkarni, S., & Jung, P. (2003).\n",
+    "\n",
+    "This notebook demonstrates how the dynamics of _astrocyte_lr_1994_ is implemented, and generates a recording of the dynamics (test_astrocyte.dat). This recording serves as a reference for the verification of _astrocyte_lr_1994_ using test_astrocyte.py in the PyNEST tests.\n",
+    "\n",
+    "## The problem\n",
+    "\n",
+    "The equations governing the evolution of the astrocyte model are\n",
+    "\n",
+    "$$\\left\\lbrace\\begin{array}{rcl}\n",
+    "    \\frac{d[\\mathrm{Ca^{2+}}](t)}{dt} &=& J_{\\mathrm{channel}}(t) - J_{\\mathrm{pump}}(t)  + J_{\\mathrm{leak}}(t) \\\\\n",
+    "    J_{\\mathrm{channel}}(t) &=& \\mathrm{r_{ER,cyt}} \\cdot \\mathrm{v_{IP_3R}} \\cdot m_{\\infty}(t)^3 \\cdot n_{\\infty}(t)^3\\cdot h_{\\mathrm{\\mathrm{IP_3R}}}(t)^3 \\cdot \\left([\\mathrm{Ca^{2+}}]_{\\mathrm{ER}} - [\\mathrm{Ca^{2+}}](t)\\right) \\\\\n",
+    "    J_{\\mathrm{pump},k}(t) &=& \\frac{\\mathrm{v_{SERCA}} \\cdot  [\\mathrm{Ca^{2+}}](t)^2}{K_{\\mathrm{m,SERCA}}^2+[\\mathrm{Ca^{2+}}](t)^2} \\\\\n",
+    "    J_{\\mathrm{leak}}(t) &=& \\mathrm{r_{ER,cyt}} \\cdot \\mathrm{v_L} \\cdot \\left( [\\mathrm{Ca^{2+}}]_{\\mathrm{ER}}(t) - [\\mathrm{Ca^{2+}}](t) \\right) \\\\\n",
+    "    m_{\\infty}(t) &=& \\frac{[\\mathrm{IP_3}](t)}{[\\mathrm{IP_3}](t) + \\mathrm{K_{d,IP_3,1}}} \\\\\n",
+    "    n_{\\infty}(t) &=& \\frac{[\\mathrm{Ca^{2+}}](t)}{[\\mathrm{Ca^{2+}}](t) + \\mathrm{K_{d,act}}} \\\\\n",
+    "    [\\mathrm{Ca^{2+}}]_{\\mathrm{ER}}(t) &=& \\frac{[\\mathrm{Ca^{2+}}]_{\\mathrm{tot}}-[\\mathrm{Ca^{2+}}](t)}{\\mathrm{r_{ER,cyt}}} \\\\\n",
+    "    \\frac{dh_{\\mathrm{IP_3}}(t)}{dt} &=& \\alpha_{h_{\\mathrm{IP_3}}}(t) \\cdot (1-h_{\\mathrm{IP_3}}(t)) - \\beta_{h_{\\mathrm{IP_3}}}(t) \\cdot h_{\\mathrm{IP_3}} (t) \\\\\n",
+    "    \\alpha_{h_{\\mathrm{IP_3}}}(t) &=& \\mathrm{k_{IP_3R}} \\cdot \\mathrm{K_{d,inh}} \\cdot  \\frac{[\\mathrm{IP_3}](t) + \\mathrm{K_{d,IP_3,1}}}{[\\mathrm{IP_3}](t)+\\mathrm{K_{d,IP_3,2}}} \\\\\n",
+    "    \\beta_{h_{\\mathrm{IP_3}}}(t) &=& \\mathrm{k_{IP_3R}} \\cdot [\\mathrm{Ca^{2+}}](t) \\\\\n",
+    "    \\frac{d[\\mathrm{IP_3}](t)}{dt} &=& \\frac{[\\mathrm{IP_3}]^{*} - [\\mathrm{IP_3}](t)}{\\tau_\\mathrm{IP_3}} + \\Delta_\\mathrm{IP3} \\cdot J_\\mathrm{syn}(t)\n",
+    "\\end{array}\\right.$$\n",
+    "\n",
+    "where $[\\mathrm{IP_3}]$, $[\\mathrm{Ca^{2+}}]$, and $h_{\\mathrm{IP_3}}$ are the state variables of interest."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "622a93cb",
+   "metadata": {},
+   "source": [
+    "## Technical details and requirements\n",
+    "\n",
+    "### Implementation of the functions\n",
+    "\n",
+    "* The reference solution is implemented through [scipy](http://www.scipy.org/).\n",
+    "\n",
+    "### Requirements\n",
+    "\n",
+    "To run this notebook, you need:\n",
+    "\n",
+    "* [numpy](http://www.numpy.org/) and [scipy](http://www.scipy.org/)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "76d05384",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from scipy.integrate import odeint"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "17cecb4d",
+   "metadata": {},
+   "source": [
+    "## Reference solution\n",
+    "### Right hand side function"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "967b18fb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def rhs(y, _, p):\n",
+    "    \"\"\"\n",
+    "    Implementation of astrocyte dynamics.\n",
+    "\n",
+    "    Parameters\n",
+    "    ----------\n",
+    "    y : list\n",
+    "        Vector containing the state variables [IP3, Ca, h_IP3R]\n",
+    "    _ : unused var\n",
+    "    p : Params instance\n",
+    "        Object containing the astrocyte parameters.\n",
+    "\n",
+    "    Returns\n",
+    "    -------\n",
+    "    dCa : double\n",
+    "        Derivative of Ca\n",
+    "    dIP3 : double\n",
+    "        Derivative of IP3\n",
+    "    dh_IP3R : double\n",
+    "        Derivative of h_IP3R\n",
+    "    \"\"\"\n",
+    "    IP3 = y[0]\n",
+    "    Ca = y[1]\n",
+    "    h_IP3R = y[2]\n",
+    "\n",
+    "    Ca = max(0, min(Ca, p.Ca_tot * (1 + p.ratio_ER_cyt)))\n",
+    "    alpha = p.k_IP3R * p.Kd_inh * (IP3 + p.Kd_IP3_1) / (IP3 + p.Kd_IP3_2)\n",
+    "    beta = p.k_IP3R * Ca\n",
+    "    Ca_ER = p.Ca_tot + (p.Ca_tot - Ca) / p.ratio_ER_cyt\n",
+    "    m_inf = IP3 / (IP3 + p.Kd_IP3_1)\n",
+    "    n_inf = Ca / (Ca + p.Kd_act)\n",
+    "    J_ch = p.ratio_ER_cyt * p.rate_IP3R * ((m_inf * n_inf * h_IP3R) ** 3) * (Ca_ER - Ca)\n",
+    "    J_pump = p.rate_SERCA * (Ca**2) / (p.Km_SERCA**2 + Ca**2)\n",
+    "    J_leak = p.ratio_ER_cyt * p.rate_L * (Ca_ER - Ca)\n",
+    "\n",
+    "    dCa = J_ch - J_pump + J_leak\n",
+    "    dIP3 = (p.IP3_0 - IP3) / p.tau_IP3\n",
+    "    dh_IP3R = alpha * (1 - h_IP3R) - beta * h_IP3R\n",
+    "\n",
+    "    return dIP3, dCa, dh_IP3R"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6f9246f9",
+   "metadata": {},
+   "source": [
+    "### Complete model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "6276334c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def scipy_astrocyte(p, f, simtime, dt, spk_ts, spk_ws):\n",
+    "    \"\"\"\n",
+    "    Complete astrocyte model using scipy `odeint` solver.\n",
+    "\n",
+    "    Parameters\n",
+    "    ----------\n",
+    "    p : Params instance\n",
+    "        Object containing the astrocyte parameters.\n",
+    "    f : function\n",
+    "        Right-hand side function\n",
+    "    simtime : double\n",
+    "        Duration of the simulation (will run between\n",
+    "        0 and tmax)\n",
+    "    dt : double\n",
+    "        Time increment.\n",
+    "    spk_ts, spk_ws : list\n",
+    "        Times and weights of spike input to the astrocyte\n",
+    "\n",
+    "    Returns\n",
+    "    -------\n",
+    "    t : list\n",
+    "        Times at which the astrocyte state was evaluated.\n",
+    "    y : list\n",
+    "        State values associated to the times in `t`\n",
+    "    fos : list\n",
+    "        List of dictionaries containing additional output\n",
+    "        information from `odeint`\n",
+    "    \"\"\"\n",
+    "    t = np.arange(0, simtime, dt)  # time axis\n",
+    "    n = len(t)\n",
+    "    y = np.zeros((n, 3))  # state variables: IP3, Ca, h_IP3R\n",
+    "    # for the state variables, assign the same initial values as in test_astrocyte.py\n",
+    "    y[0, 0] = 1.0\n",
+    "    y[0, 1] = 1.0\n",
+    "    y[0, 2] = 1.0\n",
+    "    fos = []  # full output dict from odeint()\n",
+    "    delta_ip3 = 5.0  # parameter determining the increase in IP3 induced by synaptic input\n",
+    "\n",
+    "    # update time-step by time-step\n",
+    "    for k in range(1, n):\n",
+    "        # solve ODE from t_k-1 to t_k\n",
+    "        d, fo = odeint(f, y[k - 1, :], t[k - 1 : k + 1], (p,), full_output=True)\n",
+    "        y[k, :] = d[1, :]\n",
+    "\n",
+    "        # apply synaptic inputs (spikes)\n",
+    "        if t[k] in spk_ts:\n",
+    "            y[k, 0] += delta_ip3 * spk_ws[spk_ts.index(t[k])]\n",
+    "\n",
+    "        fos.append(fo)\n",
+    "\n",
+    "    return t, y, fos"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "246ef086",
+   "metadata": {},
+   "source": [
+    "### Parameters"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "3fb81ba9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "astro_param = {\n",
+    "    \"Ca_tot\": 2.0,\n",
+    "    \"IP3_0\": 0.16,\n",
+    "    \"Kd_act\": 0.08234,\n",
+    "    \"Kd_inh\": 1.049,\n",
+    "    \"Kd_IP3_1\": 0.13,\n",
+    "    \"Kd_IP3_2\": 0.9434,\n",
+    "    \"Km_SERCA\": 0.1,\n",
+    "    \"ratio_ER_cyt\": 0.185,\n",
+    "    \"delta_IP3\": 5.0,\n",
+    "    \"k_IP3R\": 0.0002,\n",
+    "    \"rate_L\": 0.00011,\n",
+    "    \"tau_IP3\": 7142.0,\n",
+    "    \"rate_IP3R\": 0.006,\n",
+    "    \"rate_SERCA\": 0.0009,\n",
+    "}\n",
+    "\n",
+    "\n",
+    "class Params:\n",
+    "    \"\"\"\n",
+    "    Class giving access to the astrocyte parameters.\n",
+    "    \"\"\"\n",
+    "\n",
+    "    def __init__(self):\n",
+    "        self.params = astro_param\n",
+    "        self.Ca_tot = astro_param[\"Ca_tot\"]\n",
+    "        self.IP3_0 = astro_param[\"IP3_0\"]\n",
+    "        self.Kd_act = astro_param[\"Kd_act\"]\n",
+    "        self.Kd_inh = astro_param[\"Kd_inh\"]\n",
+    "        self.Kd_IP3_1 = astro_param[\"Kd_IP3_1\"]\n",
+    "        self.Kd_IP3_2 = astro_param[\"Kd_IP3_2\"]\n",
+    "        self.Km_SERCA = astro_param[\"Km_SERCA\"]\n",
+    "        self.ratio_ER_cyt = astro_param[\"ratio_ER_cyt\"]\n",
+    "        self.delta_IP3 = astro_param[\"delta_IP3\"]\n",
+    "        self.k_IP3R = astro_param[\"k_IP3R\"]\n",
+    "        self.rate_L = astro_param[\"rate_L\"]\n",
+    "        self.tau_IP3 = astro_param[\"tau_IP3\"]\n",
+    "        self.rate_IP3R = astro_param[\"rate_IP3R\"]\n",
+    "        self.rate_SERCA = astro_param[\"rate_SERCA\"]\n",
+    "\n",
+    "\n",
+    "p = Params()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "854edc71",
+   "metadata": {},
+   "source": [
+    "### Simulate the implementation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "c83f56f0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Parameters for the simulation\n",
+    "simtime = 100.0\n",
+    "resolution = 0.1\n",
+    "spike_times = [10.0]\n",
+    "spike_weights = [1.0]\n",
+    "\n",
+    "# Simulate and get the recording\n",
+    "t, y, fos = scipy_astrocyte(p, rhs, simtime, resolution, spike_times, spike_weights)\n",
+    "data = np.concatenate((np.array([t]).T, y), axis=1)\n",
+    "\n",
+    "# Save the recording (excluding the initial values)\n",
+    "np.savetxt(\n",
+    "    \"test_astrocyte.dat\",\n",
+    "    data[1:, :],\n",
+    "    header=\"\"\"\\ntest_astrocyte.dat\\n\n",
+    "This file is part of NEST.\\n\n",
+    "This .dat file contains the recordings of the state variables of an\n",
+    "astrocyte, simulated using the ODEINT solver of SciPy, with the\n",
+    "implementation detailed in\n",
+    "``doc/htmldoc/model_details/astrocyte_model_implementation.ipynb``.\n",
+    "This data is used as reference for the tests in ``test_astrocyte.py``.\\n\n",
+    "       Times                     IP3                       Ca                     h_IP3R         \"\"\",\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "df8ff655",
+   "metadata": {},
+   "source": [
+    "-----------------------------\n",
+    "### License\n",
+    "\n",
+    "This file is part of NEST. Copyright (C) 2004 The NEST Initiative\n",
+    "\n",
+    "NEST is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 2 of the License, or (at your option) any later version.\n",
+    "\n",
+    "NEST is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

models/aeif_cond_alpha.h

@@ -100,7 +100,7 @@ and
 
  \tau_w \frac{dw}{dt} = a(V-E_L) - w
 
-For implementation details see the
+For the reference implementation of this model, see
 `aeif_models_implementation <../model_details/aeif_models_implementation.ipynb>`_ notebook.
 
 See also [1]_.