MFS Notebook.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.integrate as sp_int\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Set up the constants\n",
    "C_m  = 1 # membrane capacitance (uF/cm^2)\n",
    "g_Na = 112.5 # maximum conductances (mS/cm^2)\n",
    "g_K  = 225\n",
    "g_l  = 0.25\n",
    "g_d  = 0.39\n",
    "theta_m = -22   # half-maximum, a parameter in the original paper\n",
    "                # NB: theta_m is not given by parameters from supplemental information in paper\n",
    "E_Na = 50 # reversal potentials (mV)\n",
    "E_K  = -90\n",
    "E_l  = -70"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# equations governing opening/closing rates.\n",
    "def m_inf(V): return (1 + np.exp(-(V - theta_m)/11.5))**(-1)\n",
    "def h_inf(V): return (1 + np.exp((V + 58.3)/6.7))\n",
    "def tau_h(V): return 0.5 + 14*((1 + np.exp((V + 60)/12))**(-1))\n",
    "def n_inf(V): return (1 + np.exp(-(V + 12.4)/6.8))**(-1)\n",
    "def tau_n(V): return ((0.087 + 11.4 * ((1 + np.exp(V + 14.6)/8.6)**(-1))) *\n",
    "        (0.087 + 11.4 * ((1 + np.exp(-(V - 1.3)/18.7))**(-1))))\n",
    "def a_inf(V): return (1 + np.exp(-(V + 50)/20))**(-1)\n",
    "def b_inf(V): return (1 + np.exp((V + 70)/6))**(-1)\n",
    "tau_a = 2\n",
    "tau_b = 150"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# define the currents of the model\n",
    "def I_K(V, n): return g_K * (n**2) * (V - E_K)\n",
    "def I_Na(V, h): return (m_inf(V)**3) * h * (V - E_Na)\n",
    "def I_d(V, a, b): return g_d * (a**3) * b * (V - E_K)\n",
    "def I_l(V): return g_l * (V - E_l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# define function that is time derivative of voltage: f = dV/dt\n",
    "def f(t, y):\n",
    "    # set external current: make sure it is 0 for large times.\n",
    "    if 0 < t <= 10:\n",
    "        I_e = 3\n",
    "    else:\n",
    "        I_e = 0\n",
    "    # set noise current in input\n",
    "    D = 0\n",
    "    I_noise = np.sqrt(2*D*np.random.normal(0, 1))\n",
    "    # set the variables that are to be integrated\n",
    "    V, n, h, a, b = y\n",
    "\n",
    "    # define DV/dt ('_dot' denotes time differentiation)\n",
    "    V_dot = 1/C_m * (I_e + I_noise - (I_K(V, n) + I_Na(V, h) + I_d(V, a, b) + I_l(V)))\n",
    "\n",
    "    # enter the equations controlling the gating variables.\n",
    "    n_dot = (n_inf(V) - n) / tau_n(V)\n",
    "    h_dot = (h_inf(V) - h) / tau_h(V)\n",
    "    a_dot = (a_inf(V) - a) / tau_a\n",
    "    b_dot = (b_inf(V) - b) / tau_b\n",
    "\n",
    "    # since y = [V, n, h, a, b] we return [V_dot, n_dot, h_dot, a_dot, b_dot]\n",
    "    return [V_dot, n_dot, h_dot, a_dot, b_dot]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# enter intial values for V, n, m, h\n",
    "V_0 = -72\n",
    "n_0 = n_inf(V_0)\n",
    "h_0 = h_inf(V_0)\n",
    "a_0 = a_inf(V_0)\n",
    "b_0 = b_inf(V_0)\n",
    "y_0 = [V_0, n_0, h_0, a_0, b_0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/Thomas/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:2: RuntimeWarning: overflow encountered in exp\n",
      "  from ipykernel import kernelapp as app\n",
      "/Users/Thomas/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:5: RuntimeWarning: overflow encountered in exp\n",
      "/Users/Thomas/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:7: RuntimeWarning: overflow encountered in exp\n",
      "/Users/Thomas/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:8: RuntimeWarning: overflow encountered in exp\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHB1JREFUeJzt3X2QFfWd7/H3F4YHgSyIMyCPAiIqqLgyEIziBhejMVfJ\nJmbFXQ25m12yrPF6s8n1xrIqld2NtcbNeu/ummwViSZGcyVUFh+iUbNgSiOmwAFReVQUcHgQGECF\nmMwg871//Posx3HOPPU50326P6+qUz2nu0/394f4mR+/0/1rc3dERCT7+iRdgIiI9A4FvohITijw\nRURyQoEvIpITCnwRkZxQ4IuI5IQCX0QkJxT4IiI5ocAXEcmJmqQLKFZbW+sTJkxIugwRkaqydu3a\nJnev62y/VAX+hAkTaGhoSLoMEZGqYmY7u7KfhnRERHJCgS8ikhMKfBGRnFDgi4jkhAJfRCQnFPgi\nIjmhwBcRyQkFfkwvvACrVyddhYhI51J141U1mjUrLPVoYBFJO/Xwy6SlJekKREQ6psAvkwMHkq5A\nRKRjCvwYinv1hw4lV4eISFfEDnwzu8nMtpjZRjO7s2j9rWa2zcy2mtnlcc+TRsW9egW+iKRdrC9t\nzWwuMB+Y7u7NZjYiWj8VWABMA0YDK8xsirsfj1twmijwRaSaxO3hLwbucPdmAHffH62fDyx192Z3\n3w5sA2bFPFfq7N9/4mcFvoikXdzAnwLMMbPVZvaMmc2M1o8BGov22xWtyxT18EWkmnQ6pGNmK4BT\n29l0W/T54cBsYCawzMwmdacAM1sELAIYP358dz6aOPXwRaSadBr47j6v1DYzWwwsd3cH1phZK1AL\n7AbGFe06NlrX3vGXAEsA6uvrq+r2pQMHoKYGhg9X4ItI+sUd0nkYmAtgZlOA/kAT8CiwwMwGmNlE\n4AxgTcxzpcqOHfDqq1BXB6ecosAXkfSLO7XCvcC9ZrYBaAEWRr39jWa2DNgEvA/cmKUrdLZuhbPO\nCj/PmAEDByrwRST9YgW+u7cA15fYdjtwe5zjp9GxY7BgwYn3Y8dCays0Npb+jIhIGuhO22769rdh\n/Xr4ylfC+8su0xi+iFQHBX43bN4Mf//3oYd/111w8CDceKMCX0SqgwK/G/72b2HwYPjXfw3vhw8/\nsTx6VDNmiki6KfC76Mknw+sb3whX5hQrBP/hw71fl4hIVynwu6C1FW65BSZPDkM4bRUCX8M6IpJm\neuJVFzz6KLzyCjzwAPTv/+HtCnwRqQbq4XfCHW6/HSZNgmuvbX8fBb6IVAP18DuxYgU0NMD3vx+m\nUWjPySeHpcbwRSTN1MPvxN13w4gRcMMNpfdRD19EqoECvwM7d8Jjj8Ff/iUMGFB6v6FDwUyBLyLp\npsDvwJIlYbloUcf79ekThnUOHqx8TSIiPaXAL+H99+Gee+BTn4LTTut8/9paaGqqfF0iIj2lwC9h\n5UrYtw++8IWu7V9Xp8AXkXRT4Jfwk5+Esfkrr+za/rW1H3zkoYhI2ijw2/Hee/DQQ3DNNWGu+67Q\nkI6IpJ0Cvx2PPRYmQ/vzP+/6ZwpDOl5VD2kUkTxR4LfjoYdCgF9ySdc/U1sbHo7y7ruVq0tEJA4F\nfhstLfDEE3DVVdC3b9c/V5hBU+P4IpJWCvw2fv1reOcduPrq7n2utjYsNY4vImmlwG/jkUfCF7WX\nXda9z6mHLyJpp8Av4g4//3kI+0GDuvdZ9fBFJO0U+EW2bYMdO+CTn+z+Z9XDF5G0U+AXWbkyLOfN\n6/5nBw8OE6yphy8iaRU78M3sJjPbYmYbzezOaN1lZrbWzF6JlpfGL7XyVqyAcePCowy7y0w3X4lI\nusV6AIqZzQXmA9PdvdnMRkSbmoCr3H2PmZ0DPAWMiVdqZR0/Dk8/DZ/+dAjvnqir05COiKRX3Cde\nLQbucPdmAHffHy1fLNpnI3CSmQ0o7JdG69eHJ1b1ZDinQD18EUmzuEM6U4A5ZrbazJ4xs5nt7PNZ\nYF2awx5C7x7g0hiDT+rhi0iaddrDN7MVwKntbLot+vxwYDYwE1hmZpPcw4wyZjYN+DbwiQ6OvwhY\nBDB+/Pju1l82q1bBGWfAqe21tIvUwxeRNOs08N295CCHmS0GlkcBv8bMWoFa4ICZjQUeAj7v7q93\ncPwlwBKA+vr6RKYec4fnn+/6VMil1NWFu3RbWqB///LUJiJSLnGHdB4G5gKY2RSgP9BkZsOAx4Gv\nu/uqmOeouG3bwlDMxz4W7zi6+UpE0ixu4N8LTDKzDcBSYGHU2/8yMBn4hpmtj14jOjpQkp5/Piwv\nuijecUaODMv9++MdR0SkEmJdpePuLcD17az/FvCtOMfuTatWwbBhcPbZ8Y5TCPy33opfk4hIuelO\nW0IP/8ILoU/MP43CF7779sWvSUSk3HIf+EeOwKZNMHt2/GOphy8iaZb7wH/xxXCVTn19/GMNGRJm\n2VQPX0TSKPeBv3ZtWF5wQXmOd+qpCnwRSafcB/66dTB6dLwbroqNHKkhHRFJp9wH/tq1MGNG+Y43\ncqR6+CKSTrkO/KNHYcuW8g3nQPiXgnr4IpJGuQ78l14KX9iWu4d/8CAcO1a+Y4qIlEOuA3/durAs\ndw8fdLetiKRPrgP/5ZfD/DejR5fvmIVr8TWOLyJpk+vA37ABzj2350+4ao8CX0TSKreB7x4C/5xz\nynvcwpCOvrgVkbTJbeDv3Bmu0il34KuHLyJpldvA37AhLMsd+IMHhykW1MMXkbTJfeBPm1b+Y+ta\nfBFJo1wH/rhxMHRo+Y89ejTs3l3+44qIxJHrwC/3cE7BmDGwZ09lji0i0lO5DPz334fNmysznAMh\n8HfvDlcCiYikRS4Df+dOaGmBqVMrc/wxY6C5GQ4dqszxRUR6IpeB/+qrYTllSmWOX7hzV+P4IpIm\nCvwKGDMmLBX4IpImuQz8rVth2LAwj04lFAJfX9yKSJrkMvBffRXOPLO8c+gU05COiKRR7MA3s5vM\nbIuZbTSzO9tsG29mR83sa3HPU06vvlq54RyA/v2hrk6BLyLpUhPnw2Y2F5gPTHf3ZjMb0WaXu4An\n4pyj3N57DxobKxv4cOLSTBGRtIgV+MBi4A53bwZw9/967IeZfRrYDvw25jnK6rXXwrI3Al9j+CKS\nJnGHdKYAc8xstZk9Y2YzAcxsCPC/gb/r7ABmtsjMGsys4cCBAzHL6VzhCp0zz6zsedTDF5G06bSH\nb2YrgFPb2XRb9PnhwGxgJrDMzCYB3wT+j7sftU6+GXX3JcASgPr6+orfm1oI/MmTK3ueMWPCYw5b\nWsKYvohI0joNfHefV2qbmS0Glru7A2vMrBWoBT4KXBN9iTsMaDWz37v73WWqu8defx1GjQrTGFdS\n4UqdvXvhtNMqey4Rka6IO4b/MDAX+JWZTQH6A03uPqewg5l9EziahrAHeOMNmDSp8ucZNy4sGxsV\n+CKSDnHH8O8FJpnZBmApsDDq7afW9u0wcWLlz1MI+Z07K38uEZGuiNXDd/cW4PpO9vlmnHOUU0tL\n6HH3Rg9//PiwVOCLSFrk6k7bN98MUxb3Rg9/0KBw85UCX0TSIleBv317WPZG4EMY1lHgi0ha5Crw\n33gjLHsr8MePV+CLSHrkKvC3b4d+/U7MZllphR5+ur/GFpG8yF3gn3Ya9O3bO+c77TT43e/g4MHe\nOZ+ISEdyFfhvvNF7wzmgSzNFJF1yFfi9dQ1+gQJfRNIkN4F/5EgYWunNwNe1+CKSJrkJ/DffDMve\nnOZg+PAwZ48CX0TSIDeB39gYloU5bnqDWfgFs2NH751TRKQUBX6FnX76iev/RUSSlKvANzsxbXFv\nOf30MCWzrsUXkaTlKvBHjQo3XvWmyZPDc3Tfeqt3zysi0lauAr+3h3Mg9PAh9PJFRJKkwK+wQuBv\n29b75xYRKZaLwHdPLvALUzmohy8iSctF4B86FOa0KdwI1Zv69w/nVQ9fRJKWi8BP6pLMgsmT1cMX\nkeQp8HtB4dJMEZEk5SLwC9MqJBn4hw7B4cPJnF9EBHIS+I2N4fr7kSOTOf/kyWH52mvJnF9EBHIS\n+Lt3hzts+yTU2rPOCsutW5M5v4gIlCHwzewmM9tiZhvN7M6i9eeZ2W+i9a+Y2cC45+qpvXt7f0qF\nYqefDjU1sHlzcjWIiNTE+bCZzQXmA9PdvdnMRkTra4AHgBvc/SUzOwU4FrvaHtqzB84+O6mzh+Gk\nyZMV+CKSrLg9/MXAHe7eDODu+6P1nwBedveXovUH3f14zHP1WNI9fAi/cBT4IpKkuIE/BZhjZqvN\n7Bkzm1m03s3sKTNbZ2a3xDxPj/3ud/D222HitCSdfXa4+aqlJdk6RCS/Oh3SMbMVwKntbLot+vxw\nYDYwE1hmZpOi9RdH694DVprZWndf2c7xFwGLAMZX4FbYvXvDMunAP+ssOH48XI+f5PCSiORXp4Hv\n7vNKbTOzxcByd3dgjZm1ArXALuBZd2+K9vsFcAHwocB39yXAEoD6+vqyzxpfCPw0DOlAGNZR4ItI\nEuIO6TwMzAUwsylAf6AJeAo418wGRV/g/hGwKea5emTPnrBMQw8fNI4vIsmJdZUOcC9wr5ltAFqA\nhVFv/7CZ3QW8ADjwC3d/POa5eiQtPfwhQ8Kdvgp8EUlKrMB39xbg+hLbHiBcmpmoPXvCZZGnnJJ0\nJTB1KmzcmHQVIpJXmb/Tdu/eMJxjlnQlMH06bNoExxK7I0FE8izzgb9nT/Lj9wXnnRcuy9yyJelK\nRCSPMh/4abjpqmD69LB86aVk6xCRfMpF4Kelh3/mmeEJWAp8EUlCpgP/978P89CnpYffrx9MmwYv\nv5x0JSKSR5kO/LfeCsu09PAhDOuohy8iSch04O+PpnJL6sEn7Zk+HfbtCy8Rkd6Ui8AfMSLZOooV\nvrh98cVk6xCR/Ml04Bd60WkK/AsuCPcEvPBC0pWISN5kOvDT2MMfOjTMq6PAF5HelvnA/8hH4KST\nkq7kg2bOhDVrwMs+N6iISGmZDvx9+9LVuy+YNSvU1tiYdCUikieZDvz9+9N1hU7BrFlhuWZNsnWI\nSL5kPvDT2MM/77xwx60CX0R6U6YDP61DOgMGwPnnw+rVSVciInmS2cA/fhyamtI5pANw0UWhh9/c\nnHQlIpIXmQ38Q4egtTWdPXyASy4Jc/00NCRdiYjkRWYDP403XRW7+OKwfPbZZOsQkfzIbOCncR6d\nYrW14ZGHCnwR6S2ZD/y09vAhDOusWhW+bxARqbTMBn5hSCetPXwIgX/kCKxfn3QlIpIHmQ38/fuh\npgaGDUu6ktI+/vGwXLEi0TJEJCcyHfh1ddAnxS0cNSrchPXUU0lXIiJ5EDsOzewmM9tiZhvN7M5o\nXT8zu8/MXjGzzWZ2a/xSu+fAgRD4aXf55fDcc3D0aNKViEjWxQp8M5sLzAemu/s04DvRps8BA9z9\nXGAG8CUzmxDnXN3V1BSuhEm7K66AY8fgV79KuhIRybq4PfzFwB3u3gzg7tG1MTgw2MxqgJOAFuDd\nmOfqloMH4ZRTevOMPXPRRTBokIZ1RKTy4gb+FGCOma02s2fMbGa0/mfAb4G9wJvAd9z9UMxzdUu1\n9PAHDIBLL4XHH9f8+CJSWZ0GvpmtMLMN7bzmAzXAcGA28L+AZWZmwCzgODAamAh81cwmlTj+IjNr\nMLOGAwcOlKVRra1haoVq6OEDzJ8PO3bASy8lXYmIZFmnge/u89z9nHZejwC7gOUerAFagVrgz4An\n3f1YNMyzCqgvcfwl7l7v7vV1ZfqW9e23Q+hXQw8fQuD36QPLlyddiYhkWdwhnYeBuQBmNgXoDzQR\nhnEujdYPJvwLYEvMc3VZU1NYVksPv64u3ISlwBeRSoob+PcCk8xsA7AUWOjuDnwXGGJmG4EXgB+6\n+8sxz9VlBw+GZbX08AE+8xnYuBG2bk26EhHJqliB7+4t7n59NMRzgbs/Ha0/6u6fc/dp7j7V3f+p\nPOV2TaGHX02B/yd/Ambw058mXYmIZFWK70PtuWob0gEYOxbmzoUf/1hX64hIZWQy8KtxSAfg85+H\n11+H3/wm6UpEJIsyGfhNTdCvHwwZknQl3fPZz4absO67L+lKRCSLMhv4tbVhTLyaDBkC11wDS5dq\nbh0RKb9MBn61TKvQnsWL4d134f77k65ERLImk4FfLdMqtOejH4UZM+Duu/XlrYiUlwI/Zczgy1+G\nTZvg6aeTrkZEsiSTgV/NQzoACxaERzP+4z8mXYmIZEnmAr+1NQR+tfbwAQYOhFtugZUrw0PORUTK\nIXOBX5g4rZp7+ABf+lKYY+cf/iHpSkQkKzIX+IWbrqo98AcPhq99LTwY5bnnkq5GRLIgc4F/KHrM\nSrUHPsCNN4YpF26+GY4fT7oaEal2mQ38k09Oto5yGDwY7rwT1q2DH/4w6WpEpNplLvAPHw7LLAQ+\nhCt2Lr4Ybr0VyvRAMBHJKQV+ypnB974X7r7967/WzVgi0nMK/Cpw7rnhap3ly+GBB5KuRkSqVeYC\n/9ChMOPkgAFJV1JeX/1qGNr5m78JT8YSEemuzAX+4cPZ6t0X9O0LDz4YvsidP//El9MiIl2lwK8i\nY8fCQw9BY2N4Bu577yVdkYhUk0wG/vDhSVdRORdeCD/6ETz7bHhgSnNz0hWJSLXIZOBntYdfcN11\n8P3vw5NPwtVXw5EjSVckItUgc4F/6FD2Ax/gi1+Ee+4JE6zNnQt79iRdkYikXeYCPw89/IK/+At4\n9FHYsgXOPz/MuyMiUkqswDezn5rZ+ui1w8zWF2271cy2mdlWM7s8fqmda2kJX2RmeQy/rSuvhBde\nCPPnX3EF3HRTmDFURKStWIHv7te6+/nufj7wH8ByADObCiwApgFXAN8zs75xi+1MFm+66oqzz4bV\nq8OTsr77XTjrLPjBD8IvQBGRgrIM6ZiZAX8KPBitmg8sdfdmd98ObANmleNcHcnSxGndNWgQ/Nu/\nhd7+xInwV38FZ5wRno37zjtJVyciaVCuMfw5wD53fy16PwZoLNq+K1r3IWa2yMwazKzhQMzZwfLa\nwy82YwY8/zw88US4bv+mm2DUKLjhBvjlL9XrF8mzTgPfzFaY2YZ2XvOLdruOE737bnH3Je5e7+71\ndXV1PTnEfykEfp7G8NtjFsbzn3sO1qyBhQvDl7uXXx6eonXttXD//bBrV9KVikhvqulsB3ef19F2\nM6sBPgPMKFq9GxhX9H5stK6i1MP/IDOYOTO87roLVqwIwf/zn8OyZWGfiRNhzhyYPTtM0nbOOTBs\nWLJ1i0hldBr4XTAP2OLuxf3FR4H/Z2Z3AaOBM4A1ZThXh/I8ht+Zk06Cq64Kr9ZWWL8efv3rcMfu\nL34BP/7xiX1Hj4Zx48KQ0JgxMGIEfOQj4fUHfwBDhoTJ6QYMgP79T/xc/Cqs71vxr+pFpKvKEfgL\naDOc4+4bzWwZsAl4H7jR3Sv+kL5CD1891I716QMXXBBeN98c5thvbIQNG8Jr8+Yw3LNpUxj3j3Mn\nb9++Xf/lUPzerPSrT5+Ot3f11ZmuPnugGvYr7Fv4ub11xQp/PsV/Vm1/Lq6h1PHaO3Zn9Xd3fWfb\nyqE3nkMxc2Z45kUlxQ58d/9CifW3A7fHPX53HD4ceqH9+vXmWaufGYwfH15XXvnh7c3NIfQLr6NH\nw7rCq6Xlg+87Wl9q3eHDH3zfNpiKX62tHW/v6qsrod+Vfaplv7bBXSrMuxLihWVHvxDa+7mr9Xd3\nfWfbyqHSxx84sLLHh/L08FMjT3fZ9qZCz7u2NulKRCSOTE2tcPiwhnNERErJVOC/844CX0SklMwF\n/tChSVchIpJOmQr8t99W4IuIlJKpwFcPX0SktMwEvju8+67G8EVESslM4B89Gq7PVg9fRKR9mQn8\nwhTACnwRkfYp8EVEckKBLyKSE5kJ/MJzXBX4IiLty0zgq4cvItKxzAW+LssUEWlf5gJfPXwRkfZl\nKvD79oVBg5KuREQknTIV+EOHVv4hBSIi1SpzgS8iIu3LTOBrpkwRkY5lJvDVwxcR6ZgCX0QkJzIV\n+LoGX0SktFiBb2Y/NbP10WuHma2P1l9mZmvN7JVoeWl5yi1NDzAXEelYTZwPu/u1hZ/N7J+B6PYn\nmoCr3H2PmZ0DPAWMiXOujjQ3w5EjUFdXqTOIiFS/WIFfYGYG/ClwKYC7v1i0eSNwkpkNcPfmcpyv\nrYMHw7K2thJHFxHJhnKN4c8B9rn7a+1s+yywrlJhD9DUFJYKfBGR0jrt4ZvZCuDUdjbd5u6PRD9f\nBzzYzmenAd8GPtHB8RcBiwDGjx/fhZI/bOBA+NznYPLkHn1cRCQXzN3jHcCsBtgNzHD3XUXrxwJP\nA//d3Vd15Vj19fXe0NAQqx4Rkbwxs7XuXt/ZfuUY0pkHbGkT9sOAx4GvdzXsRUSkssoR+Av48HDO\nl4HJwDeKLtscUYZziYhID8W+Ssfdv9DOum8B34p7bBERKZ/M3GkrIiIdU+CLiOSEAl9EJCcU+CIi\nOaHAFxHJidg3XpWTmR0AdsY4RC1h4ra8yFt7QW3OC7W5e05z906nj0xV4MdlZg1dudssK/LWXlCb\n80JtrgwN6YiI5IQCX0QkJ7IW+EuSLqCX5a29oDbnhdpcAZkawxcRkdKy1sMXEZESMhH4ZnaFmW01\ns21m9vWk66kEMxtnZr8ys01mttHMbo7WDzez/zSz16LlyUnXWk5m1tfMXjSzx6L3mW4vhOnFzexn\nZrbFzDab2YVZbreZfSX6O73BzB40s4FZa6+Z3Wtm+81sQ9G6km00s1ujPNtqZpeXq46qD3wz6wt8\nF/gkMBW4zsymJltVRbwPfNXdpwKzgRujdn4dWOnuZwAro/dZcjOwueh91tsL8C/Ak+5+FjCd0P5M\nttvMxgD/A6h393OAvoQp17PW3h8BV7RZ124bo/+vFwDTos98L8q52Ko+8IFZwDZ3f8PdW4ClwPyE\nayo7d9/r7uuin48QQmAMoa33RbvdB3w6mQrLL3pq2qeAHxStzmx7AcxsKHAJcA+Au7e4+9tku901\nwEnR0/MGAXvIWHvd/VngUJvVpdo4H1jq7s3uvh3YRsi52LIQ+GOAxqL3u6J1mWVmE4A/BFYDI919\nb7TpLWBkQmVVwv8FbgFai9Zlub0AE4EDwA+joawfmNlgMtpud98NfAd4E9gLvOPuvySj7W2jVBsr\nlmlZCPxcMbMhwH8A/9Pd3y3e5uGSq0xcdmVm/w3Y7+5rS+2TpfYWqQEuAP7d3f8Q+C1thjOy1O5o\n3Ho+4RfdaGCwmV1fvE+W2ltKb7UxC4G/GxhX9H5stC5zzKwfIex/4u7Lo9X7zGxUtH0UsD+p+srs\nIuBqM9tBGKa71MweILvtLdgF7HL31dH7nxF+AWS13fOA7e5+wN2PAcuBj5Hd9hYr1caKZVoWAv8F\n4Awzm2hm/QlfdjyacE1lZ2ZGGNfd7O53FW16FFgY/bwQeKS3a6sEd7/V3ce6+wTCf9On3f16Mtre\nAnd/C2g0szOjVX8MbCK77X4TmG1mg6K/439M+H4qq+0tVqqNjwILzGyAmU0EzgDWlOWM7l71L+BK\n4FXgdeC2pOupUBsvJvyT72VgffS6EjiF8A3/a8AKYHjStVag7R8HHot+zkN7zwcaov/WDwMnZ7nd\nwN8BW4ANwP3AgKy1F3iQ8B3FMcK/4r7YURuB26I82wp8slx16E5bEZGcyMKQjoiIdIECX0QkJxT4\nIiI5ocAXEckJBb6ISE4o8EVEckKBLyKSEwp8EZGc+P+OFr1w1stCygAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1196bcf28>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# create timescale. t_interval is the time interval in which to calculate the solution.\n",
    "# t_points are the points at which the solution is stored.\n",
    "t_interval = (0.0, 100.0)\n",
    "numpoints = 1000\n",
    "t_points = np.linspace(t_interval[0], t_interval[1], numpoints)\n",
    "\n",
    "# solve coupled ODEs with scipy's solver\n",
    "soln = sp_int.solve_ivp(f, t_interval, y_0, 'RK45', t_points)\n",
    "\n",
    "# plot & show solution\n",
    "plt.plot(soln.t, soln.y[0, :], 'b')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:anaconda3]",
   "language": "python",
   "name": "conda-env-anaconda3-py"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}