{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Subthalamic Nucleus neuron\n",
"## Validation of the model implementation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Imports"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import logging\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from PySONIC.core import Batch, PointNeuron, NeuronalBilayerSonophore, PulsedProtocol, AcousticDrive\n",
"from PySONIC.neurons import OtsukaSTN\n",
"from PySONIC.utils import logger, si_format, pow10_format, Intensity2Pressure\n",
"from PySONIC.plt import CompTimeSeries, GroupedTimeSeries\n",
"from PySONIC.constants import *\n",
"from PySONIC.postpro import computeFRProfile\n",
"\n",
"logger.setLevel(logging.INFO)\n",
"\n",
"pneuron = OtsukaSTN()\n",
"ref_Vm0 = pneuron.Vm0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Spontaneous spiking activity\n",
"\n",
"The Otsuka paper predicts a spontaneous firing rate around 10 Hz. However, their corresponding figure shows a slightly smaller firing rate, around 6-7 Hz.\n",
"\n",
""
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
" 30/01/2020 17:18:18: OtsukaSTN: simulation @ A = 0.00 A/m2, 0.0 s stim, 1.0 s offset\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Firing rate: median = 4.82 Hz, last value = 4.86 Hz\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAxAAAAEYCAYAAADMNRC5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XucXVV58PHfkzuBBAIBTLgYIEAAuSYqCHIRELwgBOpLpVZFFEGtoqWtor7aVtG+iohVtFQoVAShilEQ5FaC4WpDgyZIkEsUhKADSUjIPZPn/ePsGc5MZiZ7Zs5lJvl9P5/9mX07e6+zZp2117PX2udEZiJJkiRJZQxpdgIkSZIkDR4GEJIkSZJKM4CQJEmSVJoBhCRJkqTSDCAkSZIklWYAIUmSJKk0AwhJvRIRu0XEjyNiZkTcGxGXRsSYYtuuEXFSH475fC/2/VRE3BERt0XErRExNSK2L9IzMyKWRMSvivmzIuILxfKwqmM8EBGTepnGb0TErr15TafXfyEizulhe3ve9fdcVcc8JiJujIj/LvLjnoj4aER0W/dHxM8i4sZ+nHN88b+ZFRHXRcToTtuHRMR3I+L+Ik2Ti/WHRsSDRZn6fG/TFxHXdz5XN68/tjj3LyPiRz29JiImF3k2KyK+0znfImKL4rMwKyJujojtN3b+MiLio8XfEyPi7Bocb3JEzOtmW6l878M5vxQRsyPi6Ij4z+IzN6UXr98/Io6sVXok1ZYBhKTSImIL4GfA/8vMozPzcOBB4NpilzcBh9fx/PsC7wCOz8w3A/8AXJGZLUV6jgYeBt5TLF9evHQS8On+nDszz8vMp/tzjI1oz7tanCsi3gv8DXBOZr6pyJs3AxOBb3fzml2ArYBxEbF7H0/9f4FrMvONwBzgQ522nwKMyszDgE8BFxXrvwucARwBvD4iDimbvogYBZCZK0qk71LglMw8Engc+EAP+34d+GzxXgI4udP2c4G5xfb/BD5b4vxlfBYgM3+RmZf150AR8dfAD4Hx3eyy0Xzvo9OBYzJzJnBCZh6amfN78frTgH1rlBZJNWYAIak33gbcnZkPtq3IzKuA8RGxB5UG4RkR8Y6I+HBxZ/P+iPgqQERcGREnFvMnRsSV1QePiAsj4ttR8dHizvmsiLgpIkYAfwZ2Bd4fETtl5sPA60qk+/8BfxURB3c637YRcUPnnYu7p/cX6T+vWDczIqYUPQlXR8Qvip6N9xV3+R8r7uZOiogHqo7VobcjIoZGxPei0nsyOyL+OSKGdsq7tnNtU7z3X0bEfRHxpuIYv4mIf42Iu4t9t+6U/lcB76HSiDunSMM1wM3A5cCrImKvLvLpLOCnVBrDH+4iX46IV3p62qZ3dNrtCOAXxfwtwHHdbc/MB4BpETEWGJmZT2bl101vBY7tRfqOBe4q7nbfWvRSzImID0WlF2R+RJxb7Ht0Zv6pmB8GrOriPG2mAneXeS9dbY9Kz1hb78+9EbF/p+17Ff/XuyPizojYKSI+A2wblZ6990XEV4oy9WBUekxmR8T5xWfp1xFxYQ/pB1gMHNXVhjL5HhGnFue+pyj3nXthDi623V3k/a4R8X+BnYGfR8R/UAn4ftrV+y2O8eUif+6PiHcW698HfDIiyny+JTWYAYSk3tgdeLKL9QuAXYCvULn7/DPgTODjxZ3mp6JqCFFXIuJrwLDM/AiVu73bAccVd3eHA6/NzBeo9EAcDtwfEfOBt5dI98vAB4ErI2Jk28rMXJSZp3ax/3uo3JU9EljZxfaVmXkicAPw1sw8qXjvf1kiLbsAD2TmCVQaoOdmZisd867NZ4Hbi7vl7wQuLxpwY4FrM/Mo4FngLZ3O8X+AfyvSswdwGPA1YN/MfBK4HTiw+gXFcc8Avk/ljvXpUelxapeZ97T19FRN1emlSNtLxfwyYOsetgO0FuuWVq3b4HUbSd/bgJ8X8ztTuXt9LpX8+2sq+fOh4j0sLI43HTiGSjDSnSga1mXeS1fbX1dsfwvwsWL/ascDD1EJPL4EjMvMLwGLMrNzALc7lQDq7cA/A58EXl+s61Zm3pSZy7vZvNF8B94FXJyZRwC3dfEe/h34aFEWLwW+npn/BDwPvDkzzyzez8ldvd+IeAuwW9GbeQzwGWA5cGVxrF/19P4kNYcBhKTeeJbKcKDO9gQ6D7k5k8rd77uBV1MJCqpVL+8IHEBleAqZuR5YA1wbEZdTaRQOj8p4+aWZ+f7M3BV4N/CdiNh2YwnPzFnAHcA/bWxfKg3vL1O5I7tNF9v/t/i7BPhtMb8YGNXFvp3f9yLgtRHxA+BiYOSGL2m3D/DLIv3PUmnstY2zn1P8faaL804B5lJpcF9SNIKXA209I5OoBH3VTgDGANcA11O5PpzR4Y2U64FYWhyH4u+SHrZTnKfzuq5e11P6dq4a8jUvM9cWr38yM9fQ6X8TEZ8AzgdOzMyeeiDWbyRNG3uvt1DpwfgplXK3vtP2y4EXqPRifBRY10NansrMl4pz/KkIflcB2cNrNqZMvn8SOLL4HL+BDd/DxKInECpldb8eztfV+90fmBoRM4v1w6nUF5IGMAMISb3xU+D46mEFEfEBoCUzn6LSuGirVz5IZfz9UcDBVBofq4AJxfbqsdZ/otJA3C8qQ5sOoDJO/XQq4/iHUGmIH0AlYGhrDP6Oyh3e1pLp/wzwVmBydzsUPRTvpHLn9U3A+yKic4Omp0bbKmCHYqjSNsBunba/D1iSmX9FZfz/6IgIOuZdm0eBNxbp2gkYB7xYIg0rin0XA1sWx/808KuIOJTK2PKHOr3mA8AHMvPEonfl/wAf6fCmy/VA3Eslj6Fy531Wd9uLtMzNzKXAmojYo0jrCV28rsv0RcSBVIKl9mT2kC8UQ4TeSKV364We9gXmRMTRZd5LN9uPBhYWz+t8Eeg83OhkYFZmHgv8F5VnemDDoBP6Fyh0qWS+nw18ofgcBzC90/bnis8rVIZK/a6HU3b1fucDdxXP6LyJSnDYuS6RNMD0OKRAkqpl5stR+aagiyNiOyp1yG+oNLah0pD7TET8bzH/PxHRQqXn4kEqd8GviIi/olNDIzMzIt5P5a7/UcDyiJgNrAYWUrnT+YOI2Ad4MCJeptLA+LvizmyZ9K+KiDOB+6HyDATwvephTJm5OiIWUXkYezGVYRulH2jOzOcj4nbgf4AniqnancAPI+KNRX48TuXB5uq8a3Mhlfz6C2AL4OzMXFdp6/XoNipDd75IZchPC5VG2bup9E6cS9Xwm4jYgcpwmNOr3se9ETEqIt6QmfeVff/FOa+KiA9Sudt8RnGOtoeMf0IlCL2PSoP0zOJ15wA/AIYCt1U/Z9NT+qgMX7qpTMIiYkfg81R6kG4p8vG6zPxORMwsGrHV/hb496g8f/Mo8KPiOLdRGUr0neK93kOlx+yMTq//NXBdVJ6jaWXD3q/ZwNURsY5Kg/kTxfrfRsTVVHrMSomI9wFk5pUl9n0TcEQx1KjbfC/8Crg9Il6kMsSpc15/EPhWEYCso+chVV293znA0RExi0oP5E8yc1lEPAR8NSIezcy7NvaeJDVWvDK8U5K0qSiGSD0BfDUzX65a/1oq3y70N1VDTzZ7EXFJZn682enoq6IXYFpmXtHstEja9BlASNImqHjo+Ewqd8WHUxkCM4xKz8o/Zuafm5i8AScidsnMZ5qdjr6KiJ2BZ9OLuqQGMICQJEmSVJoPKEmSJEkqzQBCkiRJUmmb7LcwnXjiifmLX/xi4ztKkiRJm4+NfpXfxmyyPRAvvLCxr/eWJEmS1FubbAAhSZIkqfYMICRJkiSVZgAhSZIkqTQDCEmSJEmlGUBIkiRJKs0AQpIkSVJpBhCSJEmSSjOAkCRJklSaAYQkSZKk0gwgJEmSJJVmACFJkiSpNAMISZIkSaUZQEiSJEkqzQBCkiRJUmkGEJIkSZJKM4CQJEmSVNqwZiegOxExB3ipWFwA/BtwCbAOuC0z/7FZaZMkSZI2VwMygIiIUQCZeXTVuoeB04CngJ9HxCGZ+b/NSaEkSZK0eRqoQ5gOBEZHxG0R8d8RcSQwMjOfzMwEbgWO7fyiiDg7ImZHxOyWlpZGp1mSJEna5A3UAGIF8DXgBOAc4D+KdW2WAVt3flFmXpaZ0zJz2vbbb9+QhEqSJEmbkwE5hAn4HfBE0dvwu4h4Cdi2avsYYElTUiZJkiRtxgZqD8T7gYsAImIiMBpYHhF7RERQ6ZmY1cT0SZIkSZulgdoDcTlwZUTcAySVgGI98ANgKJVvYXqwiemTJEmSNksDMoDIzDXAGV1sOrTRaZEkSZL0ioE6hEmSJEnSAGQAIUmSJKk0AwhJkiRJpRlASJIkSSrNAEKSJElSaQYQkiRJkkozgBjkLrjgAj7xiU80OxmbjO9///tcc801zU6GNiGPP/44X/ziF1m+fHmzk6JBYuXKldxyyy2sXr262UnRIDBjxgy+8pWvNDsZ2sxEZjY7DXUxbdq0nD17drOTUXeVH+aGZcuWsdVWWzU5NYPbunXrGD58OACb6udCjTd69GhWrlzJeeedx8UXX9zs5GgQOOOMM7j22ms599xzufTSS5udHA1wbe2AOXPmcNBBBzU5NRokor8HsAdiE7F27dpmJ2HQqw4aDCBUKytXrgTg4YcfbnJKNFhce+21AFx11VVNTokGk0WLFjU7CdqMGEBsIlpbW5udhE2K+alas0ypt9avX9/sJGgQsY5RIxlAbCK80PRfda+DFbFqzc+oest6SL1hHaNGMoDYRDjkpra8cKvW/Iyqtywz6g3LixrJAGITYcXRf/ZAqJ78jKq3LDPqDcuLGskAYhNh12X/VVe+5qdqzTKl3rJBqN6wjlEjGUBsIrzQ1JYVsWrNz6h6yzKj3rC8qJEMIDYRNnj7z69xVT1ZptRblhn1huVFjWQAMYjZ4K0fAzLVmp9RSfXkdUuNNKzZCehKRAwHrgAmASOBLwJ/BG4EHi92+05mXteUBA5AVhz95zMQqicDCEn1ZB2jRhqQAQTwbuDFzPzriNgOmAP8E/D1zLyouUkbOGzw1pb5qXry4q6yhgwZYh2kXrOOUSMN1CFM/wV8rmp5HTAVeFtE/DIiLo+IMZ1fFBFnR8TsiJjd0tLSqLQ2jQ3e+rEiVq35GVVZEdHsJGgQso5RIw3IACIzX87MZUWQ8CPgs8CvgL/LzCOBp4DPd/G6yzJzWmZO23777Rub6Cbzdwv6z4BM9WRQqrIMINQX1jFqpIE6hImI2AX4CXBpZl4TEdtk5pJi80+Af21e6gYGG7z1Y36q1ry4qywDCPWFdYwaaUD2QETEjsBtwD9k5hXF6lsj4nXF/LHAQ01J3ABiAFFb5qfqyTKlsgwg1BfWMWqkgdoDcQEwDvhcRLQ9C/FJ4BsRsQZ4Hji7WYkbiKw4+s+vxVU9WaZUlgGE+sI6Ro00IAOIzPw48PEuNr2h0WkZyLxjXj/mp2rNi7vKGjJkQA4O0ABnHaNGspYaxAwgasv8VD1ZplSWPRDqC+sYNZIBxCbCiqO2zE/VmncHVZYBhPrCOkaNZAAxiHnHvLZ8BkLSQGAAob7wuqVGMoDYRBhA9J8BmerJi7vKMoBQX3jdUiMZQAxiNnjrx/xUrVmmVJYBhPrCmxRqJAOIQay6svCXqPvPgEz15MVdZfktTOoL6xg1krXUJsIGb21ZEavWLFMqyx4I9YXtADWSAcQg5h3z2jI/VU+WKZVlAKG+8CaFGskAYhCzwVtb5qfqyYu7ynIIk/rCOkaNZC21ibDBW1vmp2rNi7vKsgdCfeF1S41kADGIece8tvwdCNWTn1GVZQChvvC6pUYygBjEDCBqy/xUPXlxl1RP1jFqJAOITYQN3toyP1VrXtxVls9AqC+8bqmRrKUGMe+Y15b5qXqyTKkshzCpL7xJoUYygNhE2DipLSti1ZplSmUZQKgvrGPUSAYQg5h3zGvL/FQ9eXFXWQ5hUl943VIjWUsNYtUNktbW1iamZNNgAKF6MoBQWfZAqC+sY9RIgyaAiIghEfHdiLg/ImZGxORmp2kgscFbW+anas2Lu8oygFBfWMeokQZNAAGcAozKzMOATwEXNTk9Tecd89rydyBUT5YplWUAob6wjlEjDaYA4gjgFwCZ+QAwrbnJaT4DiPoxP1VrlimVZQChvrCOUSMNpgBiLPBS1XJrRAyr3iEizo6I2RExu6WlpbGpazIrjv4zIFM9eXdQZRlASBroBlMAsRQYU7U8JDPXVe+QmZdl5rTMnLb99ts3NnVNYIO3tsxP1ZMBhMryW5jUF9YxaqTBVEvdC7wVICIOBeY2NznNZ4O3fqyIVWt+RlWWPRDqC+sYNdKwje8yYPwEOD4i7gMCOLPJ6RlQrDj6z4BM9WRQqrIMINQX1jFqpEETQGTmeuCcZqdjILHBWz/mp2rNi7vKMoBQX1jHqJEG0xAm9cAGb/8ZkKmeLFMqywBCfWEdo0YygBjE/CXq2vJ3IFRPlimVZQChvrCOUSMZQAxi3jGvH/NTUrP4LUzqCwMINZK11CbCBm//GZCpnry4qyx7INQXXrfUSAYQg5gN3voxP1VrlimVZQChshx6q2YxgBjEDCBqy4pY9WSZUlkOYVJf2A5QI1lLbSKsOPrPgEz1ZAChsuyBUFlet9QsBhCDmBVH/ZifqjXLlMoygFBZtgPULAYQg5gVR22Zn6oneyAk1ZrXLTWLAcQmwoqj/3wGQtJA4DMQ6gvbAWoka6lBzDsP9WN+SmoWhzCpLNsBahYDiE2EFUf/WRGrHmwMqrcsMyrL65aaxQBiEKuuOFpbW5uYkk2PFbFqxeEo6i3LjMoygFCzWEsNYtWVhRVH//kMhOrBu8nqLcuMyjKAULMYQAxiVhy1ZX6qHmwMqrcsMyrL65aaxQBiELPiqC3zU/VgY1C9ZZlRWV631Cy9DiAiYsuIGFqPxKh3rDhqy/xUPdgYVG9ZZlSW1y01y0YDiIgYEhFnRMTPI+LPwHxgYUQ8EhFfjYg9a5mgiNg6Im6MiLsj4v6IOKxYf2pEPBkRM4vpqFqedzCy4qgtn4FQPdgYVG9ZZlSW7QA1y7AS+9wF3AF8GpiXmesBImJb4BjgKxHxk8y8ukZp+iRwZ2Z+IyL2Bq4FDimmv8/MH9foPIOeFUdtmZ+qB79RR71lmVFZXrfULGUCiOMyc23nlZm5CPgx8OOIGF7DNF0MrK5K36pifipwcEScB/wK+IfMXFfD8w46Vhy1ZX6qHrybrN6yzKgsr1tqljK3OS6OiDf0tENXAUYZEXFWRMyrnoA9M3NlRLwKuJpKzwfA7cDfAEcCWwHndHG8syNidkTMbmlp6UuSBhUrjtoyP1UPNgbVW5YZleV1S81SJoB4HLgoIn4fEf8SEQfV6uSZeXlmvqbT9D8RsT9wJ3BBZt5d7H5FZj6VlU/LT4GDuzjeZZk5LTOnbb/99rVK5oBlxVFb1XnoMxCqFYejqLcMIFSW7QA1y0avbJl5SWYeBhwFLAL+IyIejYj/GxF71TpBEbEv8F/AGZl5S7EugN9ExM7FbscCD9X63IONv0RdW1bEqgcbg5IaweuWGqn0rbHM/ENm/ktmHgycAUwHHq1Dmr4MjAIuKb5t6adFr8MHgBsi4m5gNPDvdTj3oGKDt7bMT0nSYOJ1S81S5iFqAIoHpU8E/pJKD8DdwD/WOkGZeXI3628Dbqv1+QYzK47aMj9VD/ZAqLcsMyrL65aaZaMBREQcD7wLeDvwIPBD4OzMXF7ntGkjrDhqy9+BkCQNJrYD1CxleiAuAK4Bzi++ulUDhBVHbZmfkgYCeyBUltctNctGA4jMPAYgIqZFxPeAScXrorI5D6hrCtUtK47aMj8lSYOJ1y01S+lnIIAfAH8HzAUspQOAFUdtVeeh+SlpIMhMeyTULdsBapbeBBAtmfmzuqVEvWbFUVs+A6F6sPGn3ur8mzSWIXXHdoCapTcBxOeLIUx3AqvbVmbmDTVPlUqx4qgt81PSQNC5LvLHCNUdr1tqlt4EEGcCU4DhvDKEKQEDiCax4qgt81P14N1j9ZZ1kcqyrKhZehNAHJiZ+9ctJeo1K47aMj8lDQSdhzBJ3fG6pWbpTb/oAxGxb91Sol6rrjhaW1ubmJJNg89AqB7sgVBv2ShUWV631Cy96YE4AnhvRCyg8gyEX+PaZF5kasv8lDQQ+I1w6gvLihqpNwHEiXVLhfrEBm9tedFWPdgDod6ybldZlhU1y0YDiIiIrPjDxvapbdK0MVYctWV+ShoIvJmhsiwrapYyz0DcFRF/ExG7Vq+MiBER8aaIuAp4b32Sp55UP/dgxdF/jiVVPdgDod7yZobKWrt2bfu8ZUWNVGYI04nA+4FrI2I3YAkwChgK3AZcnJkP1y+J6o4VR22Zn5IGAu8qqyyvW2qWjQYQmbkKuBS4NCKGA+OBlZm5pN6JU8/WrVvXPm/F0X/mp+rBb0hTb1U3Cu0NVU8MINQsvXmImsxcCyysU1rUSzZ4a8uKWPXQuTHokCZtzJo1a9rnrYvUE69bapbe/A6EBhgrjtryrp/qoboxaLlSGatXr26ft25XT7yRqGYZcAFEVDwbETOL6cvF+pMi4n8i4v6I+GCz0zkQGEDUlvmpWstMG4PqNXsgVJbXLTVLr4YwNcgewP9m5kltK4pnLy4GXgssB+6NiBsz8/kmpXFAWLFiRfu846z7b+XKle3zVsSqheq7g2C5UjkGnSrLYFPN0qcAIiKGAfsBIwEy81c1TNNUYKeIuAtYCXyiOM8Tmbm4OP89wBuB/6rheQedZcuWtc9bcfTf0qVL2+fNT9VCdZAPliuVs3z58vZ5y4x64nVLzdLXIUzXA2+h0og/oq8nj4izImJe9QQ8D3w5M48BLgSuBsYCL1W9dBmwdRfHOzsiZkfE7JaWlr4ma9BYtGhR+7wVR/8tXry4fd6x6qqFzvWQn1NtzIoVK+wNVWkvvvhi+7xlRY3U1yFMj2TmV/p78sy8HLi8el1EjAbWFdvviYidqAQMY6p2G0Pl9yg6H+8y4DKAadOmbfItwPnz57fPW3H0n/mpWnv00Uc7LFuutDHV9RB4M0M987qlZulrD8TaiLg9Iq6JiGtqmiL4PHAeQEQcCDwN/BbYMyK2jYgRwJHA/TU+76CyePFibrrppvZlK47+WbhwITfffHP7ss+UqBauvPLKDst+TrUxl1/e4Z6aZUbdWr58Odddd137smVFjdTXHohXZebxNU3JK74CXB0Rb6PSE/G+zFwbEZ8EbqUS9FyRmc/W6fwD3h//+EdOP/10Fi9ezBZbbMHKlSutOPph7ty5nHHGGbz88svt66q/2ULqrZdffpnPfOYz3HDDDYwePbr9WQg/p+rOihUr+NKXvsSll17K0KFD229iWGbUlaeffpr3vve9PP300wwbNox169ZZVtRQfe2BGB0RfxkRb42It9YyQZm5ODPflplHZeaxmTm/WH9jZr42M6dm5rdrec7BYNWqVdxxxx186EMfYq+99uK+++5j55135qqrrgK8yPTGunXrmD9/Pt/5znc48cQTOeCAA5g3bx577703//qv/9q+j9Qbixcv5tZbb+VjH/sYu+22G9/85jcZNmwYV155JePGjQP8nOoVmckf//hHfvrTn/KRj3yESZMmceGFFxIRfPvb32bvvfcGLDOqyEz+8Ic/cPXVV3P66aez5557MnPmTHbccUeuv/56wLKi7rW0tHD77bdz/vnnM2XKlJocc6M9EBHxLeCazLyvavVdVL4ZafuapEJA5VuVnnvuORYuXMhzzz3HM888wyOPPNI+VX+136mnnsq3vvUtFi6s/DD45lhxrF+/nlWrVrU/dLhy5cr2+RUrVrBkyRL+9Kc/dZieeOIJfve733X46rtRo0Zx1llnceGFF/Lf//3fgD0Qqlyw165dy4oVK1i8eDGLFi1i0aJF7fMtLS38/ve/5/e//z0LFixgwYIFHV7/ute9jksuuYRDDz2Uc889F9g8P6ebsnXr1rFmzRpWr17dYVq6dCkvvfTSBtOf//xnnnnmGZ555hmefvrpDl+EATB16lS+9rWvcfTRR3PJJZcAlplN1fr16ztct5YvX86LL77YPr3wwgu8+OKLPP300zzxxBM88cQTLFnyyqOfEcHpp5/ORRddZA/nZqS1tZU1a9awfPlyli1bxtKlS1m2bFmH+UWLFvHcc8+1T0899VR7W7GWygxhehy4KCImANcB12bmVTVPSY29+OKLXHnllWRm+0NobfNdLZfZpz+vWbduXYeLSucLzJIlSzb4ysfODj74YN785jfznve8h3333ReAP/3pTwDMmTOHCy64gL/9279lu+22q1k+9lZ1o6urafny5Rs09Lv6W2afVatW9Tmdu+66K69//es5/vjjmT59OuPHjwdg+PDhANx77718/etf5wMf+ABjx46tSd7UQ2bS2tpKa2sr69ata59vbW3doAzW62/btH79+vap83LZdd3t09raytq1a1m3bh3r1q1rn+/r39WrV7eXs56m3lyQR44cydSpUzn88MM5/fTTOeSQQ4gIgPa/n/vc5/jgBz/IIYcc0uv/db211VNr167d6NTa2tr+/9nYfH/2rV7f3VTv7WvXru0ySFi9enW/G2zbbLMN06ZN49BDD2X69OkcfPDBG5SZCy+8kHPOOYfDDz+83//jRmorT21TT3ncl/9LPV5bj3OuW7dug+vXihUr+nT9GjduHIcffjhHHXUUp59+OrvssgsATz75JABPPPEEX/jCF/jwhz/MDjvsUNP/Z61Ul4vqOqV6uboeKDP1dv+u6pzO/8fu6qBav2bNmjUdptWrV2+wrnp9X+ucLbfckv3335/DDz+ck046aeMvKCHKfsNDRLwa+MtiGgVcA/wwMx+vSUpqLCIG3VdXbLHFFkycOJGJEycyYcIEdtppJ6ZMmcJ+++3HfvvtxzbbbLPBa37zm99w4IEHti8fc8wx3HHHHQwZ0r8fGV++fDnPP/88Cxcu5Pnnn+f5559n8eLFLFmyZIPppZdeYvny5SxfvpwVK1a0j91thFGjRjF69Gi22GKLDf6OHTuWHXfcscM0adIkpkyZwlZbbdXl8W677TZOOOGE9uVTTjnFYLdqAAAWaElEQVSFG264of1i3heZydKlS3n22Wd57rnnWLRoUZf52NZoXbVq1Qbzq1ataq9Yq4MF7zjV17Bhwxg9ejTjxo1j3LhxbLvttmy77baMGzeO7bbbjle/+tXstttu7LbbbkyaNIkRI0Z0eZwdd9yRP//5zwCMHTuWRx55hJ133rlfactMlixZwrPPPktLSwuLFy/eYFq2bFl7g6W68VK9vGbNmvaypd4ZMmQII0eObJ9GjBjByJEjGTNmDFtvvfUG0/jx49lll13apx133LHbumX//fdn3rx5QCU4feihh9hvv/36neZVq1bx3HPP0dLS0mU9tGzZsg5BdHV9tHLlSlavXt0hMOhusm7auLbrVdu07bbbst122zF+/Hi22247tttuOyZOnMiee+7JHnvswQ477NBleVmwYAG77757+/K0adN44IEHGDp0aL/St379elpaWnjuued48cUXN6hflixZ0uNNwRUrVrB69eoNbuCo7yKCESNGMHr0aMaOHcuYMWMYM2ZMh/lx48a1tyUnTpzILrvswqRJkzq3C/veqCmUfog6M/8A/AvwLxFxMHAF8AWgfyW0Trbddlve/va3ExHtE9Djcpl9+vqaoUOHMnbs2C4vKmPHjmWbbbZhq6226nVDtXOgcNddd3HZZZdxzjnnbPS1q1evZu7cufz6179u7yJ9/PHHeeqppzr8SF1vDR8+vEOl2DZtueWWXTb0u/u7sX1GjhzZ70Cps2HDOn4kZsyYwYwZM5g+ffpGX7to0SLmzJnD/Pnzeeyxx3jsscdYsGABzz777EZ7l/pjyJAhDBs2jKFDhzJ06FCGDRvGkCFD2vOm8x3NevxtO9+QIUM2WO5qXZl9Oq8bPnw4w4YN6/XfrtaNGDGivYx1NbWVs87loa+qP9dLly7lU5/6FFdfffVGX5eZLFiwgEcffZT58+czf/58Hn/8cZ599lmeffbZDr8XUAttedTT1JanQ4YMYejQoe3/n3rNty23zXc19Wf7xl47bNiwDkFC9VSr8tGV6rpt9erVfPzjH+f2228vdY147rnn+O1vf8tjjz3G/PnzeeKJJ9rLTOdhU/XUVlZ6yt++/l9q/bp6nrPztXCLLbZg1KhRNbt+dT7O7NmzueKKK/jgBz+40deuXbu2Q/3y2GOP8dRTT/Hss8+ycOHCujT4hw4d2qE+6TxffQ0rO1XXG719Xef/YfXxyqzrz2vabjx0NXW1bejQof26oVlLvemBGA6cSKUH4ljgbirDmWbUL3l9N23atJw9e3azk1F3v/3tb9vvSo0ZM4Zly5YxZswYHnnkkfbuzTYrVqzg7rvv5tZbb+WXv/wl8+bN63as/8iRI3nVq17FhAkTmDBhAq961avYbrvt2Hrrrdlmm206TFtvvTVbbbVVe+XYNgxoMJo1axZHHnlkh3U77bQTjz76KGPGjOmwftmyZdx8883cfPPNPPDAA/zud7/r9rijR49mp512YsKECWy//fbteTdu3Lj2PNxyyy0ZNWpU+8Wl+m9bRdLWoKm+OA2UykTdGzp0aPsd2ZEjR7J69Wruuusujj766A77ZSaPPPIIN910E/feey8PPPAAL7zwQrfH3Wqrrdhpp53Ycccd23tJqqexY8d2CNyrA/DRo0czatQoRowY0X4BtywNHIcccghz5swBXikz119/Pe985zs32HfBggXceOONzJo1iwceeIA//vGP3R532LBhTJgwob3MVNfl48aNY6uttuoQTLfVQW3zo0aN6hBIdje13QBQ/T3zzDPsuuuuHdaNHz+eefPmseOOO3ZYv27dOu69915uvvlm7r//fmbPnt3jjYi2XpC261Z1/dJ247Nz3VJdx4wcOXKD4MByMSDUvwciIo4H3gW8DfgV8EPg7Mxc3t+Tq/+q7zycf/75zJkzhxkzZnDmmWdyww030Nrayq233sr111/PLbfc0mHcZUQwZcoUDj74YKZMmcLkyZOZPHkye+yxB9tuu+1m+SGvDn4+97nPccsttzB79mze//73c8UVV7B06VJuvvlmZsyYwR133LHBw9gHH3ww++67L3vvvTdTpkxhjz32YOedd2bMmDGbZX6qoi14GD9+PB/96Ef5whe+wLvf/W7uuOMO9txzTx588EFuuOEGZsyY0T6euc0OO+zAAQccwJQpU5gyZQp77rknu+yyCzvttNOAfj5H/VPdC3zxxRfz4Q9/mLPPPpvJkydz0EEHMXfuXH70ox/xk5/8pH2oU5uxY8e2l5m9996bvfbaq73MjB8/vuY9t2qu6v/neeedx9y5c7nzzjt517vexQ033MCIESO4/fbbmTFjBjfeeGOHX68G2H333XnNa17TXse0XbcmTpzIqFGjGv12NEhstAciIu6i8rzDjzOzcX2f/bS59EA8/vjj7LXXXkDlR6uOP/54XvOa17B48WIiYoNfMZ06dSonnHACxx13HNOmTdvgrvrm7qGHHmLatGkAXHXVVRx00EEcccQRXQ7piggOP/xwTj75ZI455hgOOOCAQd37ovppCx732Wcffv3rX3PMMcdw7733AjBixIgOgej48eM56aSTOP744znssMN49atfbfC5GRo3blz7t+6sX7+eU045hZ/97GcA7b//02bs2LG85S1v4c1vfjOHHXYYe++9t0HCZmThwoVMnDgRgG984xuceuqpTJ06lZaWFoYOHdr+5RRtJk+ezCmnnMLRRx/N61//+vYvEdFmpf49EJl5TH9PovqpvkhMmjSJiRMncvfdd/ORj3yEWbNmMXLkSF73utdx2mmncdppp/X7wc1NXfWY5t12240DDjiAmTNn8vGPf5z77ruPLbfckqOOOop3vOMdvOMd79ige1jqydZbb83w4cO55ZZb+NjHPsYPfvAD1qxZw+67787JJ5/M9OnTecMb3sDQoQPy0TI1UOev7Lzuuus4//zz+d73vsfKlSuZMGECJ510EqeddhpHH310tw/wa9O35ZZbts/vuuuu7LLLLtxzzz186EMfYubMmQwZMoSpU6cyffp0pk+fzj777ONNCfVb6WcgBpvNpQeieuzjggULmDRpUvu2tWvXtj+wo3KefPJJJk+eDFTytjrgam1tdVyv+qStzBx//PHcdttt7evbvlK2ugEgwStlZsqUKTz66KPt69euXcvy5cu7/FY+bZ4ys/1m4qxZszjiiCPat61cuZKhQ4caYKqzfjdk7OMc5Noeyp0wYcIGD00PHz7c4KGXdt11V0aPHs3WW2/NTjvt1GGbDyyrr8477zwAPvaxj3VYP2zYMIMHdenzn/88ABdccEGH9cOHDzd4UAcR0f517vvss0+HbVtssYXBg+rCHohNwPPPP8/QoUPZfnt/GLwWXnzxRYYMGcK4ceOanRRtIlatWsXChQvZbbfdmp0UDRKtra3MnTuXAw880BsX2qjFixezatUqJkyY0OykaHDod6ViACFJkiRtPhzCJEmSJKlxDCAkSZIklWYAIUmSJKk0AwhJkiRJpRlASJIkSSpto79E3WgR8SngxGJxG+BVmfmqiPgkcBbQUmz7UGY+1ow0SpIkSZurARdAZOZXgK8ARMRNwD8Umw4B3pOZDzUrbZIkSdLmbsAOYYqIU4HFmXlrsWoq8OmIuCciPt3EpEmSJEmbraYGEBFxVkTM6zS9ttj8aeAfq3b/IXAO8CbgiIh4exfHOzsiZkfE7JaWls6bJUmSJPXTgPwl6ojYF7gkM48vlgMYm5kvFcsfBrbLzH/u7hj+ErUkSZK0gU32l6iPA26pWh4LzIuIrYpg4k2Az0JIkiRJDTZQA4i9gafaFoqehwuAu4BZwCOZeXOT0iZJkiRttgbctzABZOZHulj3feD7TUiOJEmSpMJA7YGQJEmSNAAZQEiSJEkqzQBCkiRJUmkGEJIkSZJKM4CQJEmSVJoBhCRJkqTSDCAkSZIklWYAIUmSJKk0AwhJkiRJpRlASJIkSSrNAEKSJElSaQYQkiRJkkozgJAkSZJUmgGEJEmSpNIMICRJkiSVZgAhSZIkqTQDCEmSJEmlNT2AiIjpEXFN1fKhEfFgRNwbEZ8v1g2JiO9GxP0RMTMiJjcvxZIkSdLma1gzTx4RlwAnAA9Xrf4ucBrwFPDziDgEmASMyszDIuJQ4CLg5AYnV5IkSdrsNbsH4j7g3LaFiBgLjMzMJzMzgVuBY4EjgF8AZOYDwLQmpFWSJEna7DUkgIiIsyJiXqfptZl5HZBVu44FllYtLwO2Lta/VLW+NSI26D2JiLMjYnZEzG5paanHW5EkSZI2aw0ZwpSZlwOXl9h1KTCmankMsAQY3Wn9kMxc18V5LgMuA5g2bVp23i5JkiSpf5o9hKmDzFwKrImIPSIiqDwfMQu4F3grVB6yBuY2L5WSJEnS5qupD1F34xzgB8BQ4LbMfDAi/gc4PiLuAwI4s5kJlCRJkjZXUXlWedMzbdq0nD17drOTIUmSJA0k0d8DDKghTJIkSZIGNgMISZIkSaUZQEiSJEkqzQBCkiRJUmkGEJIkSZJKM4CQJEmSVJoBhCRJkqTSDCAkSZIklWYAIUmSJKk0AwhJkiRJpRlASJIkSSrNAEKSJElSaQYQkiRJkkozgJAkSZJUmgGEJEmSpNIMICRJkiSVZgAhSZIkqbSmBxARMT0irqlaPjYi7o+IX0bEjyJidLH+ZxFxb0TMjIhbmpdiSZIkafM1rJknj4hLgBOAh6tWXwocmZl/iogvAx8AvglMBvbLzGx8SiVJkiRB83sg7gPO7bTu6Mz8UzE/DFgVETsC2wA3RsQ9EfH2RiZSkiRJUkVDAoiIOCsi5nWaXpuZ1wEdehQyc2HxmunAMcB/AiOAi4BTgFOBiyNihy7Oc3ZEzI6I2S0tLfV+W5IkSdJmpyFDmDLzcuDysvtHxCeAvwBOzMxVEfE88N3MXAf8OSLmAHsDf+50nsuAywCmTZvmUCdJkiSpxpo9hGkDEfEZ4I3AcZn5QrH6OOD6YvtWwGuAR5uTQkmSJGnz1dSHqDsrnnX4PPC/wC0RAXBdZn4nIk6IiAeA9cAFVcGFJEmSpAZpegCRmTOBmcX8n6g879DVfuc1LlWSJEmSujLghjBJkiRJGrgMICRJkiSVZgAhSZIkqTQDCEmSJEmlGUBIkiRJKs0AQpIkSVJpBhCSJEmSSjOAkCRJklSaAYQkSZKk0gwgJEmSJJVmACFJkiSpNAMISZIkSaUZQEiSJEkqzQBCkiRJUmkGEJIkSZJKM4CQJEmSVJoBhCRJkqTSmh5ARMT0iLimavnUiHgyImYW01ERMSQivhsR9xfrJjczzZIkSdLmalgzTx4RlwAnAA9XrT4E+PvM/HHVfqcCozLzsIg4FLgIOLmhiZUkSZLU9B6I+4BzO62bCrw/ImZFxEURMQw4AvgFQGY+AExrbDIlSZIkQYN6ICLiLOATnVafmZnXRcTRndbfDswAFgDfBc4BxgIvVe3TGhHDMnNdp/OcDZxdLL4cEY/V6C30x3jghWYnQgOW5UM9sXyoO5YN9cTyoZ7My8zX9OcADQkgMvNy4PKSu1+RmUsAIuKnwGlUgocxVfsM6Rw8FOe5DLisn8mtqYiYnZn2mKhLlg/1xPKh7lg21BPLh3oSEbP7e4xmD2HqICIC+E1E7FysOhZ4CLgXeGuxz6HA3OakUJIkSdq8NfUh6s4yMyPiA8ANEbES+C3w70ArcHxE3AcEcGYTkylJkiRttpoeQGTmTGBm1fJtwG1d7HpOg5JUawNqSJUGHMuHemL5UHcsG+qJ5UM96Xf5iMysRUIkSZIkbQYG1DMQkiRJkgY2AwhJkiRJpRlA1ElEDImI70bE/RExMyImNztNaryIGB4R3y9+GPFXEfGOiJgcEfcU674TEUOKfT9f7HNfRLyu2WlXY0TEDhHxTERMsWyoWkR8uriGPBQRZ1k+1Ka4tlxT/M9nWX+oTUS8PiJmFvOly0R3+3bHAKJ+TgFGZeZhwKeAi5qcHjXHu4EXM/ONwFuAbwFfBz5brAvg5Ig4BDgKeD3wl8C3m5ReNVBEDAf+DVhZrLJsCIDiR1bfABxO5f+/C5YPveKtwLDMfAPwT8CXsHxs9iLi74HvAaOKVb0pExvs29O5DCDq5wjgFwCZ+QDgD7psnv4L+FzV8jpgKnB3sXwLcByV8nJbVjwNDIuI7RuaUjXD14DvAs8Vy5YNtTmBym8e/QS4EbgJy4de8Tsq/+shwFhgLZYPwZPAqVXLvSkTXe3bLQOI+hlL5Re027RGRNO/NleNlZkvZ+ayiBgD/Aj4LJVvP2v7+rNlwNZsWF7a1msTFRHvA1oy89bq1ZYNFcZTufH0TipfY/4DYIjlQ4WXgUnAfCq/l/VNrD82e5n5YyrBZJvelImu9u2WAUT9LAXGVC0Pycx1zUqMmicidgHuAr6fmdcA66s2jwGWsGF5aVuvTdf7qfxA5kzgIOA/gR2qtls2Nm8vArdm5prMfAxYRccLuuVj8/YJKuVjL+BA4CpgRNV2y4egd+2NrvbtlgFE/dxLZYwiEXEola5obWYiYkcqP4z4D5l5RbF6TjG+GSrPRcyiUl5OKB6+35VKwPlCwxOshsnMIzPzqMw8GngYeA9wi2VDhXuAE6NiIrAlcKflQ4XFvHIXeREwHK8t2lBvykRX+3bLITX18xMqdxfvo/IwyplNTo+a4wJgHPC5iGh7FuLjwDcjYgTwKPCjzGyNiFnA/VQC+480JbVqtr8F/t2yocy8KSKOBH7FK//3BVg+VHExcEXxvx9B5VozG8uHOurNNWWDfXs6sL9ELUmSJKk0hzBJkiRJKs0AQpIkSVJpBhCSJEmSSjOAkCRJklSaAYQkSZKk0gwgJEmSJJVmACFJkiSpNAMISRIRsU1EfLhq+b46nWeLiLg7Iob28zgjIuKXEeEPokpSgxlASJIAtgHaA4jMfEOdzvN+4IbMbO3PQTJzDXAncHpNUiVJKs0AQpIE8BVgj4h4OCK+GhEvA0TEpIiYHxHfi4h5EfGDiDguIu6NiMcj4nVtB4iId0fEr4pj/Fs3vQx/Bfy0N8eOiC0j4ucR8etiv7agYUZxPElSAxlASJIAPgU8mZkHZebfddo2GbgEOACYApwBHAGcD1wAEBH7UOkNODwzDwJa6dS4j4gRwO6Z+fveHBs4EXguMw/MzNcAvyjWzwNe27+3LUnqLQMISdLGLMjMuZm5HngEuDMzE5gLTCr2ORaYCvxPRDxcLO/e6TjjgSV9OPZc4LiI+JeIeGNmvgRQDINaExFjavheJUkb4cNnkqSNWV01v75qeT2vXEcCuCozP93DcVYCo3p77Mz8XURMBd4KfDkibsvMfyr2Gwms6sV7kST1kz0QkiSAZUB/7uTfCfxFROwAEBHbRsSrq3fIzMXA0IjoHET0KCImAisy82rga8AhxfrtgJbMXNuPdEuSeskeCEkSmfli8fDyPOCWPrz+txHxWeC2iBgCrAU+Avyh0663UXnG4Y5eHH5/4KsRsb447rnF+mOAm3ubVklS/0RlqKkkSfUXEQcDn8zMv67BsW4APp2Zj/U/ZZKkshzCJElqmMycA9xVix+SA2YYPEhS49kDIUmSJKk0eyAkSZIklWYAIUmSJKk0AwhJkiRJpRlASJIkSSrNAEKSJElSaQYQkiRJkkr7//UDtzDjhz9IAAAAAElFTkSuQmCC\n",
"text/plain": [
"