{
 "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_astro, h_IP3R]\n",
    "    _ : unused var\n",
    "    p : Params instance\n",
    "        Object containing the astrocyte parameters.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    dCa_astro : double\n",
    "        Derivative of Ca_astro\n",
    "    dIP3 : double\n",
    "        Derivative of IP3\n",
    "    dh_IP3R : double\n",
    "        Derivative of h_IP3R\n",
    "    \"\"\"\n",
    "    IP3 = y[0]\n",
    "    Ca_astro = y[1]\n",
    "    h_IP3R = y[2]\n",
    "\n",
    "    Ca_astro = max(0, min(Ca_astro, 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_astro\n",
    "    Ca_ER = (p.Ca_tot - Ca_astro) / p.ratio_ER_cyt\n",
    "    m_inf = IP3 / (IP3 + p.Kd_IP3_1)\n",
    "    n_inf = Ca_astro / (Ca_astro + p.Kd_act)\n",
    "    J_ch = p.ratio_ER_cyt * p.rate_IP3R * ((m_inf * n_inf * h_IP3R) ** 3) * (Ca_ER - Ca_astro)\n",
    "    J_pump = p.rate_SERCA * (Ca_astro**2) / (p.Km_SERCA**2 + Ca_astro**2)\n",
    "    J_leak = p.ratio_ER_cyt * p.rate_L * (Ca_ER - Ca_astro)\n",
    "\n",
    "    dCa_astro = 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_astro, 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_astro, 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 = p.delta_IP3  # 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\": 0.0002,\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_astro                  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.12.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}