Hodgkin-Huxley.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hodgkin Huxley Model"
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAATIAAABoCAYAAABypHzgAAAStklEQVR4Ae2defB+1RzH35aISJQsiWRJUYRQliwVkRRKE2WXLNnGTtnXrBGRJbJMtorsskRFaRDJFmOQbWyDGf5hXn3Pmd/5nc79Ps/d7/M878/MnXufc+49n899n3vfz7ln+XwkixEwAkbACBgBI2AEjIARMAJGwAgYASNgBIzAAiNw+gLbbtONgBEwArqCpDONgxEwAkZgkREwkS1y7dl2I2AELkPAROYHwQgYgYVHwES28FXoGzACRsBE5mfACBiBhUfARLbwVegbMAKrh8CDs1s2kWWA+KcRMALTRqBEWqW0ad+FrTMCRmClELiRpF2SOy6RViktucSHRsAIGIFxEXi+pCMTE0qkVUpLLvGhETACRmBcBE6VdIfEhBJpldKSS3xoBIyAERgXgUskbZKYUCKtUlpyiQ+NgBEwAuMhQP/Y2Zn6EmmV0rLL/HNoBKi8gyUdEbZDMwPumORxzqZJ/nZZ3rZJ3mZZ3q5JHoeHJ/n7ZHn7Jnm2R+oKn+tLummGtX9uQOAhkt684edlRyXSKqVll/nnkAjcWdJ/JR0l6SBJh0jaPzNg55BOHttVkvxtsjxelChXy/J2ihlhf0CSf7csb88kz/ZIXeGztaRPStorw9s/1xC4e9Y/dnVJx0s6R9LRCUgmsgSMsQ+vIekESduPbYj1D4oAZPY/SbcdVOviKjtR0sVzjGQu7h3aciOwoAjsVmh5L+it9G72UyRdKmn3RJNbZAkYPjQCRmD6CPC5+S9J9PlGMZFFJLw3AkZgIRC4qqTzMktNZBkgY/68lfvHxoTfuhcIgbQ1htkmsglVHiNh6QzmCZlmU3pGgAGeZ/asY5mLN5FNqHZNZBOqjIFNuY2k4wbWuUzqTGQTqk3mkHkIfkIVMqApJrJ2YJvI2uHnq41AJwjcXNLzOilpNQsxka1mvfuujcBSIWAiW6rq9M0YgdVEwEQ2oXpnsXi6NnJCptkUIzBpBExkE6oej1pOqDIGNuUGwUnAwGqXRp2JbEJVWUVkeEZ49px2npKdt5+kw7I0/5weAlMbtWQ+40MXbLtwwewF35sM8Cim3nEGULfmGqY0IfZASW+QdD1JO4bKunWwiLSHSYouee4d0q8piZBZLLB9eki7haQ9BrkTK6mLwNSI7Jfh2eH5GXKDjJrqY0Jx02vb6G2q86WSvlX3Qal5Pm6PTqp5TevTYc4Se0Yiw8ncbyQ9Q9IvwqLZb0h6tKTvhv618yVhPG5OGM7/USCyp4UbeoukV7W2dLkL2ErSnZJbxKlkdH74wCS9y8OpEdlHu7y5GmV9rca5XZ46hl4aIX3jDGfACSx/HF1SIoveMk+WxMP/WUkcYzD9BBjNi/fqYPVzApFBdK+X9ApJ3x/9jqZtAN4V4r8YeJ0mCXKDxHjwcu+wXd3NFbsqqINy+n7Bqkwcg1CwZQy9QxDZpPoOUyKLkybfF/wx0Qd2e0lflrR3IDJWB5CPQF58Wn5HEi6y+bwkvJalGoFIZC+X9LEQ+IJWLvi9e0VWX5jIqp+PrnKWlsi2kLR5ASVcUB8bWl2RyN4rCZfXuPzFVTKfmNcK7k34Z4fYviDpokBk9JedK+ksSU8q6HDSBgQgsr9L+n1htv29Npy21Ecmsv6rd2mJrGrUchakW1accJ0snQEAyM6yPgIQ2ZmSri3p56Eltv4V7XP5E2Ot7VTERNZ/TZjI+sd4pTXET0tAoAXGwAnxFPoUd/avoTtGXxWax9BrIuvzjXLZurEkPuejPCIbxYzpXe5NZGtojkEoaB5D79ISGZ+CfGJYVg8BE9lanY9BKGgeQ+/SEtnqvb6+44gAMUeHmOUd9c3aD9lHxohwXEXww3B8v1kGdpAP5rlepi71LVHvYyV9M9jQl95JTb/oG1iXbwRyBIYksheHkXdmyMfZ+cx7jKtUctu6+v1CSZ8KqwFSvfcJCq4kiSk4zAjAX1xXkuplcjv3zf32oXcUIruyJMCzGIGxERiSyJhSdN/shktp2Smtfz5XUt7yi2lMbfqMpD9JOrS1po0LiDrS1JjWtd5RiKzp9IsUEB8vJgKEOLvhhExfZSJ7VJiHeUGYAN11tUTSSsslrQ+9JrIUZR/3jsAqd/aXWl+ltK4roYpQjpF0RiAzfAR2LUPqXQoiYwE6jhrTjX/+9eS62flMCu1ahtARbW6CQbx2yL2JbGO0xyKyt0k6WxJzCRlV7ENKRNaX3lGIjFGrLj8v6KAEoHSL7n+qKoiO1/T8J1Sd2CJ9CB3RvCYYxGuH3JvINkZ7DCJjZQWjpowm9ik5kfWpdxQi6xM8lz1tBJhDSCtgKrJqfWSQyanByUI+ANB1naRE1rdeE1nXtddBebRY+5pv04F5S1XEKhFZJJMdJKUk01eFRh1D6DWR9VWLDcuFwHCps03D631ZPQRWhchwOEpLDBJDIsmEn73s0DGU3lGIbLfgmqcX9Ba4UJr6JrH+KnBXSQ/KNiaB5mnb9mRCqT+slNa1+tjBHkmM8ocgsiH1jkJknkd2+Ud1VUiMmAvR++/lUWiWwmg1Hdf4nyNKU5XgkBPnm+lGsJv0N8cxLkRVOVXpuI86OGx4Ms5nyZdIq5RWVX6TdD7rSh37OZHh/STaTh/mJk2UJdfMqzde0la/iSwiOeJ+VUgMiLsetbxdeFFZDvOi4FdtrKrEizGOPp8cliD9UbrMw0i0p0RapbR4ftt97JvCe3LesZ8T2SXBbmz/oiQChjSVOnqjjrb6RyEypkbk/1bxhtL94yV9OyyhYBnFVDYCovxZ0t9CdBgWxP41pBGVp46dX5X0PUmfD2Ux94yRPT556pTT9blEvcGWrqVrIuOlwwV6FF7EWXMI47ld71mvyKx1BBsIiJPO0SqRViktFNFqF8mkqmM/JTImwzKnLArThshvInX0xvK70D8KkcUbmLUnNBxRfnD2N8XtjZIIa8fL+fYObOy6vDaYRVtm1VHdfAYxulrTR2wBgtJMRXC5TswI9vzZvSQzrERapbTssto/UzLh4pS0YmFpGp/Bvwt/phAaBNzEwWZdvdGWLvTnRPa6MOGdP5TRo6kR0o0QcNEFydT25yVE9oEO7Oy6vDZ4RVviwzbF/T0kfSQxjBdprDgNvEi/CiEKcddDVKrXJLZxWCKtUlp2Wa2fOZlwcUpasbA07bWSHikJtzuEUyOGw0GSjgsn049IP+R60kRvLK+k/2XBXT0DMfN85uZExlcOcpdQF+HnOLvtwvBt08CgfV9HywIvHgRBIXhwW31dl9fGnmjLODU/n1a6J34q6Z4hhB3kG+NyzldCd2fdUtJXkuJwW4PLnFRKpFVKS6+pewyJ5uH2UtKK5aVpdG/cLGRsJunH4RMZJ4x0/L9SEp9/2PrBZApHLIt9E73x+pJ+RjxpxLyzcD/xunSfEhmkS1cNQphISLkXgfW376VkF7pqCOwp6cTw+dDXtIl5MIXI7p+ciI+x3M9YibRKaUkxnRympBULTNMYqY1C0B5efvr6GCD4eFgFAJEdJemEGv1nqY5YfimtpB8i+4Sk98QLZ+xTIoN8ITJGqQn0Td30Ip5+0QusC1Eof2B5S2Uow3mJaGWkG0Gc098c5yN8XdlXIq1SWlf6Yjkl8iilxfPZQ2RMieCPAl9ld5V0SphzN2+82JKOUlqqNx5DZAx68Yl5SExcZ58SGa1DAnfTEqNVR14v0heRsbSHKNl15Kl1Tm55bhP76qocQkddm9Lzux61TMtucrwqM/tzcp5FKExroWsHeUyIQwvpcjyvg4WSjlJaULPRjtYtfXY4YT1yo5zyj5TIymf0kErHIFHCuxaG4V9Qs1Car0NJE/vq2jaEjro2peevMpEx1w0SSAdhaAHmXmNTvLo4pgVV0puTWxe60jKG1DsKkaU329XxSSHYLEPHkNOsPriHh/N+G/ZP7MqQinLq2ldRzLrJQ+hY14A5Mumg5599KjJki2yrZPDqB2Fw6HEDuH3nE41+KAZ/Lkz00trpU4bUuzRERoU0aY24Rdbnozz9socksojGmC8dLcAxpG+9Y2LaOZ4Mu9b1dcX6sqGkiX11bRtCR12bpny+iWyY2llKImMolxfOYgTGRsBENkwNLCWR9TVqOUyVWEsbBPBO0dsExQaGmcgagNbgEhNZA9B8yXQRWOVRy1grY/bnzCIUgudGdz7sqa8uZJZeRnJTvXVjeoyCqVtkXTwai1nGMhHZppKYUHtOsuHkYJaM8tIFo2YRCt5bWF4VNwbQupD19LIi4yeJTnTX9ZA8CqaELmOzrB4Cy0RkLEhmtvssYT4VL3K6sUY0/c3xPrMKqplfVy8TYHEf1VaY9JrfW+l+4zw2WmNvaql0FCJrabMvHwmBJq5dSqbmi5pL5wyV1qaPDH/0pwXvtHioxYX7IguODy5NWpe4whlC0POzRC+rB+qKiawuYkt0Pp9GFyWOCfFpj/uZKmHd27JJGyL7kCScKTKplY31iDuvAxBrAPP4AHutc35XWfPqxf9fPh2Jibs42ORTr+4aWZYYHZBtpVYnrbd8Avtbw/Kkd0niuZwloxDZFmHt1izjnN8vAjxoeLalj4dW0h6Jt4EdJR0W1tvR8QqJ8dmBb6plkjZEdrEk3N5EgcT4nKuSzZP4AMywZ7Y9i5r7llQv/vur9J4VCJklU2yswsBTMK50WMPLM8Ai8uPD6oBZdhPDgFYrG3rZ5/eLGywWpfOZGfVuHYLw4OMv9yBSpXMUIluvsx/25wUC7HncYec3xoTYtP8NMOfpgM3LWYXfENkFkvgnZvlKJDJ8U+Fb6/BAcjwkSBctMv7EWGs7FWlKZCy/+U/yOcSfAZ9m6xFZvOdRXrqgnNZPlfAnhVfbuPHeQCb0eeENg/pnMTfOK8+vKqQivUovboOivrin8x9yO6PGBPdRMK0iMm6ATkEY/2hJvw6uRErY4G+o5JudJn7al4PvM5zBWS6PQCQyWhWMvtHpiv8nPkX4jHiHpH9kfwyXL6VeyjJ19ud3PqtFFs/npaNFNoZUEUqVLRAZEZVwfgiRfTpMyfhc1QUV6XX1fji8x2dK4s9vlkyKyIj2kvbRwP7EeYT4YgfgMZIODIE+8CbANzjMDbER6IHWBTfOv+PXw6gSRMY/6LEhys3+s1BZkfxIZNzuA4KfeYgMx3pgxSJnhsXTFm5baJadyM4N8Rvi81rCi5fuL6WMAdLqEgrvF4K3muifDEI7QlKduV519eLHH6F/bJ5lh6MQGaRSYtnTsxYYLxqjKPgYo4WGQEq0sohcRIvsS4HoIDvAhtSYPU4QiC0D4FzzLElEZeJTE4KzrPl74tMyClMJIDI+/cCPETm8icSHKp7XZr/MRDYvLotEZPPe06zz6hLZrPLy/FGILDci/s7DUMHEfPKkRMbLlhIZ5HVyaKoTjAIiI9wcYcIQyoDIuI6IOzSV+Q3JWaoRoIM4xiPoskXGnxMdx1ORpn1kbexfNSKj35soTX3KpIiMkTLCaOG1lf4aIqHQ8Q8ZQVA0bekEhMhgeM5nz8REPjmJach5jHoQiIKh7fcH4mK2MBGu6cAkXBdTDyxGwETW/zNAl88felYzCpHh0I1/+5LsEvq5mEcSPWdyLv02hNeC2Oi7YViWEFa7B4JiaJd/e6L/0FlNq4zYjPRV0KfG9AJc5uILnFiUFiMAAiay/p+DpSWyqlHL/iG1hrERoF+zTidx3/aayPpGeG0Qru8WGQ0gggsP2m1hIuv/4Zmqhql19rM0hm6Mobd/jqCTe2Ri6tD3ylcSXUJ9CsFScFs/aJeRiazPKp122VMjsv3CZOA2wYybXEvcgibXcU30u9/keuYHNrmu7TVMtu5bBu8yovk3pc+LvgF2+RsQmBqRbbBscY76nsqwOEjYUiMwEgLMIZxnguNI5i2EWhPZQlSTjTQCRmA9BExk66HjPCNgBBYCARPZRKoJB3Tr+W2aiJk2wwhMEgET2USqBQd08zhLm4i5NqNDBHYKE5g7LHLlijKRTaTKCWgQV9VPxCSbMRAC1Dtrai3NETCRNceu0ytZSsRyIZYqWVYHARaf/9ut8dYVbiJrDWE/BTCnLJ1tzNqsKJBemrdDzAh7JjTGfBaPp8LE25i3b5oR1mLGPPaphwfbs7ZWtWt88MuO3zPL/AjQBZP7+WdCbJ6GU1LLyAjwgKczh9NKwdtrmpfP3iV6TcyP4aXi7eD+J+bhtjkVfG7FPPbp0gbbs+aTrGt8qhwFpPXi440R2Dvx84/7dzYcKMTjuKfv0WIEjIARMAJGwAgYASPQJwJ0t9A3lm5jeO7o8x5dthEwAkbACBgBI2AEpo0Ag155Bz/9ZxYjYASMwMIgQNCc2LEf9+tFalqYG7OhRsAIGAEjYASMgBEwAkbACBgBI2AEjIARMAJGwAgYASNgBIyAEWiEwP8B5U1yCmZsBxcAAAAASUVORK5CYII="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As with the leaky integrate and fire model, the Hodgkin-Huxley model of a neuron is realized by considering a current balance law, treating the neural membrane as an electric circuit. Some models consider the fact that the neuron is spatially extended, which requires some knowledge of advanced linear algebra and/or partial differential equations (see most of Ch.3 of Gerstner et al (2014)). We simply consider a local model here treating the neuron as a single piece of membrane, so the applied current $I(t)$ is split into capacitive current $I_C(t)$ and the sum of all ion channel currents $I_k(t)$ ($k=1,2,....$), so\n",
    "$$\n",
    "I(t) = I_c(t) + \\sum_k I_k(t).\n",
    "$$\n",
    "The Hodgkin-Huxley model has three ion channels as shown here:\n",
    "<img src=\"attachment:image.png\" width=\"300\">\n",
    "a sodium channel (Na), potassium channel (K), and leak channel (L) sometime considered a chloride channel. Noting the capacitance formula $I_C = C du/dt$, we have\n",
    "$$\n",
    "C \\frac{du}{dt} = - I_{Na}(t) - I_K(t) - I_L(t) + I(t)\n",
    "$$\n",
    "It is important to note that each of these ion channel terms is really the average of many ion channels. As in the LIF model, the leak channel's current is proportional to the difference between the neuron voltage $u$ and the resting potential $E_L$ so Ohm's law ($V=IR$) implies $u-E_L = I_L R$ or $I_L = g_L (u-E_L)$ where $g_L = 1/R$ is the *conductance* of the channel.\n",
    "\n",
    "The other two channels follow the same physical laws, but have a probability of being open between zero and one that depends on the voltage $u$, and this probability evolves according to dynamic *gating variables* $m$, $n$, and $h$. The sodium channels is controlled by $m$ and $h$, so $1/R_{Na} = g_{Na} m^3 h$, so $m$ describes the channels opening and $h$ controls channel inactivation (blocking). The potassium channel conductance is $1/R_K = g_K n^4$, so the sum of current is\n",
    "$$\n",
    "\\sum_{k} I_k = g_{Na} m^3 h(u - E_{Na}) + g_K n^4 (u - E_K) + g_L (u - E_L),\n",
    "$$\n",
    "where the variable prefactors ($m,h,n$) determine the fraction of ion channels that are open. The gating variables evolve according to *voltage-gated* transition rates, meaning the rate of transition between open $\\alpha_x$ and closed $\\beta_x$ changes with the voltage. Putting these equations together with the voltage differential equation yields the system\n",
    "$$\n",
    "C \\frac{du}{dt} = I(t) - g_{Na} m^3 h(u - E_{Na}) - g_K n^4 (u - E_K) - g_L (u - E_L)\n",
    "$$\n",
    "$$\n",
    "\\frac{dm}{dt} = \\alpha_m (u) (1-m) - \\beta_m (u) m\n",
    "$$\n",
    "$$\n",
    "\\frac{dn}{dt} = \\alpha_n (u) (1-n) - \\beta_n (u) n\n",
    "$$\n",
    "$$\n",
    "\\frac{dh}{dt} = \\alpha_h(u) (1-h) - \\beta_h (u) h.\n",
    "$$\n",
    "This fully models the voltage dynamics of a neuron up to the needed parameters. We build code using classic parameters (rather than those in Gerstner et al (2014)) from Hodgkin-Huxley shifted so that the resting potential of the neuron is close to zero. First, we present the resting potential $E_x$ and conductance $g_x$ parameters:\n",
    "$E_{Na} = 115$mV, $E_{K} = -12$mV, $E_L = 10.6$mV \n",
    "$g_{Na} = 120$mS/cm$^2$, $g_K = 36$mS/cm$^2$, $g_L = 0.3$mS/cm$^2$,\n",
    "where the conductance is given in millisiemens per area (square centimeters). Membrane conductance is assumed to be $C = 1\\mu$F/cm$^2$. Note, in the code we will take ratios of conductance and capacitance, which will yield mS/$\\mu$F$=$1/ms (1/units of time). The form of the voltage-depedent gating variable switching rates is somewhat complicated and will be given in the python code below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, to walk through the Hodgkin-Huxley simulation code, we start by importing plotting and numerical computation tools along with the ode integration package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from scipy.integrate import odeint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Subsequently, we define the conductance, voltage reversal potential, and capacitance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "gK = 36.0   # average potassium channel conductance per unit area (mS/cm^2)\n",
    "gNa = 120.0 # average sodium channel conductance per unit area (mS/cm^2)\n",
    "gL = 0.3    # average leak channel conductance per unit area (mS/cm^2)\n",
    "\n",
    "EK = -12.0  # potassium potential (mV)\n",
    "ENa = 115.0 # Sodium potential (mV)\n",
    "EL = 10.6   # leak potential (mV)\n",
    "\n",
    "Cm = 1.0    # membrane capacitance per unit area (uF/cm^2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then the strength of the constant external current:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "Id = 10    # membrane capacitance per unit area (uF/cm^2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then the end time and time vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 100.0    # end time (in milliseconds)\n",
    "tvec = np.linspace(0, T, 10000)    # time vector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, the voltage-gated ion channel rate functions. Note, the *def* syntax gives functions which must return something:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# potassium ion-channel rate functions\n",
    "def alpha_n(Vm):\n",
    "    return (0.1-0.01*Vm)/(np.exp(1-0.1*Vm)-1)\n",
    "\n",
    "def beta_n(Vm):\n",
    "    return 0.125*np.exp(-Vm/80)\n",
    "\n",
    "# sodium ion-channel rate functions\n",
    "def alpha_m(Vm):\n",
    "    return (2.5-0.1*Vm)/(np.exp(2.5-0.1*Vm)-1)\n",
    "\n",
    "def beta_m(Vm):\n",
    "    return 4*np.exp(-Vm/18)\n",
    "\n",
    "# chloride ion-channel rate functions\n",
    "def alpha_h(Vm):\n",
    "    return 0.07*np.exp(-Vm/20)\n",
    "\n",
    "def beta_h(Vm):\n",
    "    return 1/(np.exp(3-0.1*Vm)+1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What motivates the above forms of the voltage-dependent rates? Consider for instance both $\\beta_n$ and $\\beta_m$, which are both exponentially decaying in the voltage ($V_m$ in the code). This means that when the voltage is low (when the neuron is quiet) the rate of off transitions is higher, which keeps the neuron from spiking spontaneously. Thus, an external current is needed to move the neuron from equilibrium."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also compute the steady state values of the gating variables by solving for the equilibria of the gating variable equations when $V_m = 0$mV (about its resting potential)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# n, m, and h steady-state values\n",
    "def n_inf(Vm=0.0):\n",
    "    return alpha_n(Vm) / (alpha_n(Vm) + beta_n(Vm))\n",
    "\n",
    "def m_inf(Vm=0.0):\n",
    "    return alpha_m(Vm) / (alpha_m(Vm) + beta_m(Vm))\n",
    "\n",
    "def h_inf(Vm=0.0):\n",
    "    return alpha_h(Vm) / (alpha_h(Vm) + beta_h(Vm))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These will be useful in setting initial conditions. Lastly, we define the right hand sides of the differential equations as above, again as a function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute derivatives\n",
    "def compute_derivatives(y, t0):\n",
    "    dy = np.zeros((4,))\n",
    "    \n",
    "    Vm = y[0]\n",
    "    n = y[1]\n",
    "    m = y[2]\n",
    "    h = y[3]\n",
    "    \n",
    "    # dVm/dt\n",
    "    GK = (gK/Cm)*np.power(n,4.0)\n",
    "    GNa = (gNa/Cm)*np.power(m,3.0)*h\n",
    "    GL = gL/Cm\n",
    "    \n",
    "    dy[0] = (Id/Cm)-(GK*(Vm-EK))-(GNa*(Vm-ENa))-(GL*(Vm-EL))\n",
    "    \n",
    "    # dn/dt\n",
    "    dy[1] = (alpha_n(Vm)*(1-n))-(beta_n(Vm)*n)\n",
    "    \n",
    "    # dm/dt\n",
    "    dy[2] = (alpha_m(Vm)*(1-m))-(beta_m(Vm)*m)\n",
    "    \n",
    "    # dh/dt\n",
    "    dy[3] = (alpha_h(Vm)*(1-h))-(beta_h(Vm)*h)\n",
    "    \n",
    "    return dy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to set initial conditions as an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# state (Vm, n, m, h)\n",
    "Y = np.array([0.0, n_inf(), m_inf(), h_inf()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And solve the ODE system using *odeint*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# solve ODE system\n",
    "# vy = (Vm[t0:tmax], n[t0:tmax], m[t0:tmax], h[t0:tmax])\n",
    "Vy = odeint(compute_derivatives, Y, tvec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally plot the neuron potential as a function of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2deXhdVbm43y8nU9N5prS0pQMzlI4goCigF1EBuSAgQ1vAouIFuU7o/TlduYqKA07IUCgggoooIAIKgopDR6C0FNrS0tLSCTqlTTOe9fvjnIQ0e+199k7O8K1kvc/Tpzk7Jzlf9l5rfeP6lhhj8Hg8Ho8HoKzUAng8Ho9HD14peDwej6cNrxQ8Ho/H04ZXCh6Px+NpwysFj8fj8bRRXmoBusKQIUPM2LFjSy2Gx+PxOMXixYvfNMYMtX3PaaUwduxYFi1aVGoxPB6PxylEZF3Y93z4yOPxeDxteKXg8Xg8nja8UvB4PB5PG14peDwej6cNrxQ8Ho/H04ZXCu0wxrB6ay2765tKLYrH4/GUBKdLUvPJ9r2NXHHXQpas3wnA//vA4VzxznEllsrj8XiKi/cUsvzs6dVtCgHg+kdXsHJLbQkl8ng8nuLjlUKWFzfuClx7asXWEkji8Xg8pcMrhSzz124PXFv0WvCax+PxdGe8UoigxZ9K5/F4ehheKUTQ3OKVgsfj6Vl4pRDB3sbmUovg8Xg8RcUrhQj21Hul4PF4ehYFUwoicoeIbBWRZe2uDRKRP4vIquz/A7PXRUR+JCKrRWSpiEwplFxJqG9uKbUIHo/HU1QK6SnMA07vcO064CljzETgqexrgPcDE7P/5gA3F1Cu2DQ0pUstgsfj8RSVgikFY8zfgI41nWcBd2W/vgs4u931u02GfwMDRGREoWSLS32T9xQ8Hk/Potg5heHGmE0A2f+HZa+PBF5v974N2WsBRGSOiCwSkUXbtm0rqLD13lPweDw9DC2JZrFcs9aDGmNuNcZMM8ZMGzrUesRo3mhsSWP8XgWPx9ODKLZS2NIaFsr+39pHYgNwULv3jQLeKLJsVhqavbfg8Xh6DsVWCg8DM7NfzwQeanf90mwV0vHArtYwUzFIp8O9gbpGn1fweDw9h4K1zhaR+4B3A0NEZAPwVeAG4NcicjmwHjgv+/Y/AmcAq4E6YHah5LKRjggR1TU2M6h3ZRGlSU5zS5q0gcpyLdFAj8fjKgVTCsaYC0O+darlvQa4qlCy5CKqx5F2T+FPyzfz9UdeYmttPaccNoybLphMdUWq1GJFUtfYzJ9f2sKQPlUcd/AgylP6ldkrm2tZtbWW6WMHMbxfdanFyUlzS5p/rXmLilQZU0YPdMJg2LCjjmUbd3P4iL6MGdy71OLkxBjDwtd20JI2TB49QP28i4s/ZAdIR6QNNCuFXXVNfP63S9lZlzkp7onlW/juE6/w5Q8eUWLJwnlrTwOn3/R3ttU2APCBY0Zw0/nHqlYMd//rNb768HKMgcpUGfNmT+eECUNKLVYo9U0tnH/Lv3hhQ6Yd/IkTBjN35nTVi9ZTK7bwiXuX0NicJlUmfO+8SZw92VqAqAJjDJfNW8jTr2QqII8Z1Z/75xxPTaX7S6remVhEIj2FBr2tLhat296mEFr528rClul2lZufebVNIQA8unQT/3z1rRJKFM2+xha+/djLtA6RxpY0X3/kpdIKlYNfzl/fphAA/rH6LR5+QUXdRihff+QlGrNFHS1pw9ceWU5LRK6v1DyxfEubQgBYumEXc/++toQS5Q+vFIjOKexV7Cm8vr0ucG3V1j0lkCQ+Sy2HGf1VsSLbuLMuMAZe2VJLc4veqrRlbwTv8XPrd5RAknjU1jexvsNY3lnXxObd9SWSKDcvbtwZuPbnFVtKIEn+8UqB6OojzbuaK8vt4YCov6fULLAcZrRQ8WFGYeHDN/c0FlmS+DyxbHPg2n0LXre8Uwf7wu5xO49SG8+uDnq3SzcElbGLeKUAkW6qZqXQFGKt1imW2YbmcyvClEJtfZP1ugZce/77QuStVdyluLzMtt+2e+CVAhBlWNcr3rwWNpn2Ks6D2IgK35WaMCt2j+J77NqCFTaO9zToVbwpceseJ8ErBaIXpQbFVlfYgqXZwrLRrDjcFa549Y6LMscWrLAeY7sVj2PHbnEivFIgOnwUtvBqoDEkfOSap6C5yiQsfKjZU6hyYE9Ce1w0bspT3VcruDV6CkRkTkHxQTthCWXNC5YNzd5YmBejWfFq3o9gIyw3tk/xcbipMvvS2aB4vYiLVwpAVEhbc/vssLCXZu/GRliIRgNhilfz+d0u7F5uT9g+Ic3jIqx7suawYlzcGj0FImrzmubqo7BSec3ejY1Gxcn8MCdS8+SvDNkdrrVUOWyB3deod1yERRc0e5Bx8UqBXCWpegdmmKegWWYbTYpLUsMMBs3GQtjd1GosuGjchM09Hz7qJkRVH7k4MDW73TbCEuYaCLOuNY+LMCNHq7EQKq/iMGhYvzTN3k1cvFIgh1JQPDDDJpPmxG1YKZ/WthEuLlhhMms1FkLDR0rlhQgPUrGxEBevFIjukqr5IbuWaDbGhCb1tW4SdDFE59q4cDHRHD4u9MocF68UyOEpaJ78jsViI3eOK51MoZNf6T2GqPCRTpnDxoVWJQYRYUXF60VcvFIgV0mq3oEZngTVOTBd3CQYmgRVPC5CF1mlMofnbXSOY3DTu4mLVwqACa3XUD75HYsdR7YTUWp5+/BR4XExbxPqpSude0nwSoFcYQ33Jr/WgelimM41xQsuho/cu8ehJamKZY6LVwq4acEChJX3a538keEjrTI7OPldU2TdSSloNW6S4JUC0TkFrS43uJfsiqzyUroAuHaPwb0Nd2G2glZ5wb2y3yR4pUB4nTToTna5Fjt2MXzkYj26a+PCtXAXuKnI4uKVAuFtASAzYMO6OJaa0MmkdMGK6jGl1cJys1zSfn2fUsUbpsSaWvTOPR8+6ubkahSmVfu7NjBdPAs7PHykU15wMHzk4LhwzSBLglcKRFcfgeJF1jEX1sXNa+GTX+eYAPfCMVH9ENXOPcfucRJKohRE5FoRWS4iy0TkPhGpFpGDRWS+iKwSkV+JSGWx5InapwB6H7R7k99BizBE5sbmtMpW1C5a3ZE5PaUyu2aQJaHoSkFERgJXA9OMMUcBKeAC4NvAD4wxE4EdwOXFkinXufFaH7RrpXzRC5ZOizBqbDQo9BacTOZHjAutJeGudRNIQqnCR+VALxEpB2qATcApwAPZ798FnF0sYaImEuh90K5tXnNyn4JjMkd6Y44tsKB47jnmpSeh6ErBGLMRuBFYT0YZ7AIWAzuNMa3HFm0ARhZPpujvq51METX0US55qYi2YpXeY8dkdnEvSOQ+IaUyu2aQJaEU4aOBwFnAwcCBQG/g/Za3Wu+6iMwRkUUismjbtm15kSm3p6DzQUctAD60kR9ci9G7aHVHn3yo7x5D1OY1nfc4CaUIH50GrDXGbDPGNAEPAicAA7LhJIBRwBu2HzbG3GqMmWaMmTZ06NC8CJTLptY6mSLbcyiUOarkXGvbCNcUmWvhLnDvHkO4d6N1HCehFEphPXC8iNSIiACnAi8BTwPnZt8zE3ioWALlCrWotVYcix+7eOxplCLTKHPUWNa6YLnmjYF7e0GSUIqcwnwyCeUlwItZGW4FvgD8t4isBgYDc4slU1QYBhRbWI5NpugwgU6L0LU8iIv32LW8Dbi5fyUu5bnfkn+MMV8Fvtrh8hpgRgnEyRk+UmthObbpx7UFFhwM0TnojUWOY6WLbOixskrHcRL8jmbcLUl1LX7cnXY0g06ZXaw+ivR4lfaYigofaaz8S4JXCribU3DN8nYxtOFaHiSy6aDSBda1cQzhYzltMo38XMYrBXLvU9BodYN7k8m1BRbcU2SRVrfSUEx3SuaDTpmT4JUC7jbEc23BimxnoFBecC/kFaV41fZrcqwkNWocg85xkQSvFIjREE+p5o8amxp7xrjm2YB7/ZpyLVjubWrUNy6iQnSg18CJi1cKxPAUfCw2L7iYBHWtXNLF3fmuebyuNtCMi1cKuBsjdG0yRZdL6uzXFHmPFY6LXAeVaZTZtRBd7vCRvrmXBK8UiKP5dT5k1zavRVmxmWNP9SmFyNbZCsdFLk9BYwWSa+M4V/hIo+JNglcKuDmRwMHNa7ksLIWTybV9Ci5asa5tuDO5vDGF4yIJXikQI6egcGCCe5PJxaqNbpdTUDgunKs+cnSza1y8UiDO5jWdD9k1tzvXgqUyHONa3sZBxdvdxrFGmZPglQLuVhO4ZmHl3g+i7z67tuHOScXrWKI5ZxhUocxJ8EqBGPsUlD5kH+8uPFG5b433OJfi1bg737WDgXInmvXJnASvFHBzIoGDFpaL8W7HwkcuWrHuhY+iv6+1q3JcvFLA3Rihc5vXHLzPrt1jF5Ogrt1jFxVvErxSIN4+Bec2Vmmc/DlL+RTKHNlKRJ+8LubHohvi6bvHLoZBk+CVArmrj8C9BUBjKCa3FatPZvdCG+6F6KLmX2ZTo66556LHmwSvFMgdIwSdD9q16iMXz61wLbThohXrmrHgouJNglcKxPMUXJtMGpNdufvy6LvH7ine6O9rW2DBPUWWcxwrkzcpXikQz1PQVoFkjIkOHymTF+LU0GuUOfx79c36jl50MbThmiJz8R4noTzqmyIyCrgAeCdwILAPWAY8CjxmTK4uIG4QZ1pre9A5J5JjVjfou8cQLbMxmVxTdUWqiBJF4+KOZtdkds2zSUqopyAidwJ3AI3At4ELgU8CTwKnA8+KyLuKIWShiRc+UjYwXVxgHZxMOQ+tUSazi6cI5jYWdMmca7nQeMBVEqI8he8ZY5ZZri8DHhSRSmB0YcQqLrkGJegLH8Wxuo0xiEiRJMqNi20ucntkLfSnojjCxKC7JfNBX+LWRYMsCVE5hdOz4SMrxphGY8zqAshUdOKEhdVZhDnESRvUnU/g2uQH90JeuRYsbcYNuGcsuObZJCVKKYwE/iUifxORT4jIkGIJVWxcLEmN491oW2RdnEyuxY9zt2DQJS84eI8dy4EkJVQpGGOuJRMe+jJwDLBURB4TkUtFpG+xBCwG3TF8BPrOls5dyqdLXnDPinXxaNmch1wpu8eunr8Sl8iSVJPhr8aYTwAHAT8ErgW2dOVDRWSAiDwgIi+LyAoReYeIDBKRP4vIquz/A7vyGflGn7WS+z3qZHbQU3AtRu9aJQ84GKJzzLNJSqx9CiJyNPC/wE/JVCN9qYufexPwuDHmMGASsAK4DnjKGDMReCr7uijkcgdB38B0MXyUa4HVWLXhWptkF6uPcnmQ2vav5DQUlHnoSQmtPhKRiWT2KFwItAD3A+8zxqzpygeKSD/gXcAsyCSsgUYROQt4d/ZtdwHPAF/oymfFxcXNa7kWK9CnyJwMHzlmebtmdYN7Jx/mNhT03eMkRJWkPgHcB5xvjHkxj585DtgG3Ckik4DFwDXAcGPMJgBjzCYRGWb7YRGZA8wBGD06PxWxuQ7ZAX3WSixPQdlkcjN8FP19bYtsLiWmzbgB90JeuYzIphZDS9qQKtNTDp6EqETzOGPM/7QqBBHpl437DxKRQV34zHJgCnCzMWYysJcEoSJjzK3GmGnGmGlDhw7tghhvE6v6SFmYIE4Zrb7J5OCC5djxli5WH7lWquxiuDkJOXMKInKliGwBlpKx6hcDi7rwmRuADcaY+dnXD5BREltEZET2M0cAW7vwGcmIU32kLE6Yy7oCfQPTxdCGawtWLiXW2JKONXaKSc6TDxt1KbJ4XrqucZGEOInmzwJHGmPGGmMOzv4b19kPNMZsBl4XkUOzl04FXgIeBmZmr80EHursZySl++5T0DaZor+vMXyU+2AgXeMi3tkgumR2TvHGMciUzb0kRDbEy/IqUJfnz/0v4N5sq4w1wGwyCurXInI5sB44L8+fGYqLC6yT4aOcfYR0yQvu5UHihTbS1FQWQZiYdLecAuiTOQlxlMIXgX+KyHygofWiMebqzn6oMeZ5YJrlW6d29nd2hTjOtIvhI22LrGsWIbgX8orT2USbzDkbzGlTvN08fBRHKdwC/AV4EdD1dPJEnIfsmssN+qzYXCWp2qo2jDHOhbzihI+0JfTd8xTcm3tJiKMUmo0x/11wSUqIk6EYB62VuDL3roozLAtPrHHhpLHglsz67nHu92jz0pMQJ9H8tIjMEZEReSpJVYeL1pWLG+5cO7fCxQU2zhn32qzYnKXKykK3sfI2yhRZEuKYZB/N/v/FdtcMmU1o3YJ4iSO3JhLokznWLmxFCf048roY79ZmxboWonNx7iUhp1IwxhxcDEFKiYvho3hlcbpkdq1qw8VxEcsbUzYucuYUHJMX9I2LJEQdx3lS1A9mdzgflX+Rio+LYQIXFyzXdoK6qHjjhI/UbQZz7MjTeHNPl8xJiPIU/lNEvgM8TmYX8zagGpgAvAcYA3ym4BIWgXixbl0POV5Jqi6ZXXO7XZMX3DRw3Cv7de8eJyFUKRhjrs2eaXAumY1kI4B9ZNpc32KMebY4IhaeOPsUtD1kNyd/7vdoine7Fu4CR8NHjikFF9vWJyEyp2CM2QHclv3XbYnzkJvThqaWNBWpWEdQFJx4Tfx0DUzXqjZcC3dBXCtWm3cT/f365jTGGER07F+Ju2vcVXSscCUmbn8wTQtA9w1t6JHZNXnBTe8m1yLbkjY0xdmqXSRc83iT4pUC8RJHoGsBcNOKzf0eTTLHmvyKPBtwtSTVLQ+yx1Yf9STixGFB14N20iJ0zPJ2TV6IZyxo2tQYp5UI6BrLLo6LJMQ5T6FGRL4sIrdlX08UkQ8WXrTiEechgx+YXcXvaC48rm3EjOula6qkc82zSUocT+FOMt1R35F9vQG4vmASlYA4td2gazK52MTPtbr/OAtsc9rQHHcAFQHXFFmcxDjoktk1zyYpcZTCeGPMd4AmAGPMPkBHGUCeiOspaHK74yyw6nrGOGbFxgnFgK7WHPEqvBTJG1sp6JE5Xk5Bj7xJiaMUGkWkF9lyfhEZT7tzFboDcY8n1KT943Xw1DUw4yxYmpKgboYVc79Hlbwxh6gmg8y1MGhS4jTE+yqZXc0Hici9wInArEIKVWziurCaBmacBatF3d4KtyaTi8aCa7ttXVS8sTrRKjPIkhCnId6fRWQJcDyZsNE1xpg3Cy5ZEbFp/jIJWl26Bmb8yaRHKeR+jy7FG+99mkIFrlmxbuYU3PJ4k5JTKYjIlOyXm7L/jxaR/sA6Y0xzwSQrIrYFtndlObUN+/95ugZmvPfVN6XpW11YWeLiWsWUk6XKjp2nYOIWeSiyvF3zeJMSJ3z0M2AKsJSMp3BU9uvBIvJxY8yfCihfUbC5gzVVKYtS0DMwnVywHJtMca1YTVVeroWPuqunoGmtSEqcuMJrwGRjzDRjzFRgMrAMOA34TgFlKxq2h2w7ElJTaMPJBcu1yhgXS5UdK/uNGwbVFI6Jl1PQI29S4iiFw4wxy1tfGGNeIqMk1hROrOJiG5h9LEpBl7US732qFizHKmNcTDQ3O1YuGdfj1WSQuZa3SUqc8NErInIzcH/29fnAShGpIrt3wXVsnkJNZSpwTdPAbIlpxmoanPHOgNAjb3Pse6xnkY3ljSnavxI/fOTYPW7S1dk1CXE8hVnAauDTwLXAmuy1JjKH7ThPbE9B0WRqjtk1UtNkcs2KddNTyH3/NIU23OxQHO99DYpCoUmIU5K6D/he9l9H9uRdohJg9xRs4SM9D9nFBSuOd6NpwYqjxECbB5lb5qYWQ0vakCorvRUbe9e4qrkX30uvrghGHLQTpyHeRBF5QEReEpE1rf+6+sEikhKR50TkD9nXB4vIfBFZJSK/EpHKrn5GXGzPuHeV7vCRiwtWHO9GlxJzT/HGHRdaZG6K2TdKk7HQ5KAiS0Lchng3A81kwkV3A/fk4bOvIXO0ZyvfBn5gjJkI7AAuz8NnxMIW1+xt8RQ0LbAuLliu9YxxbYEF98aFa/ICtMQM3WpaL5IQRyn0MsY8BYgxZp0x5mvAKV35UBEZBXwAuD37WrK/84HsW+4Czu7KZyTBNjBrlFcfxV6wFMU1YzXxU3SP44cJ3LrHoOc+u6h4XZQ5CXGqj+pFpAxYJSKfAjYCw7r4uT8EPg/0zb4eDOxst0N6AzDS9oMiMgeYAzB69OguipHBllPoYwkfaXrIsRcsRcnxphgyNzanSacNZQri3XGT+VoWWEiyYOlQZPE9BR3yQvyqNE3jIglxPIVPAzXA1cBU4GLg0s5+YPaAnq3GmMXtL1veah0txphbsxvppg0dOrSzYuyH1VNQHj5y0VqJ63ZriR87GdpwTObYOQUl8oKbcy8JcTyFscaYhWQqjWYDiMh5wPxOfuaJwJkicgZQDfQj4zkMEJHyrLcwCnijk78/MdbeR1ZPQY+14mJcM4kVW1O0MoNwXLO6wS6zSLDVupad7q4pMUhg3CiSOQlxPIUvxrwWC2PMF40xo4wxY4ELgL8YYy4CngbOzb5tJvBQZz8juUzBa9ZEs6JQjG3y23dh61mwbAtApaWDqxZFZpPXVsapafLbwop9FJdX28Zxb8vGUS3ygl1m+7jQI3MSQj0FEXk/cAYwUkR+1O5b/chUIuWbLwD3i8j1wHPA3AJ8hhVb9ZFtgdViXUH4hrs9HZr4aVlgIWQBqErRWLf/5NGifMMWrN31ervn2vIgvav0dvy1ydunupy9HcaAlpAi2HMKfarK2bVv/wYPWsZxUqLCR28Ai4Ezs/+3UktmZ3OXMcY8AzyT/XoNMCMfvzcptg00tuojTQ/Z6ilUl8Pu/a9pahthVWTV5eyo238yaVmwbFZ33+qKgFLQpHjjhkK1yGxbYDPNKPc/3LFBkdUdZpB1VAqaFFkSQpWCMeYF4AUR+UV3OTchDLunYJ9IWvqZ2BYs7Z1drQuANbShQ+Ywz6YjWuSFJM0ddSyyceVVNY5t3o1yIzIJUeGjF3n7XObA940xxxROrOJij3WnKC+T/RaGtMm0CKgsL71SsC1YfZXvrbDd577V3WHB0iEvRHiQHdAyLppiLrBa5IVwj7cj3bH30QeLJkWJsYWPysqgV0XwoJ19TS1Ulpf+eEt77FhvmADCLG+9VmFYfL4j2kN0mr2xJG3rtXjptv023clTCF3dsruX1xlj1gH1wNHZf/uy17oNtqq4MhGqLM2stEwme/VRReCaFqs7nTaBKi8R9xYsm2ejRYlBeBK0I1qsWJu8VRUZL709aQONMfc0FJq4noKWcZyUOA3xPgIsAM4DPgLMF5Fzo3/KLWw5hVSZ0KsyeHu0PGhr6aHieLdNiZWXCVUVektS7dVHuid/3AVLixVrk7e8TOhlNch0KAVrTkH5ZtckxNm89j/AdGPMVgARGQo8ydt9ipzHGj4S+8DU8qBdix2H1fzb7rGWcIybyfx4IS8t48K2wGaMhWDotqGpBXoFveFiY/XGFOfGkhInOF7WqhCyvBXz55zB5pZWpMTaC13Lg7aXHupdsGwTqaKszHqP9cgcv5In7rGShca221ZzNY/Vg0zZvXTNMmtPjichjqfwuIg8AdyXfX0+8MfCiVR8bNZKRSpkwVLidsevPtKrxFIp98IEleVlVKbKAoZEQ3NaxYEqtiSo5h3CNm8sVSZUl2uWOW4VnY61IilxTl77nIicA5xEpnHdrcaY3xVcsiJis2LLQxcsHQ/aahEqToLaSg/Ly4RqxTkFW7O28lRG5o5KQcspW42WBHL/GksBgpKNVfZxodyDjFmVpkXepORUCiJyLfAbY8yDRZCn6BhjrAMzE9rQm2i2hbxs1UdaWlHb5K0M8ca03GNbhU6rzLZdzQOKJVgILWkTqKQLrfBS4vEmSzTrkDluhZcWeZMSJzfQD3hCRP4uIleJyPBCC1VMwppblYUMTC3a32YRVpWHKDIFVqFN3spy3UrBpsiqysvopTQcY73HqRB5FYwJsN/j8lSZ6qo0qzdmSYBr2tSYhJxKwRjzdWPMkcBVwIHAX0XkyYJLViTCqh8A62TSPDDDF9nSD05bM8HK8jLVOYXQe2yNd5d+XLg2JsDujVWFjAstVWk2mftZlIIWeZOSpIpoK7CZTPVRV09eU4MtMVeRbedcpTjZ1RBmxSr1bsIWLKviVRLaCF1klRoLDS1BGapClJjme1xVoVuRxfcUdNzjpMTZvPYJEXkGeAoYAnysO/U9arI84PJUuKegwSIEuxWiORwTFtrQGu6CMJlTVFvanGi+x9ZNmErusc2DrCpPqS1AaEkbe+VfT6o+AsYAnzbGPF9oYUqB7QG3egqqLawQT0FrGW2S0IYGeSEkOR6aUyi9zGH32ObxamlFbU3mh4YVdd5jzR56Z4hTknpdMQQpFbayw4q2nIJOixAirFiLhaXhcCBbuKuyPGX3bJT05UmWUyi9zGFKTHN5Z3jBhE6Zk+ZttDTxS0K32pncGax10llPQbP2D43RWy1vBQtWAgtLS7lkkmoeDd6NTd6KsOojJeM4LNGsNadgD3eVkSoT69GyWhoPJqHHK4XmkA1KQEiXVB0POYlVqGEBCAsT2D2F0ssLYd6N3jxIuGdj93g1tOZoDFlkbeNCQzWPXYllZNW8rykJPV4p2DyFyghPQctDtsWENcc2rZ5CKsyzKb28kDC0oUDmMM+mPFWmthV1eE5BZ6I5zBgD1Ho3SenxSiGsxQUoDx+FDE7bph8Niiw8FqtTXgjfW2G1YhWECaIWLK37QezGTUiuScG4CDPGQPe+piT0eKUQFoeFMM1f+ofckjaB9gAiutsDJFlgNSxWEL4AaK1Ksz3n1vtrC4VqCMckS46XflxEegpKNzUmpccrBdsC1PpwtbbvrWtsDlyrqUghorfdd51l0aypLKeqvIyOxRmNLWlrT5xiY7vPvavK1Val7W0IytDaIVVr3b/tvoUnmksvr035t64XWjc1JqXHKwXbQ+vVNpF0WoTWBTbbkEtryMu6wFZmFZlSC2uv5T73rixXWy5ZZ5EhalxoNRZ6V5WrDSvaDYXseqF0U2NSerxSsLvc4eEjDbHjvQ32BRb0VkDYrNjWBUujzM0t6UBoUSQjq1pvLHJc6FS8ezETrmoAACAASURBVKwyl6sNg9rkbVO8ikt/k+CVgk0ptIaPHPIUelW2LrBKrdgQTwF0ejdWqzsyRFf6e2zzbN4eFzrDRzZFVlMVlmhWoHit3mOrp6BT5qQUXSmIyEEi8rSIrBCR5SJyTfb6IBH5s4isyv4/sBjyWJVClHWloB492lPQOTBtC9bbnoI+mesiPButVuy+CMWrUZGl08aqfDWH6KxzL8JT0GBEJqUUnkIz8BljzOHA8cBVInIEcB3wlDFmIpnme0VprxGZaHbIU9C+YLkW2ojybLRuXrMqXsXGwr6mFjrun6uuyOwO1jqO7cn8iDCognGRlKIrBWPMJmPMkuzXtcAKYCRwFnBX9m13AWcXQ56onEKVJXHUkD3JrJTsdcwihLDQht48SFi1FCg2FmyhmIiwYqnvsW0c91GcZ4LwijTQW5iSlJLmFERkLDAZmA8MN8ZsgoziIOTMBhGZIyKLRGTRtm3buiyDtfoo+3DLyiRUMZQSmwvbtmApLaOtrQ9fADRu+tm9rylwrY/icBfALpvM1a2KTN8iay0+aFViSk+3sybGq8INslKvFZ2hZEpBRPoAvyXTlnt33J8zxtxqjJlmjJk2dOjQLsthDR+1e7gaF6wddcHJPyB7OLvW8s6ddY2BawNrKgGdCbrIe6zQ6ga7zG33WKHMtjHReliN3nEcvMetMmv1IJNSEqUgIhVkFMK9xpgHs5e3iMiI7PdHkDnpreDsawpq/vZWisbKmB17bQtsdjJZlVjprZXtFpkHRMpcasUbocQUWt0QpnjDFVmpx4VtTAzsnbnHFSmhQ7smmtPG2uq+mNhkHtRb97hISimqjwSYC6wwxny/3bceBmZmv54JPFQMeXbvCyqFfu1OUdJoYVkXrN7hVnep2xnsa2wJuNEVKXk7HKPQKrQp3gG9IyxCBZPf7t3o9RSsC2xWiYnoTDZHGQtax0VSSuEpnAhcApwiIs9n/50B3AC8V0RWAe/Nvi44u+uDE6lf9dvnrWpMHkWFCXSGu2xeQmXb4SMa20YkD8WU1oJtSRvrWH475KWvMibKuAGde26ivBvNrfaTEOc4zrxijHkWCDuK6NRiygKw25IA7dfr7duiMUFnDx/pdWGtE6mmneLV6CkkDsW0lPSUrZ11jYHyzr5V5W3NHTUeZrR9b1CJDaqJVgqlPEbUGGMdF4MiPIVSj+PO0ON3NNdaKjZyeQql1v5bausD1wb3iXBhSzz5t+wOyjukT1Xb1/ZNP6W9x5t27Qtca5VZ4ylbm3YF7/HQfm/fY43jeKtlXAzq014p6DJwtu9tDJy/0qsiFdkrzSsFB7G53H3bKQVtccKWtGGzZQEY0b8a0Dn539gZXGBHDujV9rXGMIFtkR058G2ZtZ1bkfse6wsfbbDIfGA7mauUVaVttN3jdmNC21rRWbxSyBE+0rZgbattCFgr/arL2xSZbV9FY0tpN9xt3GlRYvtNfl0LbDpt2GSRuf2Cpa3rqE2JHdg/eoEtZSgGYOOO4CI7Kocis53LUSxs8uZUvF4puMXehuZAJ8yOB9XYLELbwTzFYuPOusC1kQNr2r4W0bfhbsMOi8wDqtu+1rbpZ2ttQ+Awlb5V5bkLEEq4ALy+PXiPR+x3j3V5Cs0taTZbwkftLW9tXm8uT0Fj2W9n6NFKYWttQ+Da0L5V+yULrRZWCSfTq9v2Bq61t1YgrD1H6WReuaU2cG30oN5tX2s7pP0Vm7yDa/Z7rS2paJN57OC377E2T+G1t/YGDlIa1LuybUcz6BvHq7bsCVwbNVCvZ9NZerRS2GZRCsP6Vu332h7aKN1kenlTcPIfMrzPfq81WViNzWnWWBTZYQf0bftam2ezcnPwHh/aTl6we5AlldmiFNrLrM1TeMkyjg8dvv891jSOAVZsDjZeaC+zNsXbWXq0UthqqeIZ2rd6v9f2yV+6ybRiU3BgHj6i336vNcm8ckstzR0swmF9q3LWo5fS6l7+xq7AtcM6KAVNZbTbahvYsnt/AydVJowbGu0plHKBtd3jwDhW5Ck0taR5xWIstJdZ07zrCj1bKey2h4/ao0n7Nzanef71nYHr8Ras0sj87zVvBa4dcWCcyV8aeY0x/HvN9sD1eIq3NDLPXxu8xxOH9dlv7GoLbcy33OOO40KTp7B0w87A8x1QU9FW9Qe61oqu0KOVgq1+vmP4yD6ZSvOgn1u/I5DMHFhTwfih+4ePNFks/3w1uGAdd/Dg/V5r8hTWvLk3kACtSAlTx+x/5pNVkZVI5n+sDt7j48ftf481LVi79jWxdEPQuDnu4EH7vdZUlfbsquA9njF20H75R20hus7So5XCureCFRvtE0cQ5naX5kH/6aUtgWvvGD+Ysg6dw6z9j0qgyHbVNfH3VcH25sePyz35S6V4H126KXBt0qgB+yVAIaSlQQlkbmpJ88TyzYHrgXusyFB4YvlmOlZIHzSoFwcN2j+Zr6UqzRjDoy++Ebh+4oQh+72uTJXRcUN7U4sJJNS106OVwmtvBROgY4f03u+1lgWrqSXNQ88HB+Y7Jwbbh2vZWPXwCxsDeyqG9KnimFED9rumxVNoSRseXLIhcP3Uw4cHrmnxFJ5+eWugjUhVeRkndRgXmkIxDywO3uP3HBo8PsXeS6j493jZxt2stFQedZQ5vBzcLW+hxyoFY4zVUxjTofRQy0N+6Pk3eHPP/jmQylQZZxw1IvBee6fU4i4AzS1pbv37msD1D00aQaqDZ6NF8T6+bDOvWcbEmcceGLimwYo1xvDzv74auH7aEcPbOtC2omUcL163nQVrg/mEsyePDFzTMi5+9szqwLWpYwYGypRBV5ius/RYpbCttiEQn6+pTDG0T8ecQuknf31TCz/5y6rA9VMPH0b/do3lWrF6CkVeAO759zpe3x7c7HPe1IMC1zR4CvVNLdz4p1cC10+aMCSwDwR0xLsfW7aZJeuDsfnzpo4KXNMQ2kinDd/648uB6xOG9WHyQQMC1zWMi8XrdvDYsmB4znaPoXvkFXqsUnjZUl42ZnDvQJdLDWGCnz3zqtWC/di7xlnfX2pPYePOfXzvTysD19996NBAhQnosAh/8pfVrH0zGE785LvHW99v39RYPJl37G3k648sD1w/YkQ/Tj4kGFLUENr45YL1LFq3I3D9EyePt3aXLfW4qG9q4UsPvhi4fkC/aj48JejZgPcUnGaZrU66Q2kn2OOaxRyYC9Zu52dPB93XkyYMYcrogZafKK2n0NDcwifvXWI9y/aaUydaf6bUFuEzr2zlp5YQwYyDB/GO8YMtP1HaqrR02nDtr58P7E0A+PRpE0Pbd5cyr7Bs4y7+9w8vBa6PG9KbsyzhOSj9TvevPbzculP8qveMty7+UHpFlg96rlLYGFQKR43sH7hm9xSK85C37q7nql8uCWz+Ki8TvvzBI0J/rlTWijGG//e7Zbxg2Utx4YzRTA5TYiETyXQ8IKAArN5ay6d/9XzgLILyMuH6s48KXWDt97jwC5Yxhm8/8TLPvBKs6jr1sGG894hgUryVUnkKW3fXc+U9i609w67/8FGUW9qQQ0iIrkjGzT3/eo37F74euH70yP589LgxoT9XagMnH/RYpfCiRSkcPSqmUijCwNxZ18glcxdYW3FcefK4QNuF9pTCUzDG8PVHXuI3lsqSA/pVc93ph4X+bFmJzidY99ZePnrbfOth7Fe9ZwKHDI+4xyWyCH/69Gpu+Wswgd+vupyvn3Vk5CE/pfAU3trTwEW3z7c2k7tg+kGcMH6I5acylOqQnQeXbOArDwdDc5WpMr51ztGBQon2eE/BUXbVNQWSoCKZeGxHSpFo3tPQzMw7F1pd16ljBnLNqYdE/nwpcgrf+9NK5v3ztcD1ipTws4unWBPi7Sn2ZHpj5z4+ett8a1PEd04cwtUhoa5W7OOisIp37rNrudGSqwH44QXHMmpgsBqmPcVOju/a18Sldyxg1dZgOecRI/rxtTOPjPz5UiRt//jiJj77mxcCniPA18480hpNaI/3FBzFlk8YN6Q3vauCp5MWeyLVN7XwsbsWWUMwQ/pU8tOPTqHSIlN7iu0p/PTp1fzEkvcA+OqHjgzNfbTHmrsp0H3eVtvAxSHW66iBvbjpgsmR1iAUv1HifQvW8w1LTB4yeYRTDgsPG7VSTANnT0Mzs+5cwPI3gr26BtRU8POLp1rlaU+xw6B/eXkLV9/3XGBjHcD50w7iwhnByrmOeE/BUayhoxALoJiJ5qaWNJ/65RL+ZekX1Lsyxe0zp3NA/2rLT+5PMT2FO/+xlu8+ESzlBPjMew/h4uPD46/tKdZk2rG3kYtvn88aS6XRAf2q+eUVxzOoXbO+MIq5Q/j3z23kS78LVsEAzDphbGgCvyPFMnDqm1q44q6FPGcpl+1bVc49lx1nrfHvSDE9hX+ufpOP/yKYvwN43xHDuf7D4fml9nhPwVHiJpmheItVOm347G9e4MkVW60y3D5zOsdaarltFGvB+tXC9Xz9Ebv1+vGTx/OpUybE/l3FOLWqtr6JmXcusIblBveu5BdXxFusoHhNBx9ftpnPhIQzzp92EF/54BGxFisojqfQ2Jzm479YbG0q2KsixZ2zp1tzdzaK1dl18brtXHH3Imsi/ORDhvLjj06mIiQZ3pFirBfNLWl+/NQqa4VfPvBKIUuYUihGWZwxhv/30DJrG4vyMuHmi6eElkbaKIbb/cgLb3CdpYYbYOY7xvCF0w+NvVhB4ev+9zW2cPm8RSzdEHz2/arLuefy45gwrI/lJ+0UQ/H+deU2/uu+JdYNZmdOOpBvnnN0oO9VFIX2FJpb0lx933PWyqjK8jJunzmNaWMHWX7STjE6uy7buItZdyykrjH4e487eBA/v3hqaPmpjUK35kinDV988EW+9+eVXHTbv9nRocVJPuhxSmF3fZN1I9iRlk1VEFYWl7/FyhjDDY+9zC/nrw98TwR+cP6xseLF7Sm02/3kS1u41lLGCZmdnl/9UHQVjI1CegoNzS3MuWcRC14LWq+9K1PcddkM66a6KAqteBes3c6V9ywK9I4CeO8Rw/neRyblzHt0pJCh0HTa8PkHlvK4pTlfeZlw80VTAg3kclFoT2HllloumTufWovFPXn0AObOmk6vyvgKAQrrKRhjuP7RFW0Vfi9s2MX5t/7L2u25K/Q4pfCSJfE1bkjvtoPvO1JeJnScey1pQ3NLfh70T59ezS1/C5YYAnzzw0fzoUn2jT1RFHLB+sfqN/mkZe8EwAePGcEN/3lMIuu1lUJ5Ck0taf7rl8/x91VvBr5XXVHGHbOmh+6fiKKQinfphp1cNm+hdQF858Qh/CRBOKM9hToYyBjDlx9axoPPbQx8r0zgpgsmW5sK5sLqKeTJUHjtzb1cdPt8dljKkY8Y0Y95s2cE+kfFoZDNKH/01Gru+Mfa/a6t3LKHj9zyL/ZZPJ3O0uOUgi10dGREmVmmPUBhFqx5/wgvMfyfMw7nwhmjO/V7C7VgLV63g4+FxF5PPWwYPzj/2MTWayuF8BTSacPnfvOCteV4RUq45ZJpHDcufliuPYVSvCu31DLzjgXWePGMsYO49ZJpicIZ7SnEwUDGGL712Mvca/F0Ab5z7iQ+cEywaWMcrJ5CHubdGzv3cdHt8617gCYM68M9l8+gf6/oEuowCtW2/o5n1/KDJ+1rxewTxib2aKJQpRRE5HQReUVEVovIdYX4DGs+IUfooBCT6ZEX3uBrIUnaq0+ZENrXKA6FWLBWb61l9p0LrLHXE8YP5qcXTemU9dpKIUIF//uHl/i9JU+TKhN+fOEUa4+guBQiTLBx5z4uDrFejxnVn7mzpnVp8tur0rqmeOc+u5ZbQzzdb5x1JOeGNI6Lg+0eN3Zxp/uuuiZm3rHAWo48elAN915xHIM7NMVMQiFyTb9/bqO1RQhkKvxmnXhwl35/R9QoBRFJAT8F3g8cAVwoIuG9HDrJMkv4KOeGlDy73Yte285nfvOC9XuzThjLte+N3pyWi3z35dlW28CsOxeyuz5ovU4ZPYDbLp2Ws+Y8F/lOKt7x7FrrZjoRuPG8Yzj9qAM6/bsh/6WHtfVNXHbnQutmuonD+nDX7BmhIc645Nu4eXTpJq5/dIX1e198/2Fc8o6xnf7dkN3pnkfl29DcwpW/WGTdTDeifzX3XnEcw/vlLvmOIt/l4P9e8xafe8C+VnzsnQcnqvCLixqlAMwAVhtj1hhjGoH7gbPy+QF1jc28ui04II46MFop5HMyrX1zb2gI5typoxKVGIaRz9Pi6hqbufyuhWzYEbSsjhjRjztnz7Bu+ktKPj2FJ5Zv5huP2i2r688+ig9P7rz12ko+PYWmljSfvHeJtVR29KAafnHFcQyMsXciF/k0bha+tp1rf/289XtXnzqRK0+2d5dNSr4qpozJJMJtpbJD+lRx7xXHBU5+6wz59BRWb61lzt32YoMLph/El844vMtrhQ1NSmEk0L4D1YbstbyxYtPuQMXMQYN6dbIFQ/IHvX1vI7PvXGAND7zviOHckLDEMIx8eQotacPV9z1vLeMcO7iGu7sQe+1IvjyF51/fyTX3P2etjPrC6YdxUUQzsyTkq82FMYavPLTMmggf1rcqL9ZrK/lKgr66bU+oYTPzHWO49rR4m+nikK+9Fd994hVryXefqnLuvmwG44bGL0eOIl/7V6K88zOOPoD/+/DRBVEIoEsp2P7CwNQWkTkiskhEFm3bFqyHjuJFy+KWy0uA/MTo65tamHP3Ims57KSDBnDTBZNDu0UmJR+egjGGb/zhJZ5cEUzSDupdybzZMxjShdhrR+z13cnu8fq36rg8pGrn4uNH8/GTO5+n6YgtrFHflDze/fO/ruG+BcFunDWVKe6YNT0v1msr1Xnwbt7c08CsOxdYmwi+94jhfKUT5chR5MNT+OX89fzsmeAJda17gJKWI0eRD09hX2MLV9y9yOqdTx0zkO9/pPMFHXHQpBQ2AO2bi4wCAqrdGHOrMWaaMWba0KHJEoUV5WWMH9p7vxOocuUToOuWd+tuZdsBI6MG9uL2S7uWQOxIPjyFO/7xmjUmX1lexm2XTgucZd1V7AtW/Mm0q66JWfMW8JZlM88phw3ja3lerFJlQkUq+PsaE5QqP7p0E99+PHgSWZnATz46OdbYTEJXN1bVN7XwsbsXWU/Um3TQAH4Uo2dUUrrqKfxt5Ta+/NAy6/e+dc7R1jPOu0JXPYWWtOHTv3rO2vts7OCavOTvctH1YHD+WAhMFJGDgY3ABcBH8/kBFx03houOG8PehmZWbNrNso27YpUkdtXy/v6fV/KHpZsC1/tVlzNv9nSG9s2fxQ1dl/fJl7ZwvSUmLwI/PP9Ypo5JXtefi654Cq2tFdZsC/YzOmpkP358Yf68sPZUl6doatnfva9vSscqGV28bkdoTP5rZx6ZeMNiHLpiLLQaNrZ+RqMH1TB3Zn4Nm1a64ims3FLLVffad4R/+rSJnDctd4O7pHTVU7jhsRU8sTzonQ+sqeDO2TNi9eXqKmqUgjGmWUQ+BTwBpIA7jDHBpuZ5oHdVOdPGDoq95b4rScUHFm+wdhBtrZOfMCy8Z39niTq0Jpe1vPyNXVwdEpP/0vsP54yjO1dznovO5m2MMXz598usTQQP7F/NHTOn5yURbqOqooyOxUIZmaPzLOvfqguNyV9+0sFc2sWqnTC6Yiz88Em7YTOwpoJ5s6fnNZTYns56Cm/uaeCyeQutu5XPnToqdhPBpHSlYOLe+eu47e9rA9dbvfOD8+ydh6FGKQAYY/4I/LHUcnSks9p//pq3+OKDS63fu+GcYxL1M0pC66E1HUMZDc3pSNdzy+56Lp+3yLoX4ZLjx3DFO/NbD92eznoKt/xtDb9aFIzJ960q587ZMxiWpyStjc7kmnbWNTJr3gK2W8Jc7z1iOF864/C8ydeRznoKv3tuAz/6S9CwqUyVceul0/KWpLXRmU2Nrfk7W0z+hPGD+WYBk7SdLZj4+6ptfOUhuw38vfMmJeoZ1VU05RTU0pnJv/bNvVz5i8XWcrKrT53If3ZhU08ckpbRtpaebrb0UXn3oUP56oe6XiobRWdyCo8v28QNjwVj8qmyzME+UafT5YOkxkJjc5or77GHuY4Z1Z+bLihsArEznsLC17bzhQfsjQ+/c+4xTC/wYpXU8m4tPV1iCXONG9qbmy+amvM8kq7QmbVi1ZZaPvkLe5jrc/9xaKda3XQFrxRikLRtxPa9jVw2b6G1QuPMSQfmtWQvjCRnCGe6Wz7Pso3BjX2HDu9bsJh8e+yH7IRPpufW7+DTv7LH5P/3rCPznkC0kWTBMsZw3YNLmb82WCc/ckAvbp85jZrKwjruST2FV7ftYc7di6zJ86tPmcDZk/NaMW4lqeV9459e4eEXgqWnA2squGPm9Jzl510l6VqxrbaB2SFhrvOmjuKT787Pfo8keKUQgyTav/XEqbWWQ1ymjhnId849pqAWdytxFwBjMq14baWnQ/pUMXfWtC7vpI1DEk9h5ZZaZoeUnl5x0sF524uQiyS5pu888QoPLgk2jMuEuaYzrG/hwlytJFFim3bt45KQlhsfPGZEl3fdxyWJzHc8u5afPh0sPW3N3+W7Ys5GkrWi9bhSW5jrHeMGF3QvQhReKcQg7uRvjWXaNnsdNKgXt16S+wjCfBG3auOGx15ua8Xb8edvnzkt57m/+SJuTuH17XVcMne+1Qs77fDhfLGAMfmOxO3i+bNnVnNzaJ38VA4ZXtgwVytx5d2xt5FL5i7gjV3BUOLk0QO48bxJRVus4noKUf2BbjjnGGYcXJyYfNwcSOZ8j4Ws2GTp2jy0Nz+/uLBhrihUJZq1Eqcypr6phavuXcI/Xw1WwfStLufOWdO71GgrKbmqNowx/ODJVda23a3nOMQ96S0fxPEUNuyo4+K589myO9gf6MgD+xU8Jt+RON1z7/rna3zncftxpdeffRQnTUx2xkBXiCPv9r2NXDJ3Pqst/YFGDezFrZcUvk6+PXE8hUdeeCO0l9jVp0woeP6uPXHu8b7GzPketn1LxQpzReGVQgxyWbF1jc3MuXsxz64OtiqoKi9j7szpBSk9jSLKUzDG8M0/rrCWvwF87UNHFqz0NIxc93jtm3u56LZ/W63XsYNrmJenHkxJsN3jfe3u8c+eeTX0/OpPvWcCF3SyNXpnydXmYlttAxffPt/ag2lIn0p+cflxed9Tk4tcMj+weAOff+AFLDlaLpwxumhhrlZylYPX1jdx+Tz7gU+tu9iLEeaKwiuFGEQtsFt31zPnnsU8b9mB2LqNvliua3ts1lxdYwv1TS186cEXrQeiQGZTz8wTxhZYuiA2t3tvYyb5tnjdDq68ZzFv7gl6CMP7VXFPCRYrwNr3aUddI80tab75x5cDB6K0cuk7xvCZ9xV3sQL7mGhVYiu31HLFXYtYvz3YhqVvVTnzZs8oyWIVVjCRS+mefuQBXH/2UUWPyZeVCdUVZQFvprahmT31zcy5Z5G1oKMyVcatl0zr1IFP+cYrhRjYLNDd+5pYsn4Hn/jFYms4I9OieVJBdqbGwdZV88WNu/j+n1daFRhk2nYXalNPLmw7NbfVNnD/gvV85aHl1gqYATUV3H1ZfrpbdoYhFkW0cnMtM+9cwD9WB8OIAOdMGZn3lhtx6V1ZTqpM9it9rGts4aHnN/I/v1tmPdinV0WK22ZOy3vLjbj0qw7OvXXb67jm/uetVUaQOZ3uh0UOJbZnSJ+qQPL4j0s3ceOfVloNm8zpdMcWNZQYhVcKMRjeLzj5H3xuIw+98Ia1tjhVJnz/I5M469jCl+yFMdxSzRJmVUGmN3uhWvHGYUjvqsCCVVvfzHUP2mvkW9sdF3ovQhS2Xbx3/Wtd6PvPmTKS73TyuNJ8kCoThvWtYlOHENw199tLe/tUlXPHrOkl8XRbGdG/V+CardtpK+85dCg3X1y8gg4bw/oGlULYOK5ICT++cDKnH1XccG0UvvooBmHlgjaFUJkq4+aLppRUIYBdkYVx9SkTSqoQION2D4sZAhrRv5pfX3l8SRUCkChkNeuEsdx47qSC7/fIxQH945W+9u9VwT2XzyipQoD48kImZHRLkRPhNuKOi6ps+wpNCgG8pxCLETEH5pA+VdxyyRSmjintRAIYMzh3SKUiJfzf2Ufzken5bwzWGUYPqglYsR059qAB3HrJ1IK2r4jL+KG5Y+xlkjnHYc67xpVU6bYyZlCNtalde8YN7c3cmdOL1msnitExQ4MfP3k8n/uPQ0sWMmpPpu1HcN9Pew7oV80tl0xlUhEr/OLiPYUYDOpdyaiBQTe2PZNG9efhT52oQiEATB87iKj5MaxvFfd97Hg1CgHI2Qvq3KmjuH/O8SoUAmR2e9ti3q0MqKngrstmcOXJ41UoBIDjc3QFPu3wYfzukyeqUAiQsbonDAvvrVRTmeL7H5nEde8/TIVCADguh3c1bcxAHv7UiSoVAnilEAsR4ZyQLf3lZcK1px3CA584gQMHRCuOYjK4T1VoCOvMSQfyp2vfVdQmW3G4YPpoaiztlwfWVHDzRVO48bxJJQ8NtKc8VcaskEqt0w4fxhOffldR2m0k4UOTDrSGN2oqU3zrnKO57dJpeTtNL19ccZK9EePUMQN57Jp3cs6U4u1DiMNJE4Yw0aLIKlLC508/VJVhY0OSnhSliWnTpplFixYV5bP2NjQz844F+204edchQ/nSGYdx2AH5O7kpn+xpaOZzv3mBx5dvxpjMbtRrTzuEdx2ia6Fqzz9ffZMvPfgir71VR9/qcj4y7SCues+EovSR7wxNLWm+8YeXeHDJRvY0NHP0yP5c9Z7x/MeRB6jxDjry0hu7+e9fP8/Lm2vpVZHirGMP5JrTJlqTuhowxnDTU6u47W9r2NvYwvihvfn4yeM5Z8ooNd5BR9a9tZfP/uYFFr62g8pUGe87cjjXvvcQxhewo2wSGyWF5gAAB1BJREFURGSxMWaa9XteKcSnuSXNv9a8xa59TRw6vC8Ti9SeoKu07qnQZGVHYYyhJW0oEylZpU5naGhuiXXAjhYam9NUpESt8upIc0saA1SUOFmfBGMMaYM65RWlFHyiOQHlqTJ14YA4uKIMWhERyi1HXWrHJYUA9nOmNVPqyq3OICK4NpTdu8sej8fjKRheKXg8Ho+nDa8UPB6Px9OGVwoej8fjacMrBY/H4/G04ZWCx+PxeNpwep+CiGwDwttSRjMECJ6K073xf3PPwP/NPYOu/M1jjDHW+nqnlUJXEJFFYZs3uiv+b+4Z+L+5Z1Cov9mHjzwej8fThlcKHo/H42mjJyuFW0stQAnwf3PPwP/NPYOC/M09Nqfg8Xg8niA92VPweDweTwe8UvB4PB5PGz1SKYjI6SLyioisFpHrSi1PIRCRg0TkaRFZISLLReSa7PVBIvJnEVmV/X9gqWXNJyKSEpHnROQP2dcHi8j87N/7KxHReVpPJxGRASLygIi8nH3W7+gBz/ja7JheJiL3iUh1d3vOInKHiGwVkWXtrlmfq2T4UXY9WyoiU7ry2T1OKYhICvgp8H7gCOBCETmitFIVhGbgM8aYw4Hjgauyf+d1wFPGmInAU9nX3YlrgBXtXn8b+EH2790BXF4SqQrHTcDjxpjDgElk/vZu+4xFZCRwNTDNGHMUkAIuoPs953nA6R2uhT3X9wMTs//mADd35YN7nFIAZgCrjTFrjDGNwP3AWSWWKe8YYzYZY5Zkv64ls1iMJPO33pV9213A2aWRMP+IyCjgA8Dt2dcCnAI8kH1Ld/t7+wHvAuYCGGMajTE76cbPOEs50EtEyoEaYBPd7DkbY/4GbO9wOey5ngXcbTL8GxggIiM6+9k9USmMBF5v93pD9lq3RUTGApOB+cBwY8wmyCgOYFjpJMs7PwQ+D6SzrwcDO40xzdnX3e1ZjwO2AXdmQ2a3i0hvuvEzNsZsBG4E1pNRBruAxXTv59xK2HPN65rWE5WC7XC8bluXKyJ9gN8CnzbG7C61PIVCRD4IbDXGLG5/2fLW7vSsy4EpwM3GmMnAXrpRqMhGNo5+FnAwcCDQm0z4pCPd6TnnIq/jvCcqhQ3AQe1ejwLeKJEsBUVEKsgohHuNMQ9mL29pdS2z/28tlXx55kTgTBF5jUxI8BQynsOAbJgBut+z3gBsMMbMz75+gIyS6K7PGOA0YK0xZpsxpgl4EDiB7v2cWwl7rnld03qiUlgITMxWK1SSSVI9XGKZ8k42nj4XWGGM+X67bz0MzMx+PRN4qNiyFQJjzBeNMaOMMWPJPNO/GGMuAp4Gzs2+rdv8vQDGmM3A6yJyaPbSqcBLdNNnnGU9cLyI1GTHeOvf3G2fczvCnuvDwKXZKqTjgV2tYabO0CN3NIvIGWSsyBRwhzHm/0osUt4RkZOAvwMv8naM/Utk8gq/BkaTmWDnGWM6JrScRkTeDXzWGPNBERlHxnMYBDwHXGyMaSilfPlERI4lk1ivBNYAs8kYe932GYvI14HzyVTYPQdcQSaG3m2es4jcB7ybTHvsLcBXgd9jea5Z5fgTMtVKdcBsY8yiTn92T1QKHo/H47HTE8NHHo/H4wnBKwWPx+PxtOGVgsfj8Xja8ErB4/F4PG14peDxeDyeNrxS8PRIst1FP9nu9YEi8kDUz3Ths84Wka/k4ffcKCKn5EMmjycMX5Lq6ZFk+0H9Idtps9Cf9U/gTGPMm138PWOA24wx78uPZB5PEO8peHoqNwDjReR5EfmuiIxt7V0vIrNE5Pci8oiIrBWRT4nIf2ebzv1bRAZl3zdeRB4XkcUi8ncROazjh4jIIUBDq0IQkXkicrNkzrpYIyInZ3vnrxCRedn3pLLvWyYiL4rItQDGmHXAYBE5oDi3yNMTKc/9Fo+nW3IdcJQx5lho8xzacxSZzrLVwGrgC8aYySLyA+BSMjvibwU+boxZJSLHAT8j03OpPScCSzpcG5h935nAI9n3XAEszO5QTgEjW70YERnQ7meXZN//28792R5PNF4peDx2ns6eQ1ErIrvILN6QaRtyTLb77AnAbzJdBgCosvyeEWTaW7fnEWOMEZEXgS3GmBcBRGQ5MBb4KzBORH4MPAr8qd3PbiXTHdTjKQheKXg8dtr3zUm3e50mM2/KyPTwPzbH79kH9A/53e1/b9vvNsbsEJFJwH8AVwEfAS7Lvqc6+zs9noLgcwqenkot0LezP5w9m2KtiJwHbefkTrK8dQUwIcnvFpEhQJkx5rfAl8m0w27lEGCZ9Qc9njzglYKnR2KMeQv4RzaZ+91O/pqLgMtF5AVgOfZjXf8GTJZ2MaYYjASeEZHnyZzV+0VoOx9jAtDpDpgeTy58SarHU2BE5CYyeYQnu/h7PgxMMcZ8OT+SeTxBvKfg8RSeb5I5YL6rlAPfy8Pv8XhC8Z6Cx+PxeNrwnoLH4/F42vBKwePxeDxteKXg8Xg8nja8UvB4PB5PG14peDwej6eN/w8/LYWPLZ1TXQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot neuron potential\n",
    "fig = plt.figure()       \n",
    "plt.plot(tvec, Vy[:, 0],linewidth=5)\n",
    "plt.xlabel('time (ms)')\n",
    "plt.ylabel('voltage (mV)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note, you can go back up and change parameters and rerun all the scripts to see how this changes the evolution of voltage in the model. You will use the accompanying python code hh_spike.py as part of HW1 to study the HH model."
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}