{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# FHN bacteria model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "The model consists of two variables $V_m$ and $W$. Their dynamics are expressed in the following ordinary differential equations.
\n", "\n", "\n $\\cfrac{dV_m}{dt} = k_{K}((V_m + V_{m,0}) - \\alpha (V_m + V_{m,0})^3 + W ) + \\cfrac{dI_v}{dt}$\n", "\n", "\n $\\cfrac{dW}{dt} = (-(V_m + V_{m,0}) +\\beta - W ) + \\cfrac{dI_w}{dt}$\n", "\n", "
\n", "where $\\beta$ is defined by:\n", "
\n", " $\\beta = 1 - 0.1 \\log{k_K}$\n", "
\n", "$V_m$ is membrane potential, $W$ is recovery variables. $V_{m,0}$ is the offset in $V_m$. $\\alpha$ is a parameter for $V_m$ dynamics property.
\n", "$k_K$ corresponds to the degree of $K^+$ concentration gradient across the membrane. This determines the time scale of the dynamics and the resting membrnae potential.\n", "$I_v$, $I_w$ are the strengths of the externally applied electrical field.
\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.integrate import odeint\n", "import os\n", "import seaborn as sns\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#parameters in the differential equation of w. \n", "#While these parameters were included in the script, \n", "#they are not effective in our simulations as they are multiplificaiton terms set as 1.\n", "a = 1.\n", "b = 1.\n", "c = 1.\n", "\n", "#This parameter correponds to /alpha in the equation. \n", "alpha = 10.\n", "\n", "#This parameter corresponds to k_K in the equation.\n", "k = 10.\n", "\n", "#offset in Vm. Note that this parameter does not affect the dynamics. \n", "vm0 = 1.5\n", "\n", "#time step size\n", "dt = 0.001\n", "\n", "#time scale of the events\n", "tscale = 60.\n", "\n", "#max time to simulate\n", "Tmax = 10.\n", "\n", "tvec = np.arange(0, Tmax, dt)\n", "\n", "\n", "#duration of electrical stimulation\n", "ees_duration = 2.5\n", "\n", "#tvec for electrical stimulaiton\n", "t_ees = np.arange(0,(ees_duration*(Tmax/dt)/(Tmax*tscale))*dt,dt)\n", "\n", "#initial parameter for v and w\n", "vw0 = [-.5,-.5]\n", "\n", "#external stimulation in v and w.\n", "#v is membrane potential and w is recovery variables\n", "Iv = 0.01\n", "Iw = -0.075\n", "I1 = np.array([Iv, Iw])/dt\n", "\n", "#parameter for I when no electrical field is applied\n", "I0 = [0,0]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#ODE\n", "def f(vw,t):\n", " v, w = vw\n", " dvdt = k*((v+vm0) - alpha*(v+vm0)**3 + w) +Iv\n", " dwdt = a*(-(v+vm0) + b - c*w) +Iw\n", " return [dvdt, dwdt]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Simulation for proliferative cells" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#parameter k_K for prolieferative cells\n", "k = 10.\n", "\n", "#calculation of beta \n", "b1 = 1-0.1* np.log(k)\n", "\n", "b = b1\n", "\n", "#simulation without electric field. This is to bring the system to the equlibrium. \n", "Iv, Iw = I0\n", "vwout1p = odeint(f, vw0, tvec)\n", "\n", "#simulation with electric field. \n", "Iv, Iw = I1\n", "vwout1pe = odeint(f, vwout1p[-1,:], t_ees)\n", "\n", "#simulation after removal of electric field.\n", "Iv, Iw = I0\n", "vwout2p = odeint(f, vwout1pe[-1,:], tvec) \n", "\n", "#simulation result combining above.\n", "vwoutp = np.concatenate([vwout1p, vwout1pe, vwout2p])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Simulation for inhibited cells" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#parameter k_K for inhibited cells\n", "k= .1\n", "b2 = 1- 0.1* np.log(k)\n", "\n", "b = b2\n", "\n", "Iv, Iw = I0\n", "vwout1i = odeint(f, vw0, tvec)\n", "Iv, Iw = I1\n", "vwout1ie = odeint(f, vwout1i[-1,:], t_ees)\n", "Iv,Iw = I0\n", "vwout2i = odeint(f, vwout1ie[-1,:], tvec) \n", "vwouti = np.concatenate([vwout1i, vwout1ie, vwout2i])\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Ploting the simulation results. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#time \n", "time = np.linspace(0, (t_ees.size + 2* tvec.size)*dt*tscale, t_ees.size + 2* tvec.size) - tvec.size*dt*tscale\n", "\n", "#Initial frame for ploting. The system is at the equlibrium (before electrical stimulaiton). \n", "frameq = np.int(tvec.size*0.9)\n", "\n", "#make a figure\n", "fig, ax = plt.subplots(figsize = (5, 4))\n", "\n", "#indicate the electrical stimulation window\n", "plt.fill_between([0,ees_duration],-4, 4, color ='yellow', alpha = 0.4)\n", "\n", "#plot the simulation results\n", "plt.plot( time[frameq:], vwoutp[frameq,0] - vwoutp[frameq:,0],'-', c = 'b', lw = 4, alpha = 0.6)\n", "plt.plot( time[frameq:], vwouti[frameq,0] - vwouti[frameq:,0],'-', c = 'r', lw = 4, alpha = 0.6)\n", "\n", "plt.xticks(np.arange(0,55,10),size = 16)\n", "plt.yticks(np.arange(-0.5,1.2,0.5),size = 16)\n", "\n", "plt.ylim(-0.5,1.3)\n", "plt.xlim(-5, 45)\n", "\n", "plt.ylabel(r'$-\\Delta V_m$', size = 18)\n", "plt.xlabel(r'time (sec)', size = 18)\n", "\n", "sns.despine()\n", "plt.tight_layout()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Plotting the simulation results in $V_m$-$W$ phase space with nullclines" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "v = np.linspace(-3, 0, 100)\n", "\n", "nulc_v = alpha*(v+vm0)**3 - (v +vm0)\n", "nulc_w = (-a*(v+vm0) + a*b1)/(a*c)\n", "\n", "nulc_wi= (-a*(v+vm0) + a*b2)/(a*c)\n", "\n", "fig, ax = plt.subplots(2,1, figsize = (5.3, 8), sharex=True, sharey=True)\n", "\n", "ax1 = ax[0]\n", "ax2 = ax[1]\n", "\n", "ax1.plot(v, nulc_v, '0.3', lw =2, ls = '--')\n", "ax1.plot(v, nulc_w, '0.3', lw =2, ls = '--')\n", "ax1.plot(vwoutp[9000:,0], vwoutp[9000:,1], c = 'b', lw = 3, alpha = 0.5)\n", "\n", "\n", "ax2.plot(v, nulc_v, '0.2', lw =2, ls = '--')\n", "ax2.plot(v, nulc_wi, '0.2', lw =2, ls = '--')\n", "ax2.plot(vwouti[9000:,0], vwouti[9000:,1], c = 'r', lw = 3, alpha = 0.5)\n", "\n", "\n", "plt.xlim(-2.45, -0.2)\n", "plt.ylim(-3,2.1)\n", "\n", "#ax1.set_xlabel(r'$V_m$ (au)', size = 16)\n", "ax2.set_xlabel(r'$V_m$ (au)', size = 16)\n", "ax1.set_ylabel(r'$W$ (au)', size = 16)\n", "ax2.set_ylabel(r'$W$ (au)', size = 16)\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.6.3" } }, "nbformat": 4, "nbformat_minor": 2 }