{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# IAF neurons singularity\n", "\n", "This notebook describes how NEST handles the singularities appearing in the ODE's of integrate-and-fire model neurons with alpha- or exponentially-shaped current, when the membrane and the synaptic time-constants are identical.\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import sympy as sp\n", "sp.init_printing(use_latex=True)\n", "from sympy.matrices import zeros\n", "tau_m, tau_s, C, h = sp.symbols('tau_m, tau_s, C, h')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For alpha-shaped currents we have:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "A = sp.Matrix([[-1/tau_s,0,0],[1,-1/tau_s,0],[0,1/C,-1/tau_m]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Non-singular case ($\\tau_m\\neq \\tau_s$) \n", "The propagator is: " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "PA = sp.simplify(sp.exp(A*h))\n", "PA" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the entry in the third line and the second column $A_{32}$ would also appear in the propagator matrix in case of an exponentially shaped current" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Singular case ($\\tau_m = \\tau_s$) \n", "We have" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$\\displaystyle \\left[\\begin{matrix}- \\frac{1}{\\tau_{m}} & 0 & 0\\\\1 & - \\frac{1}{\\tau_{m}} & 0\\\\0 & \\frac{1}{C} & - \\frac{1}{\\tau_{m}}\\end{matrix}\\right]$$" ], "text/plain": [ "⎡-1 ⎤\n", "⎢─── 0 0 ⎥\n", "⎢ τₘ ⎥\n", "⎢ ⎥\n", "⎢ -1 ⎥\n", "⎢ 1 ─── 0 ⎥\n", "⎢ τₘ ⎥\n", "⎢ ⎥\n", "⎢ 1 -1 ⎥\n", "⎢ 0 ─ ───⎥\n", "⎣ C τₘ⎦" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "As = sp.Matrix([[-1/tau_m,0,0],[1,-1/tau_m,0],[0,1/C,-1/tau_m]])\n", "As" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The propagator is" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$\\displaystyle \\left[\\begin{matrix}e^{- \\frac{h}{\\tau_{m}}} & 0 & 0\\\\h e^{- \\frac{h}{\\tau_{m}}} & e^{- \\frac{h}{\\tau_{m}}} & 0\\\\\\frac{h^{2} e^{- \\frac{h}{\\tau_{m}}}}{2 C} & \\frac{h e^{- \\frac{h}{\\tau_{m}}}}{C} & e^{- \\frac{h}{\\tau_{m}}}\\end{matrix}\\right]$$" ], "text/plain": [ "⎡ -h ⎤\n", "⎢ ─── ⎥\n", "⎢ τₘ ⎥\n", "⎢ ℯ 0 0 ⎥\n", "⎢ ⎥\n", "⎢ -h -h ⎥\n", "⎢ ─── ─── ⎥\n", "⎢ τₘ τₘ ⎥\n", "⎢h⋅ℯ ℯ 0 ⎥\n", "⎢ ⎥\n", "⎢ -h -h ⎥\n", "⎢ ─── ─── -h ⎥\n", "⎢ 2 τₘ τₘ ───⎥\n", "⎢h ⋅ℯ h⋅ℯ τₘ⎥\n", "⎢─────── ────── ℯ ⎥\n", "⎣ 2⋅C C ⎦" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "PAs = sp.simplify(sp.exp(As*h))\n", "PAs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Numeric stability of propagator elements\n", "For the lines $\\tau_s\\rightarrow\\tau_m$ the entry $PA_{32}$ becomes numerically unstable, since denominator and enumerator go to zero.\n", "\n", "**1.** We show that $PAs_{32}$ is the limit of $PA_{32}(\\tau_s)$ for $\\tau_s\\rightarrow\\tau_m$.:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAF0AAAAaCAYAAADVLFAXAAAExElEQVR4nO3ZWYhXVRwH8I+aMpnSQ7ZIliWVWVPOKDUWtpOUJePY3kNJ0UNBItFC9ZA9RAtGKhUhEZlRUblltphFDO2loi1SVkYvltqDtmgu2cPvXP63O/c/i45OA/OFy+9/zv2dc37nd3/bOf9e06ZN04P9i95dLUA78DDe7GohOhPdQel1WNUOvqfx2L4VpXPQHZQ+UttK741LsWjfi7P3KFP6bdiNa/ezLGU4AodjB97Cn/gB5xf4GtAXH6T2QrGHsue6fS10WyhT+qhEV+xPQaqgPtEpmC6s/ms8WuBrxBLsTO0bMRjH4B9ckdqD8eJeyHM4dmFWybvj8SCWY6MwlI1YhlvQP2M8oGTwKGFR3+2FcJ2FOmzBVVif+ubhgQJfI+7NtX9LtF4YVjM2dII8jWm+Bbm+Xrgfd6EfPsKr2IyhGIcLcA3OoqXSD8JwfCwspLMxDfe1wXMe3k+/6/CaisJhGL7PtU8QFv12yVyn4hedo3BoEh+0Odf3DCbjGxGSi/nnQBGyj8s6iuGlLvWtwEl4Pgn9h/iCDVWEuUyUdZuwHWtxD/oU+B7HiDaezwryfFyYox4rc+1GvCu8s4hTlSfhxZiBT5KsDcJ6f8Lt5Vt0sMgli0WIIZQ5GWswpspaW4Vn3pJ1FC19dKJD8UXazHM4GeOF1R2H3xNfH/FhrhbW9wr+xsVpoeG4Pjf/pvS0B/3TWisL/fXCfTM04tkqcwzDtyX9tUJ5UzFb5IjxGCRyw/SSMZeI8DE/tQeLPe7ElSo6qYat2Y+ipWdJdAzGYgLuTAvOx2HC+jLMFAp/SFjpzWkjtcIzrhMesycYmejqXN8hGKLyIQ4VVrq4yhy9hQENUdnrwPR7dmpvwxMid2wTsbgMTcKb3kntqajBXHzVjv38R6g8MqVP1rJ6WZNoTaINwmUW4W6VyoHI3HNyfHuCkcL182GjXnhSJssEfI5fq8wxM41ZJ0pKwiDyISzfrlWuwBpcJELottQ3MdG5be6kgHx4qRHW+qPyY/ewRH9I9FaRuf8SCbKI2kR7dVSohKfSk8cylY9OhJaFrczxnrD0olx57zlW7JnIAV+WzDMOA1SqlgEige/Gp62sX4q80kem9tIqvKOE663LCUKUQq3h544K1QF8qON19ykiVxFx+RehvOzdnJIxTaJAWJLahya6RRhdh5BXepZEl5fwDRTFf3MSsCYt3IxzOrpoJ+KRPRgzJfd7Pc7Ita/XEn1EGHtPJd5nIe/A9H5XybiqyMf0LJ6XKT07ZGTvspAxqCOLdVOcLRJ4/kC0QZSX/XBuG+NbnPqLSt+uPJFkXpAl160iLp6ESVUWG6tlnd4dMUkcFIuXaTMSfRInlozrJUrnl4ovsvDSTySYVaLyKKIs9NwhYtw8keBWi494ZOLvi6Nb2Ux3wURR/hYrpFki8d4gku9ScXWyA0cJoxuCF4oTZko/RSip2iXXaHEqXZvrW4ozxZ3DWOFmm0WcXIaX27+v/y1OE4oru6ffLS7W5uEmnI4LhZ7Wi5P0ErxeHJgpfbnWS7sRVfo/x+Vty95t0ZToglZ43khPu9Ed/sToSjSJkLuuLcaOoOxqtwcVVPPwvUKPpXcBepTeBehRehegR+ldgH8BUHr5FK3WgcEAAAAASUVORK5CYII=\n", "text/latex": [ "$$\\displaystyle \\frac{h e^{- \\frac{h}{\\tau_{m}}}}{C}$$" ], "text/plain": [ " -h \n", " ───\n", " τₘ\n", "h⋅ℯ \n", "──────\n", " C " ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "PA_32 = PA.row(2).col(1)[0]\n", "sp.limit(PA_32, tau_s, tau_m)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**2.** The Taylor-series up to the second order of the function $PA_{32}(\\tau_s)$ is:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAArCAYAAACHItq+AAASG0lEQVR4nO2df7QV1XXHPw/jDwSDikapRH5YRZDHuyApBH9UYwRT0wCaxB+N+IzNSrWVsFwxamq60NViU0N9WqUu26QkNS1tQ0CJMRJQQ2JKBfRpYmljlReTCqLyS0UIVPrH94xz7ryZe+e+e2fuve/uz1qz5s6ZuXP33Dkz++x99tmnbf78+RgtyQeBfwQ+AOwDbgW+W1eJWo+PAzOA62L2LQHmA/+Vp0CGYTQP76u3AEbd2A98AXgWKfENwA+A3fUUqsWYgP7/KAOAkzHlbRhGCQbUWwCjbmwmVB5bge3AMfUTpyXpAE4Cfgr8Gpjkyk9HDaqjgG6kyHe5z8tzltEwjAbFFHj/5avAIymPnQwcDPzKKzsKeBUpmFbmO8D1GZ17ArAJmAZ0AbNd+QxgJWpUFYC5To4CMCsjWQzDaDLMhd5/KQDPpDhuKPAt4GrggFf+ZeRSf7HmkjUXtwFPAF8Hdlb43Z8nlE8HtqFG0/2ubB/whvt8HnCvd/x44D8r/G3DMPo5ZoH3XzqI71/1ORRYBtyO3LgBhwOfA/4+G9FqymIU7JUVzyEr+TN9+O74hOUVt17vHduOFP4R6Lnc7u0bhylwwzAiNIMCr8QVbIjjgeOQVfcD4G1kSX/EO6YNKb/HUDS6z4XAu8BPqpRjObLq45Y5VZ47T1keAi6PlA0Hvo0s6R3AUhQMmJYO4Gfe9gS3fR7weOTY41DMQsBIYBVwA/APwB85GYdW8PuGYTQ5zaDAC5S3JEHW4p3ZitI0THTrucDXkLJ4HljoHXMGcAnqU+12S7vbdyYKovJd6n3hamAYUjjvAp9y28OAf67y3HnK8h/Ah4CBbns0+n9eRv/jOSgAcFEF8gQKG/QcfgDFHEwHHo0c+whqIJzmtjuAfwHuAI4F7kMelEoaEIZRLTcD61CA5WvACuRZMnKiGRR4GlfwADSm9sHsxWkKCuihugRZav+DFMCx3jE/Qf9bwVsChTISuXmr5Q1gC3C0+601bnsL8g7kSTWyvIL6q3/Lbd+H+sRvBjaixs8CZD2n5QvIagY1KEa6z78DPBU5dhFqNDzvtjvQfR1CONTst9F9bnY+gboRjMbnHFQ3pyHv3n5UL4+uo0z9lYtQ12YRcQr8emR5RV2G9SCNKxhgCnrBBi7f5TSO67YeFJBy8N2uo0n/gj8M2BNTPp/k/zVYzon53gSkKLem/P1SfBl4y1v+IKbsrBLf74ss77j1QGAEcD4wL/Kby4HfVHDOJCYD/1fmmFFAD8VW/PvIv1FUa6ahZ/7NhP3fQvdtUG4StTano2f66oT9M1AXzs9RPbwCGQln5CJda7EKxSoVEafAg7GoT2cqTjrSuIIBZgIPoxYgZOe6PQ69XO+O2Xcy+oM3IHfSPrdeBVxLTOspQwrAv0fKJpIuKh3gdTSMLMo9wNgyS9R6hOSEJSvQ8Km1wAuoIbYMKacvJsh2H8Veg4diytbHfK+cLKDAvWfQy2ipVx5YFK+hOrjTncf/zXakfPPgKvRi/THwTVfWmdNvp2E08Bfov9wG7EX3dDH6/+I4CLiJ5MDJySiQ8C9RQ97Ing2oYfrnwOAUxwcBmNsylKlV2YXeS5/0C+OGkU1CD8gvchCqHAVCV3BgTS5FLwefmcCfetvBcJyJhO7SWlh/M935lnllbSgN6Y3AIagv8jvoJT8C9WmeB1xGacuwVhyO3KlRZT3RyZWGZ4hXCK+7pVJGA/8dUz4eKfF5aDjVQuD3UH/yw6jRFmUbxS+IN912Wu9CkixHAX+CrI79wJHevnbkRn8VNcwGu89JlmKr0gbcgp7FQ4AfIa/Y2+hZnoM8Jp8HvhH57hzUn/puwrkXoHfB39ZaaKMkt6MYkLnoHpSiC3Unrc1WpJZlCdJly3HGatQCHwSMQTch6UGqhvlU5oItUN4VfAqytKOBP1Bb1y0o0cYb6E8M+AbwFSdTAbmPrkHWxGUoWvkWKu+f7CTZJV2KwMJ5zisb6uRIa4E/iqzpWmVmG4AaM8MJ61zQWg/GQe9BY593uc+VjrmuRhYIlfZXkcLe4e07G3XhgF5O24EHUKPoJORSvxcpsGalk77Vt4A29CzchiyFU4FzUZ73m4AL0P8EuucTI9+/EQXmxXEK8FHgXwm7M/oznVR3L2rJUyjO4vPIS5LEHeg5+RTlu4CMvvEW0iMXBQVRBV5wZU+jsacPIAX4FrIspySc+GIUKfs66gd8AfVLRm94pS7YAuVdwTOB1cS71Wrpuh2C+t5XEFbQ69HDthGYmvBb7yCPwbUJ5601Heia/P9jInJjbkx5jp+hVvelNZLpLifDJhSrALK+/Xvtb48nOQlKFrKArOl2dA+XoDoN6veeBfyd294OfAw1QB5Hjd07UBa7aqP2m5mb0LOwAXma4jx4q5EFfRDyugQUUN/3Cwnn/ixqICQp+Go5Bhks5YyLvSg+pFHJ6jqWACeiRlQcC5EH5TySDZVhqFE3pA7yV8po4rsQ48hb5vXAp4ONqAv9dLce4Q5cjQJHTkOuzYeQezZwHR6ElPyl6Mb9mxP0Y0hpjQGu9M5fiQs2rSt4Jupbi6OWrtsLkVswmLFrGLrG/egPLedOzctyuM8tPquovPLcCvwNeuFW26J+DNUpn/EUewlGAS+5z35wVjk6ayALKIbhBVTfJ6MkN6B4irUUuwXX0zuQspUZherLHmSBxQVABjyKrPIPe2W/ixpCSQ2gj6I6mJVrdjDyHASMRO+tDcD3vPLXKH1t9Sar63jSrc+nt6fzbtTFeS6lkw3d7mS5iuT3dSPchyNR18+rqN7tKHN83jI/g/I/AL0VeBDANhWNBfYD2ZYi072AgmdA1sylKLDkK4RBZDeg9JNzkEuyL1mk0riCj0XW88XE47tLX0EtpVKu28Eku25nI6v2h257HlKKQRRmf+NR5DEZDvwyg/O3owYiqDG0hfAF3k4YnJUXt6B6/zby+gTW3m+In+7TCPki8mYsQp6NUgT59n1LrB1N5hLHIPTO2Uh2wWs9FGfz+0P0El5CfGO+Uekhm+tY59ZnR8oXocDCWSgO5XhXHozMqJQe6n8fdqDrWoDmIzif0t15PeQr8/8iQ/N4YEuSAu+kdxR64H4NLLkpyC38IBoP67MPvYCnueP6qsDLuYJ/H1WuVxPOcRdSsJuQct5LvOv2Lu9znDI+DPXhPULYiprl1tEsZv2JuGj7WjHX+7yZYovsSvIn6TfvTyg3RBvhJCwPpDg+yBbne+I+SHLO/ROQp29zwv4sKLh1d46/mQUFt+6u8jw70XvvxEj5NW69OlJ+K73TG3dSubes4NbdFX6vWoLhWguQITMdGXhpKLh1d21Feo+gMTGciAI/DPVBv0R86tLRbh08aNehh3c38bmog4w8fQ3sSeMKnknp6RVr5bqdjhoAQfT5YBRYcwD1FVdLT4ycAY/HlH2Txho2ZDQXPdSuvgXDM/eT7qU11a19A2EQyS/IQOFvL3HOHpKvJ45vUzq3feD9S5MBslp6yO7Zr+V1bEPDaH2yDtqs9X3oobJ6AjJAv466htKQdd0JnpPBUOxC73DbKxO+OAlp/8BFNt2tLyvzgy9XLmNqnqTycd19cd3ORq7Uh912kNFsF2rAVEsXxcOWQC25mU6ensi+7hr8ptG6dFG7+haMVHiTsAstiTY0jAzCWBJQ//bBvQ8HwtiRUjEcL1JZ/2KpLINtqCH/Cuq3zJousnn2a30dA8l3BEAW96GSejIAxcVA+ritPOpO8Jzsg2IFHgSwbYj50hHoYtYgZXcYUmJrUABKvfirPnynUtftQchV/xih+yJw6w90+6sN8uqKKeskDNB7osrzG4ZPV0xZJ32rbzvc+kgUeFqqQXs5CojdSHHa4+3A+xO+EwwBLTVRSyUpbMsxysnyZLkDa0RXTFkn1T/7tbyOAej+lotvqCVZ3Ie09WQASih0MgrMThsDk0fdCWJHtkPxMLKg/ztOgQcJUYJ9geukVuOEG5mz0cvDT96yFbWMD6H8WM28882XG85gS+ssefAyCqxpIxznHccpKDhoPwr08fNMvERvKzRgM7JmxlQraEpOdetSganD0f/7adSw341cpmOQy/VJV7bWHVsP0lxHG/AlNFJnD3qvxSV7GuOO7a6hfOWopfyVMADlM7jKnetyynuWAvKoO0Pc9zcFwgZMQm7iuB8PrPOg3+od1I88Dm9QeYQzKT3wv1m4CL1sohOldLn1IsIb59OGhtMtyUyyeNpsSbWciIaLbER1+eIGkKnWS17c6dZ/TTjhi8/H0YtpEEpX+9PI/nUo/iaOA8jTdwwaVpo1gSegVNBSwa2vRQFbUwhHpCxAo3DOQP/FvCyETEGa67gBKaprkQL5BOEoG58gbuHxKuRZjO5lZ8rjayl/JQxCcVLfRd3DaZU35FN3xqI4rXcgdKEf4oR+lvgJEeLc6zegPuGlKLjsOdQgOMEdfzC9oxabkVnohRONdL8b9Xd8Fv2hK1Hyin0oqvZM1IL6p7wENSpiP5oR7Fk0DecGlG2tFjENrcad6MX0GZS160FkmQcTW4xDsSazUf6FKD9EVk9SP+tS1MCaQfYzrgUJaOahHPjr6B1n00GY4jl4L6x026cSpvr9EYqzqQdpruMC4PuEMUG/JH6s/XTUTfhgFfIExmJahVhL+SvhTeRq301lyhvyqTsTCbNCvventiOFmzSByeloXJ+fKWklGia21H1/LpqNZixS6FclnKuZ+BBSwsti9h1AST4uRAP2C8AfI/fgBDSWuJPiPnejcdhMGCm6FfUptUKXUBa8i579T6Lc5xcg1+ZMlHr4emQhxSlvUJDQ91GjN46l6GWXx0yCT6OcAPvQizhugpoCuha/UT8CWW1+nv4Tybff2CfNdSxz+1ahd1dc/R+CjJjvEY7h7wvtSDk+XO5AR63k7ws76dvMfnnUnbPwEuG0zZ8/vw9ytgwL0Bj30dTvQTSyZzIazz+O/PqOjWJOQ0r/yoT9N6PncRLpc/pnxS+QB+4er+xXKJnVYq9sB2rk+zPbNRonoYbWHBSENZXilMvXoWs9mzCBV6UciRpyC9E9riXl5G80qqk7Y9Fz8F5DNu8Aq2ZjNrLSTHnnw83I7bQLBS6tIMwnEKUDjef9NQpg2YSU8GkV/uZQlD71akx515PnkeVzQsL+O5Fb/raE/XkxCCkNvxFxDL0nCxqFrNd6NzbK8SKKXZiM4iYmePsGomdyKX1X3iCrcZ/7nVpTSv5Go9q68zngz/wCU+ClGUsYdGBkzzkoKHAayjW+H7nHjo4cNwf1We9FfUenELZKK+myOBS54W6nd2CVkT83oS6oOPYgN/169CKsF4GCeNYrm4Tqop9xskBx3oxG40bUxTcOPT+3oiDmJ7xjRqJMhEkTPKVlBQrU2lLleXzSyN9oVFN3ggj1Hv+EcfOBG0a9mBHZvgJV5DMI+0+noYCnGwijn0HW2Y8Jlf1y5FqL40pkrS9GQzn6czrcZmI3mvXt/cRH8q6heCrfetCBYoH8XN8T0eidfZHjumlcr86hqME0AuW1WIuCt/y+2Y3EZ9lsBNLI32hUU3d2EpP10/rAjUZmGMpqdCZhcoSnkDUWnVghylAUmHkoGmd8CQqwAvXHTUHKwE+rewXpZ0EzDMOoK2aBG41MF2qJBkNDTkUjA9LMU/6GW09EXUVrCLN6gZS5dSEZhtG0mAI3GpU7kJV9FmGq2oJbx2ULTGIC6nvbWu5AwzCMZsIUuNGILERJQc6lOHHHQLeuZK7hCeQzq5RhGEaumAvRaDTuJlTe0XnkgzS/Sf3fh8eUjUa5kg3DMPoVpsCNRmIRGhpyGcpKdLxbBrv961DGrnvccSej/NiXoHSKHTHnHIAiVYdTXN9XoD72tSgydAoaUtZD9cNmDMMwMscUuNFIXIOmrl2NUp0Gi69QLwK+htJzdiOl/iUUpR6XNOMuFMi2ieI5p8cjC38qmqRhIRpe9hGU1MUwDKOhsT5wo5FIM4PWXjQPfNq54B9DFrjPEajxer/b3gPci8YeDyac990wDKNhMQvcaEXGo/HkcdvjKT2fr2EYRkNgCtxoRcZTnMBlFEr2Aopat2QuhmE0PKbAjVaknVCBD0PjxA94+0yBG4bR8FgfuNGK+BOebAY+7G0nTWdpGIbRUJgFbhiGYRhNiClwwzAMw2hCTIEbhmEYRhPy/wISNn8kqiH2AAAAAElFTkSuQmCC\n", "text/latex": [ "$$\\displaystyle \\frac{h e^{- \\frac{h}{\\tau_{m}}}}{C} + \\frac{h^{2} \\left(- \\tau_{m} + \\tau_{s}\\right) e^{- \\frac{h}{\\tau_{m}}}}{2 C \\tau_{m}^{2}} + O\\left(\\left(- \\tau_{m} + \\tau_{s}\\right)^{2}; \\tau_{s}\\rightarrow \\tau_{m}\\right)$$" ], "text/plain": [ " -h -h \n", " ─── ─── \n", " τₘ 2 τₘ \n", "h⋅ℯ h ⋅(-τₘ + τₛ)⋅ℯ ⎛ 2 ⎞\n", "────── + ────────────────── + O⎝(-τₘ + τₛ) ; τₛ → τₘ⎠\n", " C 2 \n", " 2⋅C⋅τₘ " ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "PA_32_series = PA_32.series(x=tau_s,x0=tau_m,n=2)\n", "PA_32_series " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Therefore we have \n", " \n", "$T(PA_{32}(\\tau_s,\\tau_m))=PAs_{32}+PA_{32}^{lin}+O(2)$, where $PA_{32}^{lin}=h^2(-\\tau_m + \\tau_s)*exp(-h/\\tau_m)/(2C\\tau_m^2)$\n", " \n", "**3.** We define\n", "\n", "$dev:=|PA_{32}-PAs_{32}|$\n", " \n", "We also define $PA_{32}^{real}$ which is the correct value of P32 without misscalculation (instability).\n", " \n", "In the following we assume $0<|\\tau_s-\\tau_m|<0.1$. We consider two different cases\n", " \n", "**a)** When $dev \\geq 2|PA_{32}^{lin}|$ we do not trust the numeric evaluation of $PA_{32}$, since it strongly deviates from the first order correction. In this case the error we make is\n", " \n", "$|PAs_{32}-PA_{32}^{real}|\\approx |P_{32}^{lin}|$\n", " \n", "**b)** When $dev \\le |2PA_{32}^{lin}|$ we trust the numeric evaluation of $PA_{32}$. In this case the maximal error occurs when $dev\\approx 2 PA_{32}^{lin}$ due to numeric instabilities. The order of the error is again\n", "\n", "$|PAs_{32}-PA_{32}^{real}|\\approx |P_{32}^{lin}|$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The entry $A_{31}$ is numerically unstable, too and we treat it analogously." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tests and examples\n", "We will now show that the stability criterion explained above leads to a reasonable behavior for $\\tau_s\\rightarrow\\tau_m$" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Aug 24 12:28:16 CopyFile [Error]: \n", " Could not open source file.\n", "Error in nest resource file: /BadIO in CopyFile_\n", "\n", " -- N E S T --\n", " Copyright (C) 2004 The NEST Initiative\n", "\n", " Version: 3.3\n", " Built: May 20 2022 17:08:28\n", "\n", " This program is provided AS IS and comes with\n", " NO WARRANTY. See the file LICENSE for details.\n", "\n", " Problems or suggestions?\n", " Visit https://www.nest-simulator.org\n", "\n", " Type 'nest.help()' to find out more about NEST.\n", "\n" ] } ], "source": [ "import nest\n", "import numpy as np\n", "import pylab as pl" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Neuron, simulation and plotting parameters" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "taum = 10.\n", "C_m = 250.\n", "# array of distances between tau_m and tau_ex\n", "epsilon_array = np.hstack(([0.],10.**(np.arange(-6.,1.,1.))))[::-1]\n", "dt = 0.1\n", "fig = pl.figure(1)\n", "NUM_COLORS = len(epsilon_array)\n", "cmap = pl.get_cmap('gist_ncar')\n", "maxVs = []" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Loop through epsilon array" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Aug 24 12:28:21 SimulationManager::set_status [Info]: \n", " Temporal resolution changed from 0.1 to 0.1 ms.\n", "\n", "Aug 24 12:28:21 NodeManager::prepare_nodes [Info]: \n", " Preparing 3 nodes for simulation.\n", "\n", "Aug 24 12:28:21 SimulationManager::start_updating_ [Info]: \n", " Number of local nodes: 3\n", " Simulation time (ms): 200\n", " Number of OpenMP threads: 1\n", " Number of MPI processes: 1\n", "\n", "Aug 24 12:28:21 SimulationManager::run [Info]: \n", " Simulation finished.\n", "\n", "Aug 24 12:28:21 SimulationManager::set_status [Info]: \n", " Temporal resolution changed from 0.1 to 0.1 ms.\n", "\n", "Aug 24 12:28:21 NodeManager::prepare_nodes [Info]: \n", " Preparing 3 nodes for simulation.\n", "\n", "Aug 24 12:28:21 SimulationManager::start_updating_ [Info]: \n", " Number of local nodes: 3\n", " Simulation time (ms): 200\n", " Number of OpenMP threads: 1\n", " Number of MPI processes: 1\n", "\n", "Aug 24 12:28:21 SimulationManager::run [Info]: \n", " Simulation finished.\n", "\n", "Aug 24 12:28:21 SimulationManager::set_status [Info]: \n", " Temporal resolution changed from 0.1 to 0.1 ms.\n", "\n", "Aug 24 12:28:21 NodeManager::prepare_nodes [Info]: \n", " Preparing 3 nodes for simulation.\n", "\n", "Aug 24 12:28:21 SimulationManager::start_updating_ [Info]: \n", " Number of local nodes: 3\n", " Simulation time (ms): 200\n", " Number of OpenMP threads: 1\n", " Number of MPI processes: 1\n", "\n", "Aug 24 12:28:21 SimulationManager::run [Info]: \n", " Simulation finished.\n", "\n", "Aug 24 12:28:21 SimulationManager::set_status [Info]: \n", " Temporal resolution changed from 0.1 to 0.1 ms.\n", "\n", "Aug 24 12:28:21 NodeManager::prepare_nodes [Info]: \n", " Preparing 3 nodes for simulation.\n", "\n", "Aug 24 12:28:21 SimulationManager::start_updating_ [Info]: \n", " Number of local nodes: 3\n", " Simulation time (ms): 200\n", " Number of OpenMP threads: 1\n", " Number of MPI processes: 1\n", "\n", "Aug 24 12:28:21 SimulationManager::run [Info]: \n", " Simulation finished.\n", "\n", "Aug 24 12:28:21 SimulationManager::set_status [Info]: \n", " Temporal resolution changed from 0.1 to 0.1 ms.\n", "\n", "Aug 24 12:28:21 NodeManager::prepare_nodes [Info]: \n", " Preparing 3 nodes for simulation.\n", "\n", "Aug 24 12:28:21 SimulationManager::start_updating_ [Info]: \n", " Number of local nodes: 3\n", " Simulation time (ms): 200\n", " Number of OpenMP threads: 1\n", " Number of MPI processes: 1\n", "\n", "Aug 24 12:28:21 SimulationManager::run [Info]: \n", " Simulation finished.\n", "\n", "Aug 24 12:28:21 SimulationManager::set_status [Info]: \n", " Temporal resolution changed from 0.1 to 0.1 ms.\n", "\n", "Aug 24 12:28:21 NodeManager::prepare_nodes [Info]: \n", " Preparing 3 nodes for simulation.\n", "\n", "Aug 24 12:28:21 SimulationManager::start_updating_ [Info]: \n", " Number of local nodes: 3\n", " Simulation time (ms): 200\n", " Number of OpenMP threads: 1\n", " Number of MPI processes: 1\n", "\n", "Aug 24 12:28:21 SimulationManager::run [Info]: \n", " Simulation finished.\n", "\n", "Aug 24 12:28:21 SimulationManager::set_status [Info]: \n", " Temporal resolution changed from 0.1 to 0.1 ms.\n", "\n", "Aug 24 12:28:21 NodeManager::prepare_nodes [Info]: \n", " Preparing 3 nodes for simulation.\n", "\n", "Aug 24 12:28:21 SimulationManager::start_updating_ [Info]: \n", " Number of local nodes: 3\n", " Simulation time (ms): 200\n", " Number of OpenMP threads: 1\n", " Number of MPI processes: 1\n", "\n", "Aug 24 12:28:21 SimulationManager::run [Info]: \n", " Simulation finished.\n", "\n", "Aug 24 12:28:21 SimulationManager::set_status [Info]: \n", " Temporal resolution changed from 0.1 to 0.1 ms.\n", "\n", "Aug 24 12:28:21 NodeManager::prepare_nodes [Info]: \n", " Preparing 3 nodes for simulation.\n", "\n", "Aug 24 12:28:21 SimulationManager::start_updating_ [Info]: \n", " Number of local nodes: 3\n", " Simulation time (ms): 200\n", " Number of OpenMP threads: 1\n", " Number of MPI processes: 1\n", "\n", "Aug 24 12:28:21 SimulationManager::run [Info]: \n", " Simulation finished.\n" ] }, { "data": { "text/plain": [ "Text(0, 0.5, 'voltage V (mV)')" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "for i,epsilon in enumerate(epsilon_array):\n", " nest.ResetKernel() # reset simulation kernel \n", " nest.resolution = dt\n", "\n", " # Current based alpha neuron \n", " neuron = nest.Create('iaf_psc_alpha') \n", " neuron.set(C_m=C_m, tau_m=taum, t_ref=0., V_reset=-70., V_th=1e32,\n", " tau_syn_ex=taum+epsilon, tau_syn_in=taum+epsilon, I_e=0.)\n", " \n", " # create a spike generator\n", " spikegenerator_ex = nest.Create('spike_generator')\n", " spikegenerator_ex.spike_times = [50.]\n", " \n", " # create a voltmeter\n", " vm = nest.Create('voltmeter', params={'interval':dt})\n", "\n", " ## connect spike generator and voltmeter to the neuron\n", " nest.Connect(spikegenerator_ex, neuron, 'all_to_all', {'weight':100.})\n", " nest.Connect(vm, neuron)\n", "\n", " # run simulation for 200ms\n", " nest.Simulate(200.) \n", "\n", " # read out recording time and voltage from voltmeter\n", " times = vm.get('events','times')\n", " voltage = vm.get('events', 'V_m')\n", " \n", " # store maximum value of voltage trace in array\n", " maxVs.append(np.max(voltage))\n", "\n", " # plot voltage trace\n", " if epsilon == 0.:\n", " pl.plot(times,voltage,'--',color='black',label='singular')\n", " else:\n", " pl.plot(times,voltage,color = cmap(1.*i/NUM_COLORS),label=str(epsilon))\n", "\n", "pl.legend()\n", "pl.xlabel('time t (ms)')\n", "pl.ylabel('voltage V (mV)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Show maximum values of voltage traces" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fig = pl.figure(2)\n", "pl.semilogx(epsilon_array,maxVs,color='red',label='maxV')\n", "#show singular solution as horizontal line\n", "pl.semilogx(epsilon_array,np.ones(len(epsilon_array))*maxVs[-1],color='black',label='singular')\n", "pl.xlabel('epsilon')\n", "pl.ylabel('max(voltage V) (mV)')\n", "pl.legend()" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "pl.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The maximum of the voltage traces show that the non-singular case nicely converges to the singular one and no numeric instabilities occur. " ] }, { "cell_type": "markdown", "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.9.13" } }, "nbformat": 4, "nbformat_minor": 4 }