{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Customizing channels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sinaps as sn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To implement a channel C, it is necessary to implement:\n",
    "\n",
    " * a static method  C._I(V,t,**st_vars,**params) returning the net  current towards inside\n",
    "  ([pA] for point channel and [pA/μm2] for density channels) of the mechanism with\n",
    " V the voltage (:array) and **st_vars the state variable of the mechanism (:array)\n",
    " * If there are state variables (ions), a static method C._dS(V,t,**st_vars,**params) returning a\n",
    "  tuple with the differential of each state variable and a function S0 returning a tuple with the initial value for each state variable\n",
    " * a static method C._J(ion,V,t,**st_vars,**params) returning the flux of each species due to the current, in mM/ms for point channels, and in mM/ms/μm^{2} ofr density channels.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example we will create one additionals channel, a NMDA receptor, with 15% of the current carried by calcium ions.\n",
    "We need to:\n",
    "\n",
    "* subclass the Channel Class\n",
    "* define the _init function\n",
    "* define the static method C._I, that computes the current, in pA for point channels, and in pA/μm^2 for density channels\n",
    "* define the static method C._J., that computes the corresponding the flux of ions, in mM/ms. The conversion from pA to mM/ms requires to divide by 96.48533132838746, and by the valence of the ion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class NMDAR(sn.Channel):\n",
    "    \"\"\"Point channel with a NMDAr-type current starting at time t0, \n",
    "voltage-dependent flow of sodium (Na+) and small amounts of calcium (Ca2+) ions into the cell and potassium (K+) out of the cell.\n",
    "    \"\"\"\n",
    "    param_names = ('t0','gnmda','tnmda1','tnmda2','V_nmda')\n",
    "\n",
    "    def __init__(self,t0,gnmda=0.02,tnmda1=11.5,tnmda2=0.67,V_nmda=75):\n",
    "        \"\"\"Point channel with a NMDAr-type current starting at time t0 [pA]\n",
    "            t0: start of the current [ms]\n",
    "            gnmda: max conductance of NMDAr [nS]\n",
    "            tnmda1: NMDAr time constant [ms]\n",
    "            tnmda2: NMDAr time constant [ms]\n",
    "            V_nmda: NMDAr Nernst potential [mV]\n",
    "        \"\"\"\n",
    "        self.params={'t0' : t0,\n",
    "                     'gnmda' : gnmda,\n",
    "                     'tnmda1' : tnmda1,\n",
    "                     'tnmda2' : tnmda2,\n",
    "                     'V_nmda': V_nmda,\n",
    "                    }\n",
    "\n",
    "    @staticmethod\n",
    "    def _I(V,t,\n",
    "           t0,gnmda,tnmda1,tnmda2,V_nmda):\n",
    "        return -((t <= t0+50) & (t >= t0))*gnmda*(np.exp(-np.abs(t-t0)/tnmda1)-np.exp(-np.abs(t-t0)/tnmda2))/(1+0.33*2*np.exp(-0.06*(V-65)))*(V-V_nmda)\n",
    "\n",
    "    \n",
    "    @staticmethod\n",
    "    def _J(ion,V,t,\n",
    "           t0,gnmda,tnmda1,tnmda2,V_nmda):\n",
    "        \"\"\"\n",
    "        Return the flux of ion [aM/ms/um2] of the mechanism towards inside.\n",
    "        \"\"\"\n",
    "        if ion is sn.Species.Ca:\n",
    "            return ((t <= t0+50) & (t >= t0)) *np.maximum(-gnmda*(np.exp(-np.abs(t-t0)/tnmda1)-np.exp(-np.abs(t-t0)/tnmda2))/(1+0.33*2*np.exp(-0.06*(V-65)))*(V-V_nmda)/96.48533132838746/2*0.15,0)\n",
    "        else:\n",
    "            return 0 *V"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting up the channels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nrn = sn.Neuron([(0,1),(1,2),(1,3),(3,4)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nrn[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nrn[:].clear_channels()\n",
    "nrn[:].add_channel(sn.channels.LeakChannel())\n",
    "nrn[:].add_channel(sn.channels.Hodgkin_Huxley_Ca())\n",
    "\n",
    "nrn[0].add_channel(NMDAR(0.5,gnmda=20),0) #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Shortcuts\n",
    "Ca = sn.Species.Ca"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nrn.add_species(Ca, D=50, C0=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running the simu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialisation of the simulation\n",
    "sim=sn.Simulation(nrn,dx=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Runing the simulation\n",
    "sim.run((0,100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting the potential\n",
    "sim.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.plot.I(NMDAR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running the calcium dynamics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Running the chemical reactions part   \n",
    "sim.run_diff(max_step=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim[:].plot.C(Ca)"
   ]
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}