Fig3cd-ion-channel-discrete-model.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import msmtools\n",
    "import scipy.linalg\n",
    "import itertools\n",
    "import pyemma\n",
    "\n",
    "import matplotlib.patches as mpatches\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm import tqdm_notebook as tqdm\n",
    "\n",
    "from local import Channel, cg_transition_matrix, kchannel_params, dndt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic ion channel model\n",
    "\n",
    "$$C1 \\underset{\\beta}{\\overset{4\\alpha}{<=>}} C2  \\underset{2\\beta}{\\overset{3\\alpha}{<=>}} C3 \\underset{3\\beta}{\\overset{2\\alpha}{<=>}} C4 \\underset{4\\beta}{\\overset{\\alpha}{<=>}} O$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Verify that we can reconstruct matrix for uncoupled systems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.84 0.15 0.01 0.   0.  ]\n",
      " [0.   0.87 0.12 0.01 0.  ]\n",
      " [0.   0.01 0.91 0.08 0.  ]\n",
      " [0.   0.   0.01 0.95 0.04]\n",
      " [0.   0.   0.   0.02 0.98]]\n"
     ]
    }
   ],
   "source": [
    "# Instantiate an arbitrary channel\n",
    "c_test = Channel(44., 4.)\n",
    "lagtime = 0.001\n",
    "\n",
    "# Generate macrostate transition matrix\n",
    "full_model = c_test.transition_matrix(lagtime)\n",
    "\n",
    "with np.printoptions(precision=2, suppress=True):\n",
    "    print(full_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate microstate transition matrix and perform mapping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1) Generate close--open transition matrix\n",
    "T_singlet = c_test.transition_matrix_singlet(lagtime)\n",
    "# 2) Expand into 16x16 corresponding to all combinations of C, O for 4 gates.\n",
    "T_allstates = np.kron(np.kron(np.kron(T_singlet, T_singlet), T_singlet), T_singlet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Microstate: 0 Configuration: (0, 0, 0, 0) Macrostate: C1\n",
      "Microstate: 1 Configuration: (0, 0, 0, 1) Macrostate: C2\n",
      "Microstate: 2 Configuration: (0, 0, 1, 0) Macrostate: C2\n",
      "Microstate: 3 Configuration: (0, 0, 1, 1) Macrostate: C3\n",
      "Microstate: 4 Configuration: (0, 1, 0, 0) Macrostate: C2\n",
      "Microstate: 5 Configuration: (0, 1, 0, 1) Macrostate: C3\n",
      "Microstate: 6 Configuration: (0, 1, 1, 0) Macrostate: C3\n",
      "Microstate: 7 Configuration: (0, 1, 1, 1) Macrostate: C4\n",
      "Microstate: 8 Configuration: (1, 0, 0, 0) Macrostate: C2\n",
      "Microstate: 9 Configuration: (1, 0, 0, 1) Macrostate: C3\n",
      "Microstate: 10 Configuration: (1, 0, 1, 0) Macrostate: C3\n",
      "Microstate: 11 Configuration: (1, 0, 1, 1) Macrostate: C4\n",
      "Microstate: 12 Configuration: (1, 1, 0, 0) Macrostate: C3\n",
      "Microstate: 13 Configuration: (1, 1, 0, 1) Macrostate: C4\n",
      "Microstate: 14 Configuration: (1, 1, 1, 0) Macrostate: C4\n",
      "Microstate: 15 Configuration: (1, 1, 1, 1) Macrostate: O\n"
     ]
    }
   ],
   "source": [
    "# Construct mapping from 16 microstates to the 5 macrostates\n",
    "chi = np.zeros((16, 5))\n",
    "for microstate_number, microstate_config in enumerate(itertools.product(*[range(2) for _ in range(4)])):\n",
    "    print('Microstate:', microstate_number, \n",
    "          'Configuration:', microstate_config, \n",
    "          'Macrostate:', c_test.index2state(sum(microstate_config)))\n",
    "    chi[microstate_number, sum(microstate_config)] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Map microstate transitions to macrostate\n",
    "cg_from_full = cg_transition_matrix(T_allstates, chi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Macro from Micro:\n",
      "[[0.84 0.15 0.01 0.   0.  ]\n",
      " [0.   0.87 0.12 0.01 0.  ]\n",
      " [0.   0.01 0.91 0.08 0.  ]\n",
      " [0.   0.   0.01 0.95 0.04]\n",
      " [0.   0.   0.   0.02 0.98]] \n",
      "\n",
      "Original:\n",
      "[[0.84 0.15 0.01 0.   0.  ]\n",
      " [0.   0.87 0.12 0.01 0.  ]\n",
      " [0.   0.01 0.91 0.08 0.  ]\n",
      " [0.   0.   0.01 0.95 0.04]\n",
      " [0.   0.   0.   0.02 0.98]]\n"
     ]
    }
   ],
   "source": [
    "with np.printoptions(precision=2, suppress=True):\n",
    "    print(\"Macro from Micro:\")\n",
    "    print(cg_from_full, '\\n')\n",
    "    print(\"Original:\")\n",
    "    print(full_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Matrices match!\n"
     ]
    }
   ],
   "source": [
    "# Verify that the the transition matrices match up\n",
    "if(np.allclose(cg_from_full, full_model)):\n",
    "    print(\"Matrices match!\")\n",
    "else:\n",
    "    raise AssertionError(\"Matrices do not match!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hodgkin-Huxley Potassium Ion Channel Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "E_k = -75.\n",
    "g_k = 24 #mmho/cm^2\n",
    "\n",
    "V_clamp = 63"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = 0.1\n",
    "T = 15\n",
    "time = np.arange(0, T, dt)\n",
    "\n",
    "# initial condition\n",
    "n_0 = 0\n",
    "\n",
    "# Other parameters of the ODE\n",
    "params = (V_clamp,)\n",
    "\n",
    "ode_solution = odeint(dndt, n_0, time, params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f770083c6d8>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAfQklEQVR4nO3deXhU9b3H8fc3IWxBBSGKohRQiyDXrQFFrldxK1VAlFIWi7iCogg+0ora6r1avYpY5JKIBAihilgLCIg3FSq1VhQwLCKLFYSyRIWglZ0s5Hv/IHIhBggzk5xZPq/nmSczZ4acDyH58Mtvzjk/c3dERCT2JAUdQEREQqMCFxGJUSpwEZEYpQIXEYlRKnARkRhVozp31qhRI2/WrFl17lJEJOYtXrx4m7unld9erQXerFkz8vLyqnOXIiIxz8w2VLRdUygiIjFKBS4iEqNU4CIiMUoFLiISo1TgIiIx6pgFbmbZZrbVzFaU2z7IzP5hZivNbHjVRRQRkYpUZgSeA3Q6dIOZdQRuBM539/OAEZGPJiIiR3PM48Dd/X0za1Zu873As+5eWPaarZGPJiKJwt3Zv38/JSUllJSUHHa/otuRni8tLT14c/fDHkd6W0WPj3br27cv55xzTkS/bqGeyPNj4HIzexrYBwx1948reqGZ9Qf6AzRt2jTE3YlIVXB3CgsL2b17N3v27GHPnj0H7x9p2759+ygsLKSoqIjCwsLjvl9UVPSDAt6/f3/QX4oqd9lll0VNgdcAGgCXAm2BN8yshVewOoS7ZwFZAOnp6Vo9QiRC3J09e/ZQUFDAN998w3fffcf27dsrvJV/bseOHQdLubS09Lj2a2bUqlWLmjVrUqtWrSPeT01NpUGDBj/YnpKSQkpKCjVq1PjBLTk5ucLtlXk+OTmZ5ORkkpKSMDOSkpIOu4W6rbKvOdqtqoRa4JuB6WWFvcjMSoFGQEHEkokkoNLSUrZt20Z+fj6bN2/m66+/pqCggIKCArZt23bw/ve3ffv2HfXz1atXj5NOOungLS0tjbPPPpsTTzyRevXqkZqaSt26dalbt26l7tetW5eaNWtWaSlJ5YVa4DOAq4D3zOzHQE1gW8RSicSp3bt3s27duoO3jRs3kp+ff/D25ZdfUlxc/IM/l5qaSlpaGmlpaTRu3Jg2bdocfJyWlkajRo2oX7/+YWV94oknkpycHMDfUqrLMQvczKYAVwKNzGwz8ASQDWSXHVpYBPSraPpEJBHt3LmTzz77jFWrVvHFF18cVthbtmw57LWpqak0adKEJk2acPnllx+8//3ttNNOIy0tjTp16gT0t5FoVpmjUHof4alfRjiLSEwpLCxkxYoVLF26lJUrV7J69WpWrVrFpk2bDr4mKSmJpk2b0qJFC7p06UKLFi0Ou5188smajpCQVevlZEVi1b59+1i8eDFLlixh6dKlLFmyhJUrV1JSUgJAnTp1aNWqFVdccQWtWrWidevWtGrVihYtWpCSkhJweolXKnCRCmzZsoUPP/yQ+fPnM3/+fBYvXnxwbjotLY2LL76Y66+/nosvvpiLLrqI5s2bk5SkK1NI9VKBiwD/+te/mDdvHnPmzGHevHmsXbsWgFq1apGens6DDz5I+/btadu2LaeffrqmPSQqqMAlIZWUlLBgwQLmzp3LnDlzWLRoEaWlpZxwwgl07NiRAQMG0KFDBy6++GJq1aoVdFyRCqnAJWEUFhYyb948pk2bxsyZM9m2bRtJSUm0a9eO3/zmN1x33XW0a9dOc9YSM1TgEtf27t1Lbm4u06dP56233mLHjh2ccMIJdOnShW7dunHNNdfQoEGDoGOKhEQFLnHH3Zk/fz6TJk3ijTfeYMeOHTRs2JCf//zn3HzzzVxzzTWaFpG4oAKXuLFlyxYmTJhAdnY2X3zxBampqXTv3p2+ffty5ZVXUqOGvt0lvug7WmLa96Ptl156ialTp1JcXMyVV17Jb3/7W7p37069evWCjihSZVTgEpNKSkqYOnUqw4cPZ+nSpZx00kkMHDiQe++9l5YtWwYdT6RaqMAlpuzZs4eJEyfywgsvsH79es4991yysrLo06cPqampQccTqVYqcIkJ+/btY+zYsTzzzDNs3bqV9u3bM3LkSLp06aIzICVh6TtfolpxcTFZWVmcc845DBkyhDZt2vC3v/2N+fPnc+ONN6q8JaFpBC5Ryd2ZNWsWQ4cOZe3atVx66aVMmjSJq666KuhoIlFDwxeJOitWrODaa6+lW7du1KxZk9mzZ/Phhx+qvEXKUYFL1Ni+fTv3338/F1xwAUuWLOF//ud/WLZsGTfccIMuHiVSgWMWuJllm9nWstV3yj831MzczBpVTTxJFDNnzqR169aMGTOGe+65hzVr1jBo0CBdl0TkKCozAs8BOpXfaGZnAtcCGyOcSRLI119/zS9+8Qu6detGo0aNWLBgAZmZmTRs2DDoaCJR75gF7u7vA99W8NRI4NeA1sKUkEydOpXWrVsza9Ysnn76afLy8mjbtm3QsURiRkhz4GbWFch3908q8dr+ZpZnZnkFBQWh7E7izK5du7jjjjvo0aMH55xzDp988gmPPvqopktEjtNxF7iZ1QUeAx6vzOvdPcvd0909PS0t7Xh3J3Fm0aJFXHTRReTk5PDYY4/xwQcf6NR3kRCFMgI/C2gOfGJm/wTOAJaYWeNIBpP44u6MHj2aDh06UFhYyHvvvcfvfvc7jbpFwnDcJ/K4+6fAKd8/LivxdHffFsFcEkf27NnDgAEDePXVV+natSs5OTlaREEkAipzGOEU4COgpZltNrM7qz6WxIv169fToUMHJk+ezFNPPcWbb76p8haJkGOOwN299zGebxaxNBJX/va3v3HzzTdTWlrK7Nmzuf7664OOJBJXdCamVIkpU6Zw3XXXceqpp/Lxxx+rvEWqgApcIsrdee655+jTpw/t27dn/vz5nH322UHHEolLKnCJmJKSEu677z6GDRtGr169eOeddzTfLVKFVOASEcXFxfTp04cxY8bw8MMPM3nyZK38LlLFdD1wCVthYSE9e/Zk5syZjBgxgoceeijoSCIJQQUuYdm7dy/du3cnNzeXjIwM7rvvvqAjiSQMFbiEbPfu3XTr1o13332XrKws7r777qAjiSQUFbiEpLCwkJtuuol58+aRk5PDrbfeGnQkkYSjApfjVlJSwi233MLcuXPJzs5WeYsEREehyHFxdwYMGMC0adMYOXIkt99+e9CRRBKWClwqzd0ZOnQo2dnZPP744wwZMiToSCIJTQUulTZ8+HB+//vfM2jQIP7zP/8z6DgiCU8FLpUydepUhg0bRu/evXnxxRe1SrxIFFCByzEtWrSIvn37ctlll5GdnU1Skr5tRKKBfhLlqDZu3EjXrl057bTTmDFjBrVr1w46koiUqcyCDtlmttXMVhyy7Xkz+8zMlpvZm2ZWv2pjShB27NhB586d2bt3L7Nnz0ZrmopEl8qMwHOATuW2zQXauPv5wOfAIxHOJQErLS3l1ltvZdWqVUydOpXWrVsHHUlEyjlmgbv7+8C35bbNcfeSsocLOLCwscSR4cOHM3PmTF544QWuvfbaoOOISAUiMQd+B5Abgc8jUWLevHk89thj9OzZkwceeCDoOCJyBGEVuJk9BpQAk4/ymv5mlmdmeQUFBeHsTqpBfn4+vXr1omXLlowfP16HC4pEsZAL3Mz6AZ2BW9zdj/Q6d89y93R3T9ebYNGtqKiIHj16sHfvXqZPn069evWCjiQiRxHSxazMrBPwMHCFu++JbCQJyrBhw/joo4/44x//yLnnnht0HBE5hsocRjgF+AhoaWabzexOIAM4AZhrZsvM7OUqzilV7J133mHkyJHcf//9/OIXvwg6johUgh1l9iPi0tPTPS8vr9r2J5VTUFDA+eefT8OGDfn444+pU6dO0JFE5BBmttjd08tv1/XAE5y7c/fdd/Ptt9/yzjvvqLxFYogKPMGNHz/+4PHe559/ftBxROQ46FooCezzzz9nyJAhXHPNNbq2t0gMUoEnqP3799O3b19q165NTk6OrjAoEoM0hZKgRo0axaJFi5gyZQpNmjQJOo6IhEDDrgS0du1afvOb39C1a1d69uwZdBwRCZEKPMF8f9RJSkoKL730kk6VF4lhmkJJMOPHj+e9994jKytLUyciMU4j8ASSn5/P0KFD6dixI3fddVfQcUQkTCrwBDJw4ECKi4sZN26cpk5E4oCmUBLEW2+9xaxZsxg+fDhnnXVW0HFEJAJ0LZQEsHfvXs477zzq1KnDsmXLSElJCTqSiBwHXQslgT333HOsX7+eefPmqbxF4ojmwOPcunXrePbZZ+nVqxcdO3YMOo6IRJAKPM4NHjyYlJQURowYEXQUEYkwTaHEsbfeeovZs2fz/PPP65hvkThUmRV5ss1sq5mtOGTbyWY218zWlH1sULUx5XgVFhYyZMgQWrVqxeDBg4OOIyJVoDJTKDlAp3LbhgHvuvs5wLtljyWKZGRksG7dOkaNGqU3LkXi1DEL3N3fB74tt/lGYFLZ/UlAtwjnkjB88803/O53v+NnP/sZ1157bdBxRKSKhPom5qnu/hVA2cdTjvRCM+tvZnlmlldQUBDi7uR4PPXUU+zYsYPhw4cHHUVEqlCVH4Xi7lnunu7u6WlpaVW9u4S3Zs0aMjMzufPOO2nTpk3QcUSkCoVa4FvM7DSAso9bIxdJwjFs2DBq1arFk08+GXQUEalioRb4LKBf2f1+wMzIxJFwfPDBB0yfPp2HH36Yxo0bBx1HRKrYMa+FYmZTgCuBRsAW4AlgBvAG0BTYCPRw9/JvdP6AroVSddyd9u3bs2nTJj7//HNSU1ODjiQiERLytVDcvfcRnro67FQSMTNmzGDhwoWMHz9e5S2SIHQ1wjiwf/9+LrjgAkpKSlixYgU1augEW5F4oqsRxrE//vGPrFy5ktdff13lLZJANAKPccXFxbRu3Zq6deuydOlSkpJ0fTKReKMReJz6wx/+wNq1a5k5c6bKWyTB6Cc+hhUWFvLkk0/Srl07unTpEnQcEalmGoHHsKysLDZu3MiECRO0SLFIAtIIPEbt2bOHp59+miuuuIKrr9YRnSKJSCPwGJWVlcWWLVv405/+pNG3SILSCDwG7du3j+eff54rr7ySyy+/POg4IhIQjcBj0MSJE/nyyy955ZVXgo4iIgHSCDzGFBUV8eyzz9K+fXutMi+S4DQCjzGvvvoqGzdu5OWXX9bct0iC0wg8hpSUlPDMM8/wk5/8hE6dyi9TKiKJRiPwGPL666/zxRdf8Oabb2r0LSK6FkqsKC0t5bzzziMlJYVly5bptHmRBKJrocS4GTNm8NlnnzFlyhSVt4gAYc6Bm9mDZrbSzFaY2RQzqx2pYPL/3J3nnnuOFi1a0KNHj6DjiEiUCLnAzawJ8ACQ7u5tgGSgV6SCyf/7+9//zqJFixg6dCjJyclBxxGRKBHu7+I1gDpmVgOoC3wZfiQpb/jw4aSlpXHbbbcFHUVEokjIBe7u+cAIDixq/BWw3d3nlH+dmfU3szwzyysoKAg9aYJasWIFb7/9NoMGDaJOnTpBxxGRKBLOFEoD4EagOXA6kGpmvyz/OnfPcvd0d09PS0sLPWmCGjFiBHXr1mXgwIFBRxGRKBPOFMo1wHp3L3D3YmA6cFlkYgnA5s2bmTx5MnfddRcNGzYMOo6IRJlwCnwjcKmZ1bUDZ5VcDayOTCwBePHFF3F3HnzwwaCjiEgUCmcOfCEwFVgCfFr2ubIilCvhfffdd4wdO5aePXvSrFmzoOOISBQK60Qed38CeCJCWeQQWVlZ7Nq1i1/96ldBRxGRKKVT+qJQcXExo0ePpmPHjlx44YVBxxGRKKVT6aPQ9OnT2bx5My+99FLQUUQkimkEHoVGjhzJ2WefzQ033BB0FBGJYhqBR5mPPvqIhQsXMnr0aF20SkSOSg0RZUaOHEn9+vV12ryIHJMKPIps2LCBadOmcffdd1OvXr2g44hIlFOBR5GMjAzMjPvvvz/oKCISA1TgUWLXrl2MGzeO7t2707Rp06DjiEgMUIFHiUmTJrF9+3adNi8ilaYCjwLuTmZmJm3btuXSSy8NOo6IxAgdRhgF/vrXv7J69WpycnKCjiIiMUQj8CiQkZFBw4YN6dmzZ9BRRCSGqMADtnHjRmbOnMldd91F7dpaE1pEKk8FHrCxY8cCcM899wScRERijQo8QIWFhYwbN47OnTvrmt8ictxU4AH605/+REFBgU7cEZGQhFXgZlbfzKaa2WdmttrM2kcqWCLIzMzkxz/+MVdffXXQUUQkBoV7GOEo4M/u/nMzqwnUjUCmhLB48WIWLFjAqFGjdNVBEQlJyAVuZicC/wHcBuDuRUBRZGLFv8zMTFJTU+nXr1/QUUQkRoUz9GsBFAATzWypmY03s9TyLzKz/maWZ2Z5BQUFYewufnzzzTdMmTKFvn37ctJJJwUdR0RiVDgFXgO4GBjj7hcBu4Fh5V/k7lnunu7u6WlpaWHsLn5kZ2ezb98+7rvvvqCjiEgMC6fANwOb3X1h2eOpHCh0OYr9+/czZswYrrjiCtq0aRN0HBGJYSEXuLt/DWwys5Zlm64GVkUkVRzLzc1l/fr1Gn2LSNjCPQplEDC57AiUdcDt4UeKb5mZmZx++ul069Yt6CgiEuPCKnB3XwakRyhL3FuzZg1//vOf+a//+i9SUlKCjiMiMU4HIFejMWPGkJKSQv/+/YOOIiJxQAVeTXbv3k12djbdu3encePGQccRkTigAq8mr732Gtu3b9d1T0QkYlTg1cDdycjI4IILLuCyyy4LOo6IxAktqVYN5s+fz/Lly8nKysLMgo4jInFCI/BqkJGRQf369enTp0/QUUQkjqjAq9hXX33FtGnTuP3220lN/cGlYkREQqYCr2JZWVmUlJQwcODAoKOISJxRgVeh4uJixo4dS6dOnTj77LODjiMicUZvYlahGTNm8NVXXzFu3Ligo4hIHNIIvAplZGTQvHlzOnXqFHQUEYlDKvAq8umnn/L+++8zcOBAkpOTg44jInFIBV5FMjMzqV27NnfccUfQUUQkTqnAq8C//vUvXnnlFfr06cPJJ58cdBwRiVMq8CowceJE9uzZw6BBg4KOIiJxLOwCN7PkskWNZ0ciUKzbv38/mZmZ/Pu//zsXXnhh0HFEJI5FYgQ+GFgdgc8TF3Jzc1m3bp1G3yJS5cIqcDM7A7gBGB+ZOLFv9OjRNGnShJtuuinoKCIS58Idgb8I/BooPdILzKy/meWZWV5BQUGYu4tun332GXPmzOGee+7RkmkiUuVCLnAz6wxsdffFR3udu2e5e7q7p6elpYW6u5iQkZFBzZo1tWSaiFSLcEbgHYCuZvZP4HXgKjN7NSKpYtCOHTuYNGkSPXv25JRTTgk6jogkgJAL3N0fcfcz3L0Z0AuY5+6/jFiyGJOTk8OuXbv05qWIVBsdBx4BpaWlZGRkcMkll9C2bdug44hIgojI1Qjd/T3gvUh8rlg0Z84c1qxZw6uvJuwMkogEQCPwCBg9ejSnnnoqPXr0CDqKiCQQFXiY1q5dS25uLgMGDKBmzZpBxxGRBKICD1NmZibJyckMGDAg6CgikmBU4GHYvn07EyZMoEePHpx++ulBxxGRBKMCD8O4cePYuXMnDz30UNBRRCQBqcBDVFRUxIsvvkjHjh35yU9+EnQcEUlAWtQ4RG+88Qb5+flkZWUFHUVEEpRG4CFwd55//nlat26tBYtFJDAagYfgL3/5C8uXL2fChAkkJen/QBEJhtonBCNGjKBx48bccsstQUcRkQSmAj9Oy5cvZ86cOTzwwAPUqlUr6DgiksBU4MfphRdeIDU1VSfuiEjgVODHYf369bz22mvceeednHzyyUHHEZEEpwI/Dv/93/9NUlISv/71r4OOIiKiAq+sDRs2MHHiRO6++26aNGkSdBwRERV4ZT377LOYGQ8//HDQUUREgPAWNT7TzP5qZqvNbKWZDY5ksGiyadMmJkyYwJ133smZZ54ZdBwRESC8E3lKgIfcfYmZnQAsNrO57r4qQtmixnPPPQfAsGHDAk4iIvL/wlnU+Ct3X1J2fyewGoi7yeH8/HzGjRvHbbfdxo9+9KOg44iIHBSROXAzawZcBCys4Ln+ZpZnZnkFBQWR2F21Gj58OKWlpTz66KNBRxEROUzYBW5m9YBpwBB331H+eXfPcvd0d09PS0sLd3fVauPGjYwdO5Zbb72VZs2aBR1HROQwYRW4maVwoLwnu/v0yESKHo8//jgATzzxRMBJRER+KJyjUAyYAKx2999HLlJ0WL58OX/4wx944IEHaNq0adBxRER+IJwReAegL3CVmS0ru10foVyBGzZsGPXr1+eRRx4JOoqISIVCPozQ3T8ALIJZokZubi65ubmMGDGCBg0aBB1HRKRC5u7VtrP09HTPy8urtv2FoqioiDZt2mBmfPrpp9SsWTPoSCKS4Mxssbunl9+uFXnKGTVqFGvWrCE3N1flLSJRTddCOcSmTZt48skn6dKli9a6FJGopwIv4+7ce++9lJaWMmrUqKDjiIgck6ZQykyZMoW3336bkSNH0rx586DjiIgck0bgwNdff83gwYO55JJLGDRoUNBxREQqJeELvLS0lH79+rFr1y6ys7NJTk4OOpKISKUk/BTKyJEjmTNnDi+//DKtW7cOOo6ISKUl9Ah8/vz5PPLII9x88830798/6DgiIsclYQt8w4YN3HTTTTRr1ozx48dz4NIuIiKxIyELfPv27XTt2pWioiLeeustnS4vIjEp4ebAd+3axfXXX8+qVauYPXs2LVu2DDqSiEhIEmoEvnPnTrp06cLChQt5/fXX+elPfxp0JBGRkCXMCHzTpk107tyZlStX8sorr9C9e/egI4mIhCUhRuBz586lXbt2/POf/+Ttt9+md+/eQUcSEQlbXBf4tm3buOeee7juuuuoX78+H374oaZNRCRuhLsmZicz+4eZrTWzYZEKFa4NGzbw6KOP0rx5c8aNG8evfvUrli5dynnnnRd0NBGRiAl5DtzMkoFM4FpgM/Cxmc1y91WRClcZxcXF5Ofn849//IMFCxYwd+5c5s+fj5nRs2dPfvvb3+oMSxGJS+G8idkOWOvu6wDM7HXgRiDiBf7UU0/x2muv4e6UlpZSWlpKUVERO3fuZMeOHZSWlgKQlJTE+eefz9NPP03v3r11VUERiWvhFHgTYNMhjzcDl5R/kZn1B/oDIa/uftppp/Fv//ZvJCUlHbzVqFGDE088kfr169O0aVPOOuss0tPTOeGEE0Lah4hIrAl5TUwz6wH81N3vKnvcF2jn7ke8HmssrIkpIhJtjrQmZjhvYm4Gzjzk8RnAl2F8PhEROQ7hFPjHwDlm1tzMagK9gFmRiSUiIscS8hy4u5eY2f3AO0AykO3uKyOWTEREjiqsU+nd/X+B/41QFhEROQ5xfSamiEg8U4GLiMQoFbiISIxSgYuIxKiQT+QJaWdmBcCGEP94I2BbBONUBWWMDGUMX7TnA2U8Hj9y97TyG6u1wMNhZnkVnYkUTZQxMpQxfNGeD5QxEjSFIiISo1TgIiIxKpYKPCvoAJWgjJGhjOGL9nygjGGLmTlwERE5XCyNwEVE5BAqcBGRGBUTBR6tiyd/z8zONLO/mtlqM1tpZoODzlQRM0s2s6VmNjvoLBUxs/pmNtXMPiv7WrYPOlN5ZvZg2b/xCjObYma1oyBTtpltNbMVh2w72czmmtmaso8NojDj82X/1svN7E0zqx9tGQ95bqiZuZk1CiLbkUR9gR+yePLPgNZAbzOLtlWKS4CH3L0VcClwXxRmBBgMrA46xFGMAv7s7ucCFxBlWc2sCfAAkO7ubThwGeVewaYCIAfoVG7bMOBddz8HeLfscZBy+GHGuUAbdz8f+Bx4pLpDlZPDDzNiZmdyYPH2jdUd6FiivsA5ZPFkdy8Cvl88OWq4+1fuvqTs/k4OFE+TYFMdzszOAG4AxgedpSJmdiLwH8AEAHcvcvfvgk1VoRpAHTOrAdQlClahcvf3gW/Lbb4RmFR2fxLQrVpDlVNRRnef4+4lZQ8XcGBVr8Ac4esIMBL4NRB1R3zEQoFXtHhyVJXjocysGXARsDDYJD/wIge+CUuDDnIELYACYGLZNM94M0sNOtSh3D0fGMGBkdhXwHZ3nxNsqiM61d2/ggMDDOCUgPMcyx1AbtAhyjOzrkC+u38SdJaKxEKBWwXbou5/QgAzqwdMA4a4+46g83zPzDoDW919cdBZjqIGcDEwxt0vAnYT/K/9hymbR74RaA6cDqSa2S+DTRX7zOwxDkxDTg46y6HMrC7wGPB40FmOJBYKPCYWTzazFA6U92R3nx50nnI6AF3N7J8cmIK6ysxeDTbSD2wGNrv797+5TOVAoUeTa4D17l7g7sXAdOCygDMdyRYzOw2g7OPWgPNUyMz6AZ2BWzz6Tko5iwP/WX9S9rNzBrDEzBoHmuoQsVDgUb94spkZB+ZuV7v774POU567P+LuZ7h7Mw58/ea5e1SNHN39a2CTmbUs23Q1sCrASBXZCFxqZnXL/s2vJsreaD3ELKBf2f1+wMwAs1TIzDoBDwNd3X1P0HnKc/dP3f0Ud29W9rOzGbi47Hs1KkR9gZe9yfH94smrgTeicPHkDkBfDoxsl5Xdrg86VAwaBEw2s+XAhcAzAec5TNlvB1OBJcCnHPj5CfxUazObAnwEtDSzzWZ2J/AscK2ZreHAERTPRmHGDOAEYG7Zz8zLUZgxqulUehGRGBX1I3AREamYClxEJEapwEVEYpQKXEQkRqnARURilApcRCRGqcBFRGLU/wHy6eAEY7qzywAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(time, g_k*(ode_solution**4), 'k')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/mi/thempel/Dokumente/Papers/2019-cecam-msm/notebooks/local.py:59: RuntimeWarning: invalid value encountered in true_divide\n",
      "  alpha = (0.01*(10.-Vm))/(np.exp((10.-Vm)/10.)-1)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x7f77007c2e48>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAeE0lEQVR4nO3deXiU9b3+8fcnCQnIvgQMhLAVEFRQiAhWQVsX0FbqVlFr3Spqq+e0PdXaep3upz+ttVatiLivdfcnWtwV3CuL7AgkbAlrAFlDQpbP+SOjZ4gJmcDMPLPcr+vKlZlnnszcfCe5+ebJs5i7IyIiyS8j6AAiIhIdKnQRkRShQhcRSREqdBGRFKFCFxFJEVlBvXCXLl28d+/eQb28iCSypUvrPg8cGGyOBDR79uzN7p7b0GOBFXrv3r2ZNWtWUC8vIonsxBPrPk+fHmSKhGRmqxt7TJtcRERShApdRCRFqNBFRFKECl1EJEWo0EVEUoQKXUQkRajQRURShApdRCROdldWc++MYmau2hqT5w/swCIRkXSxs6KKRz9ezf3vr+CL8iquObEfx/TuFPXXUaGLiMRIRVUNj3y0intmFLOtvIqTBuZy3bf7M6ygY0xeT4UuIhJl1TW1PDu7lDveWs6GHRWMGZDLz08ZwNCeHWL6uip0EZEoqa11pi1cz21vLGPl5t0MK+jA3yccxci+nePy+ip0EZEo+GTFFv70r8UsXLuDgd3act8PCzl5UFfMLG4ZVOgiIgehZGs5N7/6Of9asJ4eHVpx+/lDOXNoDzIz4lfkX1Khi4gcgPK91UyeXsy9763ADH528gAmju5Lq+zMwDKp0EVEmunNxRv57UsLWbe9gjOHdufGcYfRvUOroGOp0EVEIrVu2x5+N3URbyzeyIBubXhmwihG9In+/uQHKqJCN7OxwB1AJnC/u99c7/H2wONAQeg5/+ruD0U5q4hIIGpqnUc+WsVtbyylxp0bxg7kR8f3JTsrsQ62b7LQzSwTuBs4BSgFZprZVHdfHLbaT4DF7v5dM8sFlprZE+6+NyapRUTipLhsF794dh6frdnGmAG5/HH8ERR0PiToWA2KZIY+Aihy9xUAZvYUMB4IL3QH2lrd/jltgK1AdZSziojETU2t89CHK7n19aW0bJHJ388/ivFHdY/rbojNFUmh9wBKwu6XAsfWW+cfwFRgHdAWON/da+s/kZlNBCYCFBQUHEheEZGYW7V5N9c/N4+Zq77g5EFd+fNZR9K1XcugYzUpkkJv6L8jr3f/NGAu8C2gH/Cmmb3v7jv2+SL3KcAUgMLCwvrPISISKHfn2dml/PalRWRlGredN5Szh/VI6Fl5uEgKvRToGXY/n7qZeLjLgJvd3YEiM1sJHAZ8GpWUIiIxtqOiipteXMjL89Yxsm8nbj//KPLaB78rYnNEUugzgf5m1gdYC0wALqy3zhrg28D7ZtYNGAisiGZQEZFYmb36C/7zqc9Yv72C608byNVj+gVypOfBarLQ3b3azK4FXqdut8UH3X2RmV0denwy8EfgYTNbQN0mml+6++YY5hYROWjuzpT3VvCX15eS174lz1w1iuG9YnNq23iIaD90d58GTKu3bHLY7XXAqdGNJiISOzsrqrj+2fm8tmgDpx95KDefM4R2LVsEHeug6EhREUk7RZt2ctVjs1m1pZybTh/Ej07okzR/+NwfFbqIpJVpC9Zz/bPzaJWdyeNXHMuofvE5V3k8qNBFJC24O39/azl3vL2cows6MOmiYUm3F0tTVOgikvIqqmq44bn5TJ23jnOG5fPns48gJyu409zGigpdRFLa5l2VTHx0FnPWbOOGsQO5Zky/lNhe3hAVuoikrGUbd3L5wzPZvKuSey4axrgj84KOFFMqdBFJSbNWbeXyh2eS0yKTpyeOYmjPDkFHijkVuoiknHc+38iPn5hDXvtWPHr5CHp2SszT3UabCl1EUsrzs0u54fn5DM5rx0OXHUOXNjlBR4obFbqIpIz731/Bn/61hG9+ozP3XlxIm5z0qrj0+teKSMq6463l3P7WMk4/8lBuP/+olNwtsSkqdBFJau7O7W8t5863l3Pu8HxuOWdIUp4pMRpU6CKStNydv725jLveKeL7hfncfPYQMtK0zEGFLiJJyt356xtLufvdYiYc05M/n3VkWpc5qNBFJAm5O7e+vpRJ04u5YEQB//O9I9K+zAEygg4gItJck6YXM2l6MRceqzIPp0IXkaTy2CerufX1pXzvqO78abzKPJwKXUSSxktz1/KblxZy8qCu3HreUJV5PSp0EUkKby/ZyM+fmceI3p34x4XDaJGp+qpPIyIiCe/TlVv58RNzOLx7O+6/pJCWLdLvoKFIqNBFJKEVl+3iykdn0aNjKx6+bARtk/xCzrGkQheRhLVlVyWXPTSTrAzj4UtH0Kl1dtCREpr2QxeRhFTrzpWPzmLjjgr+OXEkBZ3T4xS4B0OFLiIJx4Gist18VrKNSRcOY1hBx6AjJQVtchGRhFOytZytuyr51bjDUv6ycdGkGbqIJJSX5q6l27Y9dG3XkitP6Bt0nKSiGbqIJIyFa7dzw3PzaduyBb07t8ZMBw41hwpdRBLC1t17ueqx2XRqnc2Abm3QQaDNp0IXkcBV19Ry7ZNzKNtVyb0XD9dRoAdIoyYigbv51c/5qHgLfz7rSIbkdwg6TtJSoYtIoF6au5b7P1jJpcf15tzh+UHHSWoqdBEJzIqyXfzqhQUc07sjN50xKOg4SU+FLiKBqKiq4donPyMnK4M7Lzha282jQPuhi0gg/jxtCYvX7+DBSwvJa98q6DgpQf8likjcvbZwPY9+vJofHd+Hbx3WLeg4KUOFLiJxVbK1nOufm8/Q/PbcMPawoOOkFBW6iMRNVU0t1/3zM3C464JhZGepgqJJ29BFJG7ueqeIuSXb+MeFR+t0uDEQ0X+PZjbWzJaaWZGZ3djIOiea2VwzW2RmM6IbU0SS3bySbdz9bhFnH92D7wzpHnSclNTkDN3MMoG7gVOAUmCmmU1198Vh63QAJgFj3X2NmXWNVWARST4VVTX87Jm5dG2bw2/PPDzoOCkrkhn6CKDI3Ve4+17gKWB8vXUuBF5w9zUA7r4pujFFJJnd8trnrCjbza3nDqV9K10TNFYiKfQeQEnY/dLQsnADgI5mNt3MZpvZDxt6IjObaGazzGxWWVnZgSUWkaTyUdFmHvpwFZce15vj+3cJOk5Ki6TQGzqJpde7nwUMB84ATgP+28wGfO2L3Ke4e6G7F+bm5jY7rIgklx0VVfzi2Xn07dKaX2oXxZiLZC+XUqBn2P18YF0D62x2993AbjN7DxgKLItKShFJSn98eTEbdlTw/DXH0So7M+g4KS+SGfpMoL+Z9TGzbGACMLXeOi8BJ5hZlpkdAhwLLIluVBFJJu8vL+PZ2aVcPaYfR+siz3HR5Azd3avN7FrgdSATeNDdF5nZ1aHHJ7v7EjN7DZgP1AL3u/vCWAYXkcS1Z28Nv35xAX27tOY/vt0/6DhpI6IDi9x9GjCt3rLJ9e7fCtwavWgikqxuf2sZJVv38PTEkbRsoU0t8aLjbkUkqhau3c7976/gghEFHNu3c9Bx0ooKXUSiprqmll8+P58ubXK4cZz2aok3nctFRKLmgQ9WsmjdDib/YJgOIAqAZugiEhWrt+zmb28u49TB3Rh7RF7QcdKSCl1EDpq787upi2iRmcEfxh8RdJy0pUIXkYP29pJNvLu0jJ+e3J9D27cMOk7aUqGLyEGpqKrh968son/XNlxyXO+g46Q1/VFURA7K5BnFlGzdw5NXHkuLTM0Rg6TRF5EDVrK1nHumF/OdIXkc109nUgyaCl1EDtgfXllMZoZx0xmDgo4iqNBF5AC9u3QTby7eyHXf6k9e+1ZBxxFU6CJyAPZW1/KHlxfTt0trrji+T9BxJESFLiLN9tgnq1m5eTf//Z3BZGepRhKF3gkRaZZt5Xu58+3lnNC/CycO1JXHEokKXUSa5a53ithRUcWvTx+EWUNXqJSgqNBFJGKrNu/m0Y9X8f3hPRmU1y7oOFKPCl1EInbLa5/TIjOD/zr1a9eAlwSgQheRiMxctZVXF27gqtH96NpO52tJRCp0EWlSba3zp38toVu7HK4crd0UE5UKXUSa9MqC9cwr2cZ/nTqQQ7J1CqhEpUIXkf2qqqnltjeWctihbTlnWH7QcWQ/VOgisl9Pzyxh9ZZybhg7kMwM7aaYyFToItKoPXtruPPt5RT26shJA7sGHUeaoEIXkUY98vEqNu2s5PrTBuogoiSgQheRBu2oqOKe6cWMGZDLsX07Bx1HIqBCF5EG3ffeCrbvqeL60wYGHUUipEIXka8p21nJAx+s5IwheRzRo33QcSRCKnQR+Zq73y2isrqWn5+iQ/yTiQpdRPaxdtsenvz3Gs4dlk+/3DZBx5FmUKGLyD4mvVuE4/zHyf2DjiLNpEIXka+s3baHZ2aV8P3CnvTooOuEJhsVuoh8ZdK7RQD8+KRvBJxEDoQKXUQAzc5TgQpdRADNzlOBCl1Evpqdn3+MZufJTIUuIv83Oz9Rs/NkpkIXSXPhs/Pump0ntYgK3czGmtlSMysysxv3s94xZlZjZudGL6KIxNI90zU7TxVNFrqZZQJ3A+OAwcAFZja4kfVuAV6PdkgRiY1NOyt4ZlYp5w7X7DwVRDJDHwEUufsKd98LPAWMb2C964DngU1RzCciMfTQh6uorqnlqtF9g44iURBJofcASsLul4aWfcXMegBnAZP390RmNtHMZpnZrLKysuZmFZEo2lFRxeMfr2bckXn07tI66DgSBZEUekOXKfF69/8O/NLda/b3RO4+xd0L3b0wNzc30owiEgNP/nsNOyuruWZMv6CjSJRkRbBOKdAz7H4+sK7eOoXAU6FLVHUBTjezanf//1FJKSJRVVFVwwMfrOSE/l10vvMUEkmhzwT6m1kfYC0wAbgwfAV37/PlbTN7GHhFZS6SuF78bC1lOyu54/yjgo4iUdRkobt7tZldS93eK5nAg+6+yMyuDj2+3+3mIpJYamqde2cUMyS/PaP66VqhqSSSGTruPg2YVm9Zg0Xu7pcefCwRiZXXFm5g1ZZy7rloGKHNpJIidKSoSBpxdybPKKZPl9acevihQceRKFOhi6SRD4u2sGDtdq4a3ZfMDM3OU40KXSSNTJ5RTNe2OZw1rEfTK0vSUaGLpIkFpdv5oGgzVxzfh5yszKDjSAyo0EXSxOQZxbRtmcWFxxYEHUViRIUukgZWbt7NtIXruXhkL9q2bBF0HIkRFbpIGpjyXjEtMjO47Jt9ml5ZkpYKXSTFbdpRwfOz13Le8Hxy2+YEHUdiSIUukuIe+HAl1bW1TNQpclOeCl0khe2oqOLJT9Zw+pF59OqsU+SmOhW6SAp7/JPV7Kys5mqdIjctqNBFUlRFVQ0PfrBKp8hNIyp0kRT1/JxSNu+q5JoTNTtPFyp0kRRUU+tMeW8FQ/PbM6qvTpGbLlToIino1YXrWb2lnGtO7KdT5KYRFbpIinF37pleTN8urTllsE6Rm05U6CIp5oOizSxat4OrxugUuelGhS6SYu6ZXky3djl872idIjfdqNBFUsi8km18VLxFp8hNUyp0kRQyeUYx7VpmccEInSI3HanQRVJEcdkuXlu0gYtH6RS56UqFLpIi7ntvBdmZGVx6nE6Rm65U6CIpYOOOCl6Ys5bzCnWK3HSmQhdJAQ9+EDpF7gk6zD+dqdBFktz28ioe/2Q1ZwzpTkHnQ4KOIwFSoYskucf/vZrde2u4eowuYJHuVOgiSayiqoaHPlzJ6AG5HN5dp8hNdyp0kST27OxSNu/ayzW6gIWgQhdJWtU1tdz33gqG9uzAyL6dgo4jCUCFLpKkpi3cwJqt5VwzRqfIlToqdJEk5O5Mnl5M39zWnDq4W9BxJEGo0EWS0IxlZSxev4OrR/cjQ6fIlRAVukgSmvRuMXntW+oUubIPFbpIkvl05VY+XbWViaP7kp2lH2H5P/puEEkyk6YX0al1NhOO0SlyZV8qdJEksnDtdqYvLeOK4/vQKlsXsJB9qdBFksg904tpm5PFD0b2CjqKJKCICt3MxprZUjMrMrMbG3j8IjObH/r4yMyGRj+qSHorLtvFtIXruXhUL9q30gUs5OuaLHQzywTuBsYBg4ELzGxwvdVWAmPcfQjwR2BKtIOKpLvJ04vJzszg8uN1AQtpWCQz9BFAkbuvcPe9wFPA+PAV3P0jd/8idPcTID+6MUXS29pte3jxs7VcMKKALm10AQtpWCSF3gMoCbtfGlrWmCuAVw8mlIjs6773VgBw5WidIlcalxXBOg0dhuYNrmh2EnWFfnwjj08EJgIUFGiXK5FIbNhewZOfruHsYT3o0aFV0HEkgUUyQy8FeobdzwfW1V/JzIYA9wPj3X1LQ0/k7lPcvdDdC3Nzcw8kr0jamTS9iNpa57pv9Q86iiS4SAp9JtDfzPqYWTYwAZgavoKZFQAvABe7+7LoxxRJT+u27eGpT0s4rzCfnp10eTnZvyY3ubh7tZldC7wOZAIPuvsiM7s69Phk4DdAZ2BS6DSe1e5eGLvYIulh0vQiHOcnJ30j6CiSBCLZho67TwOm1Vs2Oez2j4AfRTeaSHpbu20PT88s4fuFPcnvqNm5NE1HiookqH+8U4Rhmp1LxFToIgmoZGs5z84q4fxjetJde7ZIhFToIgno7neLyDDjxyfp4s8SORW6SIJZs6Wc52aXcuGxBeS11+xcIqdCF0kwf31jKVmZxjUnanYuzaNCF0kgC0q3M3XeOq44vg/d2rUMOo4kGRW6SIJwd25+bQkdD2nBVWM0O5fmU6GLJIj3lm/mw6ItXPet/rRrqfOdS/Op0EUSQE2tc/Orn9OzUysuGqkT18mBUaGLJIDn55SyZP0OfnHqQHKydK1QOTAqdJGA7aio4i+vLeXogg58d0j3oONIEovoXC4iEjt3vb2cLbsrefDSQjIyGrr8gEhkNEMXCVDRpl089OEqzi/syZD8DkHHkSSnQhcJiLvz+5cX0So7k1+cNjDoOJICVOgiAXl90UbeX76Zn508QBd+lqhQoYsEYPueKn7z0kIG5bXj4lG9go4jKUKFLhKA/zdtCZt3VfKXc4bQIlM/hhId+k4SibOPijbz1MwSrhzdlyPz2wcdR1KICl0kjvbsreHGFxbQu/Mh/OzkAUHHkRSj/dBF4ujmV5ewZms5/7xyJC1b6IhQiS7N0EXi5M3FG3nk49VccXwfRvXrHHQcSUEqdJE42LC9guufm8fh3dtxw1jtcy6xoUIXibGaWuenT3/G3upa7rrgaJ18S2JG29BFYuyOt5bxyYqt3HruEPrmtgk6jqQwzdBFYujleeu4850izh2ez7nD84OOIylOhS4SI/NKtvGLZ+dxTO+O/M9ZR2CmMylKbKnQRWJg/fY9XPnoLHLb5jD5B8O13VziQoUuEgO3vbGM3ZXVPHDJMXTWibckTvRHUZEY+MP4w/nByF4MPLRt0FEkjWiGLhIl28r3sqB0OwCHZGdxVE9dsELiS4UuEgWlX5Rz7uSPufyRmezZWxN0HElT2uQicpBmrdrKj5+Yw56qGqZcXEirbP0BVIKhQhc5QO7Ofe+v4JbXlpLfsRWPXXGstplLoFToIgdg084Kfv3CQt5aspFxRxzKLecOoV3LFkHHkjSnQhdphuqaWh77ZDV/e2MZFdU1/OY7g7nsm7110JAkBBW6SARqa503l2zk9jeX8fmGnZzQvwu/P/NwnZtFEooKXWQ/9uyt4eX567h3RjHFZbvp2akVky4axrgjDtWsXBKOCl2knoqqGj4q3szUuet4c/FGdu+tYVBeO+6YcBRnHJlHli7qLAkqokI3s7HAHUAmcL+731zvcQs9fjpQDlzq7nOinFUk6tydDTsq+Hz9TuaVbuPj4i18VrKNvdW1tG/VgjOP6s6ZQ3swsm8nzcgl4TVZ6GaWCdwNnAKUAjPNbKq7Lw5bbRzQP/RxLHBP6LNI3Lg7ldW17K2ppbKqlsrqGiqqatm+Zy9f7K7ii/K9bCuvomxXJWu/2EPpF+Ws2lLO9j1VAJjB4d3bccmoXhz3jS58s18XsrM0G5fkEckMfQRQ5O4rAMzsKWA8EF7o44FH3d2BT8ysg5nlufv6aAeesayMP76y+GvL61766xpe2vgDja3f3OdvZHW8ka9odP1G/wGNrR9QzmY+f2Nf0fznd6prQ0VeXdvYi+0jOyuDHh1akd+xFWcMyWPQoW0ZeGg7Dstrq10PJalFUug9gJKw+6V8ffbd0Do9gH0K3cwmAhMBCgoKmpsVgDY5WQzs1sjBG438RtzYL8qN/Qrd+Pqxff7G8zfyPM3OE6Xnb+Y/INY5szKNnKxMcrIyyGmRQU5WJtlZGXX3szJo36oFHQ/JpuMh2XRo3YK2OVnafCIpKZJCb+g7v/6EKZJ1cPcpwBSAwsLCZs496wzv1ZHhvToeyJeKiKS0SDYQlgI9w+7nA+sOYB0REYmhSAp9JtDfzPqYWTYwAZhab52pwA+tzkhgeyy2n4uISOOa3OTi7tVmdi3wOnW7LT7o7ovM7OrQ45OBadTtslhE3W6Ll8UusoiINCSi/dDdfRp1pR2+bHLYbQd+Et1oIiLSHNrJVkQkRajQRURShApdRCRFqNBFRFKENXaoeMxf2KwMWH2AX94F2BzFONGUqNmUq3kSNRckbjblap4DzdXL3XMbeiCwQj8YZjbL3QuDztGQRM2mXM2TqLkgcbMpV/PEIpc2uYiIpAgVuohIikjWQp8SdID9SNRsytU8iZoLEjebcjVP1HMl5TZ0ERH5umSdoYuISD0qdBGRFJHwhW5m55nZIjOrNbPCeo/9ysyKzGypmZ0Wtny4mS0IPXanxfjyNGb2tJnNDX2sMrO5oeW9zWxP2GOTm3quKOf6nZmtDXv908Mea3Ds4pjtVjP73Mzmm9mLZtYhtDzQMQtlGBsalyIzuzHerx+Wo6eZvWtmS0I/A/8ZWt7o+xrHbKtCP2NzzWxWaFknM3vTzJaHPsf9SjRmNjBsXOaa2Q4z+2kQY2ZmD5rZJjNbGLas0TGKys+kuyf0BzAIGAhMBwrDlg8G5gE5QB+gGMgMPfYpMIq6Kym9CoyLY97bgN+EbvcGFgY4dr8DftHA8kbHLo7ZTgWyQrdvAW5JkDHLDI1HXyA7NE6DA8qSBwwL3W4LLAu9dw2+r3HOtgroUm/ZX4AbQ7dv/PI9Dfi93AD0CmLMgNHAsPDv58bGKFo/kwk/Q3f3Je6+tIGHxgNPuXulu6+k7lzsI8wsD2jn7h973Ug9CnwvHllDvwl8H/hnPF7vIDQ4dvEM4O5vuHt16O4n1F3lKhF8dVF0d98LfHlR9Lhz9/XuPid0eyewhLpr9Saq8cAjoduPEKefu/34NlDs7gd6RPpBcff3gK31Fjc2RlH5mUz4Qt+Pxi5M3SN0u/7yeDgB2Ojuy8OW9TGzz8xshpmdEKcc4a4NbdZ4MOzXu8bGLiiXU/eb1JeCHLNEGxugblMUcDTw79Ciht7XeHLgDTObbXUXfwfo5qErlYU+dw0gV7gJ7Du5CnrMoPExisr3XUIUupm9ZWYLG/jY38yosQtTR3TB6hhlvIB9v4HWAwXufjTwc+BJM2t3sFmakeseoB9wVCjLbV9+WQNPFfX9VyMZMzO7CagGnggtivmYNRW7gWWB7ttrZm2A54GfuvsOGn9f4+mb7j4MGAf8xMxGB5ChUVZ3ucwzgWdDixJhzPYnKt93EV2xKNbc/eQD+LLGLkxdyr6/vkflgtVNZTSzLOBsYHjY11QClaHbs82sGBgAzDrYPJHmCst3H/BK6G5cLuodwZhdAnwH+HZo81hcxqwJCXXBczNrQV2ZP+HuLwC4+8awx8Pf17hx93Whz5vM7EXqNg9sNLM8d18f2vS5Kd65wowD5nw5VokwZiGNjVFUvu8SYoZ+gKYCE8wsx8z6AP2BT0O/xuw0s5Ghbdo/BF6KQ56Tgc/d/avNPWaWa2aZodt9QxlXxCHLl6+fF3b3LODLv7Y3OHbxyhXKNhb4JXCmu5eHLQ90zIjsouhxEfr+fQBY4u5/C1ve2Psar1ytzaztl7ep+wP3QurG6ZLQapcQn5+7xuzz23LQYxamsTGKzs9kkH+FjvAvxWdR979XJbAReD3ssZuo+2vwUsL2ZAEKqXvDioF/EDoiNsY5HwaurrfsHGARdX+9ngN8N85j9xiwAJgf+obJa2rs4pitiLpthnNDH5MTYcxCGU6nbo+SYuCmeL9+WI7jqfu1e37YOJ2+v/c1Trn6ht6feaH36qbQ8s7A28Dy0OdOAY3bIcAWoH3YsriPGXX/oawHqkIddsX+xigaP5M69F9EJEUk8yYXEREJo0IXEUkRKnQRkRShQhcRSREqdBGRFKFCFxFJESp0EZEU8b+aRlWaqBXhmwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the steady state open occupancy\n",
    "vs = np.arange(-100, 100)\n",
    "plt.plot(vs, kchannel_params(vs)[2]**4)\n",
    "plt.axvline(V_clamp, color='r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## model generation & convergence analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct a microstate trajectory\n",
    "msmlag = 0.0001 # in seconds\n",
    "tot_T = 100000 # Steps\n",
    "\n",
    "alpha, beta, ss, tau = kchannel_params(V_clamp)\n",
    "c = Channel(alpha*1000, beta*1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Macrostate Transition Matrix:\n",
      "[[0.8086 0.1764 0.0144 0.0005 0.    ]\n",
      " [0.0047 0.8488 0.1388 0.0076 0.0001]\n",
      " [0.     0.0099 0.8904 0.097  0.0026]\n",
      " [0.     0.0001 0.0155 0.9335 0.0509]\n",
      " [0.     0.     0.0002 0.0217 0.9781]] \n",
      "\n",
      "single sub-system transition matrix:\n",
      "[[0.9483 0.0517]\n",
      " [0.0055 0.9945]]\n"
     ]
    }
   ],
   "source": [
    "Tm = c.transition_matrix(msmlag)\n",
    "Tm_singlet = c.transition_matrix_singlet(msmlag)\n",
    "\n",
    "with np.printoptions(precision=4, suppress=True):\n",
    "    print(\"Macrostate Transition Matrix:\")\n",
    "    print(Tm, '\\n')\n",
    "    print(\"single sub-system transition matrix:\")\n",
    "    print(Tm_singlet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate a trajectory in empirical state space \n",
    "traj = msmtools.generation.generate_traj(Tm, tot_T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# is it connected?\n",
    "# in rare cases, tot_T not large enough to sample all 5 states reversibly\n",
    "msmtools.estimation.is_connected(msmtools.estimation.count_matrix(traj, 1, sparse_return=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate from sub-system transition matrix (single gate)\n",
    "singlet_traj = msmtools.generation.generate_traj(Tm_singlet, tot_T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Empirical P_open: 0.67777\n",
      "SS from eigendecomposition: 0.6664516597482566\n",
      "from independent sub-systems:  0.6664516597482637\n"
     ]
    }
   ],
   "source": [
    "# Validate that the steady state open probability is as we expect\n",
    "P_open = np.count_nonzero(traj == 4) / tot_T\n",
    "print('Empirical P_open:', P_open)\n",
    "print('SS from eigendecomposition:', msmtools.analysis.stationary_distribution(Tm)[4])\n",
    "print('from independent sub-systems: ', \n",
    "      msmtools.analysis.stationary_distribution(Tm_singlet)[1]**4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## dMSM vs. MSM comparison\n",
    "### estimate sampling and convergence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "ntrials = 10  # repeat for different realizations to assess error\n",
    "npoints = 40 # number of points on sampling axis\n",
    "min_sampling_log = 2 # start conducting relaxation experiment here\n",
    "\n",
    "p0 = [1, 0, 0, 0, 0] # prepare in closed state\n",
    "observable = [0, 0, 0, 0, 1] # K conductance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b26dbabf4d364865af6703fc90797d43",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=10), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/storage/mi/thempel/anaconda3/envs/py37_oct19/lib/python3.7/site-packages/msmtools/analysis/dense/decomposition.py:540: ImaginaryEigenValueWarning: Using eigenvalues with non-zero imaginary part\n",
      "  warnings.warn('Using eigenvalues with non-zero imaginary part', ImaginaryEigenValueWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "sampling_arr_full = np.zeros((Tm.shape[0]-1+2, ntrials, npoints)) + np.NaN  # classical MSM\n",
    "sampling_arr_dec = sampling_arr_full.copy() # dMSM\n",
    "\n",
    "relaxation_arr_full = np.zeros((150, ntrials, npoints)) + np.NaN\n",
    "relaxation_arr_dec = relaxation_arr_full.copy()\n",
    "\n",
    "# conduct `ntrials` times the same experiment\n",
    "for trial in tqdm(range(ntrials)):\n",
    "    traj = msmtools.generation.generate_traj(Tm, tot_T)\n",
    "\n",
    "    # when using dMSMs, we'd have a full channel in the simulation box, i.e. four \n",
    "    # gates that we extract trajectories of.\n",
    "    # cannot disassemble above `traj` as information about single channels got\n",
    "    # lost during empirical state coarse graining\n",
    "    \n",
    "    singlet_trajs = msmtools.generation.generate_trajs(Tm_singlet, 4, tot_T) \n",
    "    \n",
    "    for n_end, end in enumerate(np.logspace(min_sampling_log, len(str(tot_T)) -1, npoints)):\n",
    "        \n",
    "        # in empirical state space\n",
    "        ccmat = msmtools.estimation.connected_cmatrix(\n",
    "            msmtools.estimation.count_matrix(traj[:int(end)], 1, sparse_return=False))\n",
    "\n",
    "        sampling_arr_full[0, trial, n_end] = msmlag * 1000 * end  # unit: ms\n",
    "        sampling_arr_full[1, trial, n_end] = ccmat.shape[0] / 5\n",
    "        \n",
    "        if ccmat.shape[0] > 1:\n",
    "            tmat = msmtools.estimation.tmatrix(ccmat)\n",
    "\n",
    "            sampling_arr_full[2:2+tmat.shape[0]-1, trial, n_end] = msmtools.analysis.timescales(tmat)[1:]\n",
    "\n",
    "            if tmat.shape[0] == 5:\n",
    "                msm = pyemma.msm.markov_model(tmat, dt_model='0.1 ms')\n",
    "                times, rel = msm.relaxation(p0, observable, maxtime=15)\n",
    "                relaxation_arr_full[:, trial, n_end] = rel\n",
    "\n",
    "        ccmat, tmat = None, None\n",
    "        singlet_ccmats = [msmtools.estimation.connected_cmatrix(\n",
    "            msmtools.estimation.count_matrix(singlet_traj[:int(end)], 1, sparse_return=False))\n",
    "                           for singlet_traj in singlet_trajs]\n",
    "\n",
    "        sampling_arr_dec[0, trial, n_end] = msmlag * 1000 * end  # unit: ms\n",
    "        sampling_arr_dec[1, trial, n_end] = np.mean([ccmat.shape[0] / 2 for ccmat in singlet_ccmats])\n",
    "\n",
    "        if all([ccmat.shape[0] > 1 for ccmat in singlet_ccmats]): # exclude rare-sampling case that only visited a single state (triggers estimation error)\n",
    "            singlet_tmats = [msmtools.estimation.tmatrix(ccmat) for ccmat in singlet_ccmats] # subsystem tmatrices\n",
    "            tmat_full_uncg = np.kron(np.kron(np.kron(singlet_tmats[0], singlet_tmats[1]), singlet_tmats[2]), singlet_tmats[3])\n",
    "            tmat_full = cg_transition_matrix(tmat_full_uncg, chi)  \n",
    "            \n",
    "    \n",
    "            sampling_arr_dec[2:2+tmat_full.shape[0]-1, trial, n_end] = msmtools.analysis.timescales(tmat_full)[1:]\n",
    "            \n",
    "            # tmat_full is fully connected\n",
    "            msm = pyemma.msm.markov_model(tmat_full, dt_model='0.1 ms')\n",
    "            times, rel = msm.relaxation(p0, observable, maxtime=15)\n",
    "            relaxation_arr_dec[:, trial, n_end] = rel            \n",
    "\n",
    "# occasional errors due to finite sampling can be ignored -> NaNs in output arrary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "normalization factor: 24\n"
     ]
    }
   ],
   "source": [
    "# since our observable is not normed to be a conductance, we have to\n",
    "# multiply by a constant factor\n",
    "# derive from relaxation that uses most data\n",
    "normalization_factor =  g_k*(ode_solution**4)[-1] / relaxation_arr_dec[-1, :, -1].mean()\n",
    "print(f\"normalization factor: {np.round(normalization_factor, 0)[0].astype(int)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### plot convergence and relaxation experiment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.matplotlib.rcParams['axes.titlesize'] = 12\n",
    "plt.matplotlib.rcParams['axes.labelsize'] = 20\n",
    "plt.matplotlib.rcParams['legend.fontsize'] = 11\n",
    "plt.matplotlib.rcParams['xtick.labelsize'] = 11\n",
    "plt.matplotlib.rcParams['ytick.labelsize'] = 11\n",
    "plt.matplotlib.rcParams['xtick.minor.pad'] = 2\n",
    "plt.matplotlib.rcParams['xtick.major.pad'] = 3\n",
    "plt.matplotlib.rcParams['ytick.minor.pad'] = 2\n",
    "plt.matplotlib.rcParams['ytick.major.pad'] = 3\n",
    "plt.matplotlib.rcParams['axes.labelpad'] = 1\n",
    "plt.matplotlib.rcParams['axes.labelsize'] = 11\n",
    "from matplotlib import rc\n",
    "rc('text', usetex=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_segmented(x, y, cmap, ax=None, norm=None, lw=None, show_cbar=False, **kwargs):\n",
    "    if ax is None:\n",
    "        fig, ax =  plt.subplots()\n",
    "    if norm is None:\n",
    "        norm = plt.matplotlib.colors.LogNorm(x.min(), x.max())\n",
    "    points = np.array([x, y]).T.reshape(-1, 1, 2)\n",
    "    segments = np.concatenate([points[:-1], points[1:]], axis=1)\n",
    "    #norm = plt.Normalize(x.min(), x.max())\n",
    "    \n",
    "    \n",
    "    lc = plt.matplotlib.collections.LineCollection(segments, cmap=cmap, norm=norm, **kwargs)\n",
    "    lc.set_array(x)\n",
    "    if lw is not None:\n",
    "        lc.set_linewidth(lw)\n",
    "\n",
    "    line = ax.add_collection(lc)\n",
    "    if show_cbar:\n",
    "        ax.get_figure().colorbar(line, ax=ax)\n",
    "    return ax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "lw = 2.5\n",
    "ms = 8\n",
    "alpha = .8\n",
    "vspanwidthpercentage = 0.06\n",
    "reduce_max_conductance = 1.  # percentage of max_conductance to scale ylim with\n",
    "dashes = (5, 5) #(length, spacing)\n",
    "example_relaxation_curve_indices = [13, 26, 39]\n",
    "\n",
    "its_ref = msmtools.analysis.timescales(Tm)\n",
    "sampling_time_ms = sampling_arr_full[0, 0, :]\n",
    "n_lengths = len(sampling_time_ms)\n",
    "\n",
    "cmap = plt.matplotlib.cm.viridis_r\n",
    "norm = plt.matplotlib.colors.LogNorm(sampling_time_ms.min(), sampling_time_ms.max())\n",
    "\n",
    "cols = cmap(norm(sampling_time_ms))\n",
    "plot_colorscaled = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASsAAAFSCAYAAACqifm9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydeXhURda439NJCGQDkrAYVpUdUTZ3EHSQ0VHUT/0EFHH3wxlG/cZRHHXc9XMU992fG8q4oiIqioICgrIFGAUEFCGgILJD9u38/qjb4abTWXpLupN6n6eevrdu3arTXdXnVtWtOkdUFYvFYol2PA0tgMVisdQFq6wsFktMYJWVxWKJCayyslgsMYFVVhaLJSawyspiscQEVllFMSLyqojc29ByNEVE5FIRWRDB/D8VkUtCzCOiMkYbVlnFACIyXERURN73iT/KiZ/rijtbRFaKyH4R2Skic0Skq3PtTif9tT75XO/E3xn5b2MBUNXTVXVKpPIXka5OnS73ic8UkWIR2eSKGyIi34jIPhHZLSILReRo59qlTj6P+ORzjhP/aqS+gy9WWcUOO4ATRCTDFXcJsN57IiLdgNeAG4CWwKHAM0C56571zn1uxrvzsTQqkkXkCNf5hcBG74mIpAEfA08C6UAH4C6gyHXPBmC0iMS74uq9zVhlFUWIyAARWS4iB0TkbaC563IxMB0Y46SNAy4A/u1K0x/YqKpz1HBAVd9T1c2uNEuBJBHp6+TTF2jhxDc5RKSTiLwvIjtEZJeIPFVNusdFZIvTY80WkaGua8eIyDLn2nZvL0REmovIVCffvSKyVETaOdfmisiVrjyuEpEfnLpfIyIDnfibRWSDK/6/AvyKr1P54TQe80Dz0gNAVd9U1TJVLVDVz1X1O1ea34DvgT86MqUDJwAzApQlJKyyihJEpBlGGb2OecK9C5znk+w1TGMD03BWA1td15cDvUTkURE5WURSqinudVc+l1C58TYZHIX/MZADdMX0Kt6qJvlSzMMgHXgDeFdEvA+Tx4HHVTUNOBx4x4m/BNPD7QRkABOAAj9y/DdwJ6ZO0oCzgF3O5Q3AUCefu4CpInJIAF9zKjBGROJEpDeQCix2XV8PlInIFBE5XURaV5OPu+2NAT6kcu8r4lhlFT0cByQAj6lqiapOw6e3o6rfAOki0pOqT0hU9WdgOOZP9w6w05mk91VaU4GxIpKAaXhTI/B9YoFjgCzgRlXNU9VCVfU7Ya2qU1V1l6qWqurDQCLQ07lcAnQTkUxVzVXVRa74DKCb02vJVtX9frK/EnhQVZc6PeKfVDXHKfddVd2qquWq+jbwoyN3XfkFWAeMwM+DyZFnCKDA/wN2iMgMbw/QxQfAcBFpiZ+2Vx9YZRU9ZAG/auWd5Tl+0r0OTAROxjSgSqjqIlW9QFXbYJ7IJwG3+qTZDPwE3A/8qKpbwvMVYo5OQI6qltaWUERucIZp+0RkL6ank+lcvgIznFrrDPXOdOJfB2YBb4nIVhF50HlA+JNjQzXljndemOx1yj3CVW5deQ24FBiLnweTqv6gqpeqakcn/yzgMZ80BcAnwG1ApqouDFCGkLHKKnrYBnQQEXHFdfaT7nXgz8BMVc2vKUNVXQq8j2mAvngn4pvkENBhC9DZZ+K4Cs781CTMHGFrVW0F7AMEQFV/VNWxQFvgX8A0EUl2esh3qWofzBzPmRwcSvnKcbifcrtgejsTgQyn3FXecgPgPeAM4Gdvj606VHUt8Co1t5nXAyw/LFhlFT18C5QC14pIvIici5/uvqpuBIbh01uCilfQV4lIW+e8F2b+Y5FvWuBtYCQH51eaIkswD4kHRCTZmRA/0U+6VEzd7ADiReR2zNwSACIyTkTaqGo5sNeJLnPmDfs5c2P7McPCMj/5vwj8XUQGiaGbo6iSMcOzHU45l+FfidSIquYBp2CGm5UQkV5Or7Gjc94J0wPz12bmAadi3hzWO1ZZRQmqWgyci+mu7wFGY3pF/tIuUNWtfi7txSin70UkF/gMM1R80E8eBao62+neN0lUtQwYBXQDNmPmd0b7SToL+BQzGZ0DFGJ6Q15OA1Y7v/njwBhVLQTaA9MwiuoHzJ/d3zDsXeA+zMT9AcyLlnRVXQM8jHmQbQf6AUENv1R1mar6G2oeAI4FFotIHkZJrcL0oHzzUOdN8+5gZAgVscb3LBZLLGB7VhaLJSawyspiscQEVllZLJaYwCori8USE9S4vqQxk5mZqV27dm1oMSy1kJ2dvdNZ4Fortk5jg0Dq1E2TVVZdu3Zl2bJlDS2GpRZEpMZFjG5sncYGgdSpGzsMtFgsMYFVVhaLJSawysoS04jI1Y4tqWU7duxoaHEsEaTJzlk1dkpKSsjLyyM3N5fi4mJKSkoqQmlpaaVzb5yqoqqUl5dXHFcXaktTHbXtmOjZsydDhgyp8/dU1ReAFwAGDx5st2M0YqyyihGKiorYtGkTGzZsYOPGjWzfvp3ff/+9IuzZs4fc3Fxyc3PJy8ujqKhe7aKFjSuvvDIgZWVpOlhlFYXk5eXx7bffsmzZMpYuXcry5cvJycmptVfSGNi4cWPtiSxNEqusooRt27bx1ltvMXPmTObPn09xcXGN6UWEjIwMWrVqRVJSEvHx8agqpaWlFBUVUVRUREFBAcXFxRQWFlJUVBQTyq5Dhw4NLYIlSrHKqgEpLy/n448/5vnnn+ezzz6jvLy8SprOnTszePBgunTpgqqyf/9+fvnlFzZs2MCmTZvYuXNnA0heO5VtCFYf5xu/ZUtTNVpqqQ2rrBoAVWXatGncddddrF69utK1rl27cvrpp/OHP/yB8vJyFixYwKeffsr77/s1bRW1+OvF1aVnV1JSEglxLI0Aq6zqmdWrV/OXv/yFefPmVcS1bNmSsWPHcumll7Jt2zaeeOIJXn755XqdJD8aaAckYfxyJTgh3vXpPt4IvBxgGTdjPCx4MHZ5vcF7/iSwYYNfU+SWYFiwAPbsgfx8E0pKoLTUBO+x+7NfPxgzJrAy/vEP2LoVVCuH8nLzedttcETAxk39U9sr6mgKwGTM/0SBI1zxm4C1wEon/LG2vAYNGqT1SXl5uT722GOamJiojvyamZmp//rXvzQ7O1svuOACTUpKqrhWXUhJSdHWrVtrSkqKejyeivgE0B6gI0CH1pKHv7CwanOrMcwOooz5teQ5FjQrK6vS7wYsq60uG6pOI05Bgerq1aqffqq6aFHg9x9+eEB1qmPGBF5Gjx415/nFF1VuCaRO3SHWelbTMWZjv/Zz7XxVXVXP8tSJvLw8LrroIj788EMAPB4Pf/nLXxgyZAgPP/wwkyZNqnJPXFwchx56KB07dmTv3r2sXbuWwsJCcnNz8bqtGQz0cUJ3DnaTvwT+EKCMNXqe8KEUo30C5TfMU8WtwcpdxweAsjJ/JsqbAFu3wjffQHY2rFljws8/mx4KwOjR8FZ1Lg2rISmp7mnj4gLL20vHjlBUBCIHg8dz8LhFi+Dy9UNMKSt1fLpVN1EbjWzbto1Ro0aRnZ0NmDmpv/3tb7zyyis8+WRlu/txcXEcffTRHHXUUaxYsYIlS5bw008/kYzxaHoGxv9Wt1rK7BSEnNdj3D/nY7xwljih1OezjOAUFRjXMLWRuHdv7YkaA3v3wqxZMHMmzJ0LmzfXnD6YFw9TpkBZGSQnQ/Pm0KwZxMdDQkLlz/h4o1iCYc6c4O4LgphSVrXwb8eN1QLgFlWt0upF5GrgajBv2SLNtm3bGDZsGD/++CMAw4cPp127dlx77bWV0sXHxzN27FiSkpKYOnUqixZVdiwyhaqumcH0SjYA6zwe1icksKa8nJ9KSvw6G6yN1bUnqRc8niayA2zUKDOn5EtcHHTrBn36mNCrF3TpAsGYvhkwIGQxo4nGoqyGquoWEUnEOGd8Chjnm0jrcWvGjh07GDFiRIWiOv3001m2bBlz586tlO6ss84iPT2dqVOnUlrq39fmZxhlVQisadOGrYcfzvLERN7asIF1W7eaJQ9RumLd2wv29+nveEAj+4NVy2mnGWWVnAynnAInnggnnACDBgU2fGtCNAplpY5HYVUtEpFngBkNKU9RURFnnXUWa9asAeCYY45h1qxZldZRtW7dmvPOO4+33nqL3Nxcv/lkZmZy5plnMmLQIFbu28fUX3/llbffZvcify7dgiMhIYGkpCRSU1Np2bIlaWlppKWlVZy3atWK1q1bk56eXnE9JSWF5OTkipCamkpycjIJCQlNp2cUKuPGwTHHwEknQWJiQ0sTE8S8shKRZCBeVfc5w8AxmDeCDYKqcs0111QM5Xr27MmSJUsqpenXrx9JSUm8+OKLVe5PTk5m7NixXHzxxRx55JG8+OKL/O/997Nt27ag5GnVqhWdO3fmsMMOo3fv3vTt25fDDjuMDh060LZtW5o3bx5UvpYQ6dLFBEudiSllJSJPYByBtgdmi8gujJPK9xyvt3HAGox79Qbh2Wef5ZVXXgGgTZs2rFu3rtL1E088kRUrVpCfX/n9W8eOHZk0aRKXXHIJe/bs4brrruPDDz8MaItMly5dOOmkk+jfvz9HHXUURx11FJmZmaF/KYslCogpZaWq1wLX+rkUFRMda9as4YYbjCPb5ORkfO0rDRgwgIULKzvUzczM5J577uGyyy5j+/btnHXWWVXmtaqjU6dOnHHGGQwfPpyhQ4eSlZUVlu9hCQLV4N+oWepETCmraKa4uJhx48ZRWFgImLVVbrp168aKFSvoA7wOjAdOuuYa7rvvPhISErjssst48803ay3n8MMP56KLLuKcc86hf//+MbWMo9GyeDFMmADTp9uhXSQJZiVpYwjhXu18++23V7tyu127dgroOaD7nZW9+R06qO7dq6+88oo2a9asxpXfLVq00Msvv1y//vprLS8vD6vc0Q7RvoL95ZdVmzUzq7UHDFAtLKx/GWKMQOrUHWzPKgysW7eOBx54ADDrhHzf+m3fvp1rgGdc9+j48QwaPpzlK6t/F9CmTRsmTpzIn//8Zzv3FG2owt13w513mvP4eLjiCrPw0hIRrLIKEVVl4sSJFfan3IoqKSmJPXv2cDPwf970qaks/9vfOP7++6u1MNCyZUtuv/12rrnmGlqEcbuCJUyowg03wKOPmvPMTJg2DYYNa1i5GjlWWYXItGnTmD17dpX4uLg48vPzuRv4pzcyI4P/GzaMW++6y29e8fHxXHfdddx66620bt06YjJbQkDVzE+98II579QJZs+GHj0aVq4mgFVWIVBcXMzNN9/s91pZWRk3clBRaVYW56Wm8kE1dqn69+/PW2+9Rc+ePSMjrCV0VOHGGw8qqu7djaKqh61bFuuKKyReeuklfv75Z7/XrgIedI7L27VjcF4eH/isuQLTm5o8eTLZ2dlWUUU7998PDz9sjrt3h/nzraKqR2zPKkjy8/O55557/F47C3jOOS5v1YrBe/awwo9N9fbt2zN//ny6d+8eOUEt4eHll40hOYAOHeCLL6B9+4aVqYlhe1ZB8uSTT/rdAtMDmIr5YcuTkjgpN9evoho5ciSbNm2yiioWWLLEzFMBZGQYRWXXU9U7VlkFwZ49eyqWKviyEfiweXNUhPOKiljox5LCHXfcwaxZs0i0G1hjgz594KyzjP2n6dOhd++GlqhJYoeBQfDQQw+xtxojceVxcXT69FOOO/VUlvgoKhHhhRde4Morr6wPMS3hIiUF3nkHVq6EgQMbWpomi1VWAfLbb7/x2GOPVXv9vvvu49SRIynxUVRxcXF8+OGHnHHGGZEW0RIJPB6rqBoYq6wC5N5776WgoMDvtREjRvDPf/6zymLPxMREvv76a44++uj6ENFiaZTYOasA+Pnnn3n++ef9XsvMzGTevHlVFFVKSgqrVq1qMEXltcLp3vB88803V4q/4447/Ka3q+ejkyZbp8FsKHQHzAuws53jFCA91DzrIwSz6XXcuHF+Nxp7PB6Nj4+vEp+enq7btm0LuJxwkpqaqllZWWqq2uCVz/e4V69eCuicOXN06NChCmhOTk6DyO2FaN/I3AA0pTp1h1AV1SUYf30/O+c9gdmh5FlfIdCGvXz58mqtIrj993lDhw4ddN++fQGVESmGDRtWpWGPcXzEjRkzxm8j956np6dXuuYNbdu2rThOSEhQVdWePXtWSjNw4MCQZbfKyj9NpU7dIVRltcLpTa1wxa0KJc/6CoE07PLycj355JP9KioRqRLXo0cPzcvLq3P+kcZfw3bjbtgiUik+Pj6+4jgxMbHiGNDs7Gzt0qWL3/sXLlyo77//fsiyW2Xln6ZSp+4Q6gR7sarm+hiA8++iJYb5/PPP+eqrr/xeM7/9QQYNGsQ333xDs0ZoKmT37t2VzgcOHMiyZcto06YNcNA8jogwcODACl+Jluglluo01An2XSLSA6OVEZFxwC8hSxVFFBUVMcG7ehljAD65mrSnnnoqixYtiglFNXbs2EqfXnyVb1paWsVxkh8XUW47W2VlZRx77LEALF++nJSUlLDJG1HefBP87DKINRp9nQbTHfMGzOT6MiAP4xl8DXB4KHnWV6jrkOGWW26p6CL/j2Plcz5oqs/Qb8KECVFrxdPfkAFXN9977J2f8DcZW9397msnnHBCxfDXN02wEOlh4N13q4LqmWeqFhSELG990VTq1B1C/tNjPMr0BvoAcaHmV1+hLg173rx5FXNSNziKSkF3gfZxKapHHnmk1rwaCrdCxZl/uOGGGyrFTZo0yW9673yGN943je8137Lc+YYgf2SUVVmZ6s03V9SpZmWpbtoUsrz1QVOqU3cIVVGNAFq6zlsBp9Tx3h7AKcDxQGoocrjy+xZY73x2ryl9bQ17zZo12qxZM00HfcOlqPaCDnIqLiUlRRcsWFDnSrIETkSU1a+/qp5++kFF1a6d6tq1kRDf4odglVWoc1YPAftd5/uBydUlFpFUEblTRLYAM4F7MK7eN4vIpyJySgiyPAc8rao9gKcB/6s3a2Hfvn1cNX48F/Xpw+3FxfwEeGcANgNDgWxg2LBhbN68mRNPPDEEkS31RkEBLFoEf/ubsUX16acmvlcvWLgQrC2xqCfUt4HiaEoAVLXccTZaHV9iPFENUtXfKzIR8QBDgP8RkW6q+kJAQoi0BQYCpzpRbwJPiUgbVd1R/Z0V9/MH4B2gGUbL+WrxGRiDekUtWzJ9yhTOPvvsQES0NARvvQV//jMUFYGPU1nAuHB/8klo1ar+ZbMEjLh0TeA3i3wN/F1VFzvnxwKPquoJ1aRvpqo1vnapSxo/9wwCXlPVvq64NcA4VV3uirsauBogIyNjUNeuXQMpxtIAZGdnq6pWOwKwdRp71Fan1RFqz+omYLqIrAYEM9F+bi3lFQOISIaq7vJNEKiiCgSnx/YCwODBg3XZsmWRKsoSJkRkeU3XbZ3GHrXVaXWEpKxU9VsR6YOZJBfgG1Xd4y+tiDwJdBKRNap6C2a+6s+hlO9iC9BBROJUtcwZimY58RaLpREQstUFRznNBr4CikSk6iozQ0tVPQdYICK3h1qujwy/Ays5OBc+FrMFqNb5KovFEhuEpKxE5FwR+QUoAA4Auc6nP4oBVHUm5sXamaGU7YcJwF9FZD3wV+fc0sgRkatFZJmILNuxwz6bGjOhzlk9CFwALFLV8lrSvuI9UNVXRSSsLUtV1wLHhjNPS/TjO2fVwOJYIkioymq3qn5Tl4SqutDn/JMQy7ZYLE2IUJXVByJyDfA2UOiNVFU/i1oOIiI9gVuBbm4ZVPWYEOWxWCyNlFCV1X3O59OYLSjifNa0MBTgXczi0FeBshBlsFgsTYBQly4EO0FfqqoPhVK2JTRKSkrIz8/3G4qLiykpKaG0tLQieM9LSkooKSmhqKioIq6srAxVpby8vGIfV1mZeQa547xp3HHeY2/8cccdx6WXXtqwP44lKmko7zafichpqvpZA5Xf6CgtLWXNmjX85z//Yc2aNWzatInt27ezd+9e9u7dy4EDB8jPz6eoqKhCkUQjr7zyilVWFr+EpKxE5CjMBuKjgAr3wqpa2zBwNvChiJQDRTjDR1VtG4o8TYGtW7cyb948Zs+ezffff09OTg67d++m1I/n50DoDmQASUALIMEJ8a5P9/EvwHsBlnENcBhmvYy4gvf8VWBFiN/D4mL1ati/3+yLzM+HkhIoLTXBe+z+7NEDAvVr+eCDsH27yy6JE8rLzedf/xq2TeKh9qyeAW4DHgFOA/5C9eus3LwAXAYsJ4A5KxGZDJwHdAX6qeoqJ34TZoLfO8k/SVVn1TXfaEVVWbduHTNmzOD999/n+++/J9/fhlyHdKAzZtLwPwGWNQWzDaGuzCFwZTUWY7WiOr4B1idXZ4e1CaIKO3bAli3QogXauzf79u3jwIED5OXlkZubS15eXqVQWFhISUkJxcXFXH7//bT2MVtcE9906cKTU6dWGpr7Dt19j19csIBONbTJG77+mj7XXssVV1wR8s8RqrJqrqpzRMSjqtuA20RkLvCvWu7brarTgihvOvA48LWfa+d7lVcso6qsWLGCl156iXfeeYedO3dWSZMIDAIGYyweekOGc302B81PePF4PBV+5rzzRd5jgBpf3/ohmIZTgFk17LbmVu46LoaQe4ixSHl5Odt++onds2ZRsmgRcevXk7plC+127SLZ+T3eT0xkdFlZQL/PSKB1AHL8mpPDWzk5Acm+wymjujrN/u47WmzcGFCe1RGqsvL+crudIeEvQJc63DddRCZgrLLUecmDqi4AKjl3bCwUFhby6quv8uijj7J+/fpK1+KA44AzgJMxtnBqsvLeyU9ceXnNa3b/iTFOlo9RKiVOKK3mM5jd5n+sQxopLKw9UQxTWFhIdnY2SxcupOTLL8lauZI+O3bQr7ycDjXc16aoKGBPLNdiHmz5QFlCAuXx8UhCAiQkQHw8ZR4P5R4P5XFxlIlQGhdH97g4PB5PxcOttuOJLmer7k/3cbgsYYSqrN4WkQzg/4AFmP/VHXW4717n8xkCW/JQE/8W8wstAG5R1b2+CdzmRDp37hxCUeGjsLCQJ598koceeojqtou8Cfx3NffvBVZjjN+vBXKAYJ5j3wZxTyTo0KGmv2zsUV5eztKlS/n444+ZM2cO2dnZFBcXMxcYVs09u4B18fFsSU3l99atyW/ThsJOnfhH9+5kZGTQsmVLkpOTSU5OJiUlpeI4OTmZFi1akJCQQEJCAs2aNSMhIYG4uFD+VtFDqEsXHnEOPxORdMywsNY5qxCWPFTHUFXdIiKJwGMY66Pj/JQbNVszVJUZM2Zw3XXXkVNL13seB5XVMswY+FtgEeE3K+HxeEhKSiIpKYnExESaNWtGYmIiCQkJlc6bN29ecRwfH4/H4yEuLo4458nsPfae+0sTHx+PiBAXF1fxtO7Xr1+Yv1H9o6p8++23TJkyhenTp/P7779XSTMXo6zKgJz27dnXty96/PEkn3IKWYMGcYLLA43FEOrbwAWqOgRAVUuAEndcDfeNAJaq6j7nvBUwUFW/DEYOVd3ifBaJyDMYw55Ry/bt27nqqqv46KOP6pT+A8xcz6dA1WZfPSJCVlYWXbp0oUOHDrRt25a2bdvSpk2bis+WLVuSmppaERITExvlMLs+yMvL46WXXuLpp5+uMpQHGDBgACeddBInnngiQ9u1Q3/5hbjTTuOw9PQGkDb2CHUYWMkcjGNHqi6//EOYqRcvXtvtA/0nrx4RSQbiVXWfMwwcgzEXE5XMmjWLiy++uNohnz+2Yt7WVUfz5s3p27cv/fr1o1+/fvTt25fDDz+cTp06kZiYWMOdlnCQm5vLo48+yuOPP86uXQftScbFxXHaaadx9tln86c//anRDXHrm6CUlYjciLES2lJE3A/7JODfdckiQNvt3nKfwFgibQ/MFpFdwCjgPef+OMz0TbiM+oUNVeXhhx/mxhtvDDmvli1bcvLJJzNkyBCGDBnCgAEDYsKxamOjvLyc1157jVtuuYVt27ZVxB9++OFcc801jBs3jnbt2jWghI2MYFziAC0xa50+xrz984bWdbz/a+BY1/mxGCujIbnjCiQE5RAzSEpLS/XKK6+s4oMtkNChQwe9/vrr9auvvtKSkpJ6k72hIdJOToPkxx9/1OOPP75SHfXv31/feecdLS0trTc5YpFA6tQdQvrDA22AZq7zBKBNHe47HtiGWRI0GzPSOS4UWQIN9dWwi4uLdfTo0RUNuhPo56CH1kFBJSQk6OjRo/WLL77QsrKyepE32og2ZVVeXq7PPfecJiUlVdTTSW3a6NY+fbRs69aIl98YaChltRhIcp0nYwzx1eXe1sCfMMuH6tQjC2eoj4ZdUlKi5557bkWjHgb6u7Mh4T+gSdUoqeTkZL311lv1999/j7iM0U40Kav8/Hy98MILK9XV/xs1SstTU81faehQ1eLiiMrQGAhWWYU6wZ6oroWcqponIs2rSywi/wI+VFWvY4mZIZYftagqEyZM4P333wdgNMYmToJz/SuqLqxMSEjgxhtvZNKkSaTZV9dRxW+//cY555zD4sWLAejYsSNf/Nd/0eupp8zjB+APf4BGsqYpKglGw3kD5q1bG9d5W+C7GtJfhLFl9RPwImZyvHkoMgQbIv0UvvXWWyuevleBljk9qgLQ8X56U+ecc47+8ssvEZUpFiEKelY//fSTdu7c+eCwb+hQzb31Vq3YtpuSojp9ekTKbowEUqfuEKqyuhxYj9nMfJtzfFkd7muGGQI+D2zALCW6FMgMRZ5AQiSV1VtvvVXRsK9z7UXfDzrcR0m1bt1aZ86cGTFZYp2GVlbr16/XDh06VNTXFZdfrqU33XRQUWVmqi5bFvZyGzMNoqxMuQzHOI54CDgpyDyOAx4AVoUqT11DpJTV8uXLNTEx0TRsl6LaCXq0j6I6/fTTdd++fRGRo7HQkMpq/fr1mpWVVVFft956q5bfd99BRdWhg+qaNWEtsynQYMoqlAC0As4BjqrvsiOhrHbs2KGdOnVSQM93Df12gfZzKSmPx6OPPvqolpeXh12GxkZDKavffvtNu3btWlFnt99+u5Y/80xlRfXzz2ErrykRrLIKdbvNUqcyK6HVOH4QkanAg6r6ncYCZ0IAACAASURBVLOX8D+Y1euZInKrqr4YijwNiapy+eWXs2XLFoZjVsZ6MNtkTge+d9K1aNGCL7/8kuOOO66BJLXURl5eHqNGjWLTpk0A3Hbbbdx15JHw384OzYwM+OILOPTQhhOyCRLq28C/u46bY+yrba0h/UBV/c45vhj4QVVHikhHzALTmFVWzzzzDB999BFdgWmYSbkiTLdxiZOmbdu2LFu2jE6d/BlxsUQDZWVlXHTRRSxduhSAyy67jLvPOw9OPNH0qVJS4LPPoHfvBpa06RGq1YV57nMR+Rz4vIZb3MaKhmAm1lHVX0SkSg8tVvj+++/53//9X8AYI5uH2RN0JcaiJkCnTp1YtWqVXZIQ5dx99918+OGHAIwYMYLnn38e2b8fjjsOvvoK3nkHBg9uYCmbJuE21ZKGMbNdLSKSJSItMBPzbmVX7fqsaKagoIDRo0dTUlICQB5wPsZK41QnTYcOHVizZo1VVFHOrFmzuOeeewDo1asX06ZNIyEhwQz7PvvMDP1OP72BpWy6hHPOyoNRVA/XcMv/YdZmFQMLVHWNk89xwOZQZGkobrrpJn744YdKcQp84RwfcsghrF27lpSUlHqXzVJ3tmzZwkUXXYSqkpSUxHvvvUfLli0PJkhIMIs+LQ1GOOesSoGNqlrtnJWqvisiX2OsJrh9GmwGrgpRlnpn1qxZPPXUU9Veb9euHevWrbOKKsopLi5m9OjRFeZdXnjhBfr06dPAUll8CWkYqKrzXGFhTYrKdc9vqrrSeYXpjduqqjHVs9q1axcXXXRRtdczMzNZt24dqamp9ShVVbz2sN0G9W6++eZK8XfccYff9C1atKiIz8rKqnStJi87scYtt9zCt98aw84TJkyosV6jgSZbp8Gsd8DMI/9eXQgmz/oOoazJKS8vr7RB2TdkZGTonj17gs4/nKSmplYsbPTildP3uFevXgronDlzdOjQoQpoTk5ORbr09HTNy8tTQBMTE+tFfiK8zmrmzJkVv8HAgQO1oKAgbLJHiqZUp+4QrLLq4oS7gacxTk77A08ANwWTZ32HUJTVa6+9Vq2iat26te7atSvovCPBsGHDqjTsMWPGqKrqmDFj/DZy73l6enqlNKqqHo+nyj3ekJqaWulcVTUnJ6dKuroSSWW1bds2bdOmjQKakpKi69evD+j+hqSp1Kk7hPSHB+bVJS4aQ7DKatOmTRofH+9XUbVq1Up37NgRVL6RxF/DduNupCJSKT4+Pl7bt29f6Z60tLQqDTsvL6/idznssMN0x44dCmjnzp31kEMOqXT/3XffXWfZI6WsysrKdMSIERXyv/baa3W+NxpoKnXqDqEuXcgSkUzviXN8SG03iUhPEXldRL4RkSXeEKIsEaeoqIiRI0f6dTSZmprKunXryMzM9HNnbGPaV/WICElJSRxzjNm4sGHDhorfYevWrVx11VUV6Vq0aME111wTWYHrwOTJk5k9ezYAF198MRdffHEDS1S/xGKdhqqsHgNWisjzIvI8sAJ4tA73vQt8h7HUcKMrRDXjx4+nw/r1VRyMpqSksG7dOtq2bdsgcgXD2LFjK3168W3EaWlpnHzyyZXicnNzK517/dIlJCTgj7vuuosXXniB+Ph4CgsLadOmTUiyh8q3337LrbfeCkC3bt147txzoaysQWUKB42+ToPpjrkDcCTwV4wD2H51vGd5qOWGGgIdBj744IN6AWgp6EegiU5XOS0tLeotevobMuDT5VdV7dmzZ42TsRkZGVUmY3GGFdWVEx8frxMnTtSHH35YVbXCLlRdIczDwC1btmi7du0UjNnoHK+5l3HjVGPItn1TqVN3aBBFAdwPnNYQZXtDIMrqqaee0ktBS0AVNBe0vzNHtXv37jrn0xB4G643iIjecMMNleImTZrkN7377ZB3jsMb8vLyKtLX1rC7dOlSbb51kD9syiovL08HDRpUIcfC884zfwFQTU9X/emnOsvVkDSlOnWHUJXOCRh37VsxyxZ2UIelC8ApwAFgXyD31ZJnD4yj4vXOZ/ea0tdFWZWVlemYUaP0SW+DBs3HGNA77LDDNDc3t84VZAmOcCmrAwcO6PDhw82LENDF3bsfVFStWqlmZ0f8u1gMwSqrUFewvwTcg/FkHsig/wXgMmB5gPfVxHPA06o6VUTGYayQnhJMRgf27+eKkSPpsHgxDwJeGwk7MVYUMs46iw3OZldLdKOqLFy4kKuvuoritWv5B3BTs2a0+vFHk6BjR/j4YzjqqAaV01I7oSqrAlV9I4j7dqvqtBDLrkBE2mK8OZ/qRL0JPCUibVS1VtfHIsJQjH2aZhiLgO/4pFkAXCLCo9Onc9ZZZ4VLdEuEuP3229k3ZQrX/vorncrKWAxU7CUodlx1nHYavPwyHFLrC2xLFCCmVxbkzSL3AgtV9dMA77sF2I3RCRVmY9TlKSfA/AYBr6lqX1fcGmCcqi53xV0NXA2QkZExqGvXrsEUZ6lHss3wrNq31rZOY4/a6rQ6QlVWO4AMzPxTESCYybga3+GLSLnrVF33BeXHqK7Kys3gwYN12bJlwRRnqUdEJFtV62RAytZpbBBInboJdRgYlBWyYLRqLWwBOohInKqWiUgckOXEWyyWRkCoVhdygF+BFIw35l+dOL84PaBqEZFEEekVhBy/Y+xkeVfDjQVW1GW+ymKxxAahGt8bDLzHwSFgvIicV93QC7hFRJKANzCu57djLIT2BE7D+BK8AVgbhDgTgCkicjuwBxgfRB4WiyVKCXUY+DjGqemXACJyMvAkcKK/xKp6nogcg5kQvQPoCORjnL98AAxV1QPBCKKqa4Fjg7nXEru4J9g7d+7cwNJYIkmoyirZq6gAVPUrEUmu6QZVXcJBhy8WS0io6guYdXsMHjw4+LdFlqgn1InufKc3BYCIDMP0lCwWiyWshNqzuhZ4T0SKMEsQEoHzQpbKYrFYfAhVWbUCjgbaYibYtwNHhCqUxWKx+BKqsnoI42X5dwAR8QCTMVtfLDGOqlJWVkZJSUmVUFpaSnl5OeXl5RUbTb3Hvp91ueYtr3379vTo0aOBv3njRFUpLS2lqKiI4uJiCgsLyc3N5cCBAxw4cID8/HwKCgooLCykqKioop69oaSkpKI9lJaWUlZWVhHvbQvu+vUen3LKKVx44YUhyx+qshJ1LYFX1XJnQaalASkrK2Pbtm38+uuv7Nq1i99//52cnBx+/vlnfv31V3bu3Elubm5F4ywqKqJ1SQnJZWU0VyUJSHBCvOvTffw7B30j1pULMZvCPZhuuDd4z98D1ng8lDUCQ3iRory8nN9++41ff/2Vbdu2sXHjxkp1u2PHDg4cOEBBQQHpxcUklZXRAirq1Lce3Z85wNcBynMlkEnl+nTX6WvAlClTokJZHRCRY1V1MYCIHItxShwRRGQyZk6sK8bQ3yonfhNmj6F3n+EkVZ0VKTmiAVXll19+4fvvv2fevHl89dVXbP75Z1rs3YunrIyfA8zvS+C4ANLPIXBlNQEYWsP19cCq8vIaUjR+VJWdO3eyYsUKPvzwQ5YvXkxJTg4HDhzgx+JiNIDtcfOpxT26D28SuLL6O2aRZHUsADb6MQMeDKEqq5uA6SKy2jnvA5wbYp41MR2ztsvfb3q+V3k1RkpLS1m8eDHvvvsuH06fTkJODkdjfvDjgMuBwzEVOguzwjYQAn2FG0zDqe1vJrVcbyyUlZXx448/snbtWtauXcuiRYtYungxbXbu5MjSUvoAfYG/YZ7KccDrBL7KOVrq1OMJz+66kJSVqn4rIn2A4x25vlHVPWGRzH95C4BKzh0bMwUFBXzwwQc8PHkyeStWcAYwAmO4via3FJ1quFYdDwEvYxp4AVDihNJqPoPpPo90Pt0mJstdxwAjR470c2dsU1BQwLx585g7dy5LlixhyZIl5OXlcSRmy8Y1wKuYt1XVEWidigh3AKkeDwUeDyVxcZTFxSEJCZCQgCQk4ElMJC4xEU9iIhoXR0nz5vyhRQvi4+OJi4ur9OkO3jiPx8Pz5eV4PB7ECYiA63iACH8dMiS4H86HUHtWOMppZhhkCZV/i9FiC4BbVHWvb4JYWe28cuVKrrvuOubPnw/AFKp/qpYBPwFrMHuUcsDvENDj8dC2bVsyMjJo3bo1aWlppKWl0bJlS9LS0khJSaFFixY0b96cxMREv5/NmjUjLi6uUvB4PLXGeTwe06B9PAm7z70hPj7kJhkV7Nmzh3feeYcZM2bw1VdfUVBQUOn6p1Tf+y0BfuRgnW52Pj0eD61ataJ79+6ccMIJ9O3blz59+tCnTx9atmwZwW8THTSOlmG26WwRkUSMx52ngHG+iaJ9tfOcOXMYN24cv/32W6X4FRxUVpswY+BvMeZZ12A2ZrZq1YqePXvSrVs3TujWjXGHHUbHjh1p164d7du3p3Xr1mHrjlv8o6rMmzePp59+mhkzZlDsNfLnhxUcVFbrgYXAN5gNs2sBjY/n0EMP5fzzz+faCy+kd+/eFR5nmiqNQlmp6hbns0hEngFmNLBIAbFy5UpGjhzJjh3+jUR8gJm3+BhYByQnJzNs2DDOOvpo7hw4kEGDBpGVldVkhsfRhqoyffp0HnjgAZYsqdtOsjcwjgs+ATY6ca1ateLss8/m3X/8g549a5q2bprEvLJy9iLGq+o+Zxg4BmMuJurZv38/Z555Jl9/XfM7mBzgiyOPZOy55zJy5EgGDx5crT83S/2yfPlyrr32WhYuXFgRJyK1vrVb5YS4uDhOHzmSyZMn06dPn8gKG+PElLISkScwbxvbA7NFZBcwCrPlJw7TAVkD/LnhpKwbM2fO5Mwzz6yxUXfo0IErrriCMWPG0Lt373qUzlIbubm53HTTTTz33HMVddisWTOK67i8oHnz5kyaNImbb76Z5s2bR1rcxkEwLnEaQwjUyWm4KCsr01GjRlXx/eYOxx9/vL777rtaEkNONyMFYXZyGg4WL16s3bp1q6iv+Ph4TUpKqrFOvSEuLk7vuOMOLSgoqBdZo5FA6tQd7IxrPZKfn8+RhxzC2I8+wp8/lQEDBjBz5kwWLlzI+eef32jejDUmXnnlFYYOHcpPP/0EQNeuXTmktJSX8vOp7X3cueeeyy+//MKdd95pe1NBYP8N9cTmzZs577DD+KSsjC5AF+BkoBhIS0vjkUce4bLLLrNv7KIUVeWmm25i8uTJAMTHx9O7d29aff890zA7+VOBszBrx9y0bNmSN954gz/96U/1K3QjwyqreuCHH37gf/v0YR5mjxaYN0AeYMyYMTzxxBO0adOm4QS01EhZWRlXX301L7/8MgCZmZlkZmbS5/vveR2zrw4O1qlbWZ155pm89tprtG7dun6FboTYx3iEWbVqFbf36cMMjKIqwxiZvyw+nmdfeYU33njDKqoopqysjPHjx1coqkMPPZQWLVpw0tq1vIFRVEWY7U5/xazw9zJ58mRmzJhhFVWYsD2rCLJ69Woe6tePtzCvKYuA0cDKLl34z6ef2jd8UY6q8te//pU33jBOx3v37k1eXh7nb9nCI06aA5jX0fNc9yUnJ/PJJ58wbNiwepa4cWOVVYTYvn07dxxxBG9jFFUucA6Qf/zxZH/0ERkZGQ0roKVW7r33Xp599lkAevToQV5eHqdt3lyhqHZhVqG73aoecsghzJkzxz6IIoAdBkaA/Px8Lj3kEN7AKKo84I9AxgUX8NVXX1lFFQN88MEH3H777YDZRyoiHL95M88613cBw6msqPr27Ut2drZVVBHCKqswo6qclpHBNFWaYYZ+5wCHX3wxb7zxBomJiQ0soaU2Vq9ezfjxZjdmWloa7du3p+O6dbyO+cMcwPSo3PaIBg0axDfffMMhh/hblGIJB1ZZhZk/HnkkbxYWkoyZTB8LdLrsMl599dUmvxE1Fti7dy9nn302ubm5iAj9+/dn95IlvIuZTC/ELE9w96iOPPJIZs+eTVpaWoPI3FSwyiqM3HbbbSxftYp1zvn1gOe883jxxRft+qkYYeLEiWzYsAGAESNGMH/+fAqADc71S4C5rvQ9e/Zkzpw5tGpVkzUqS1gIZtl7Ywjh3pqxcOHCg9svQMeADujfXwsLC8NaTlODetxu884771TU4dFHH11pm0xz0PN9ts506dJFt27dGlKZTZFA6tQd7NvAMHDgwAFOPPHEivNS4Mu2bVn31Vd2jipG2LZtGxMmTACgdevWrF+/vtL1QmCa6zwrK4v58+fbOap6xI5NQkRVqyzqTEhIYNmyZXZoECOoKldccQW7d+8GzKT6vn37qk2fmZnJ/Pnzo9rabGPEKqsQOeKIIygqKqoUt2zZMjp1CsYSevjxNSUMcPPNN1eKv+OOO/ymb9GiRUW817ifN+TnG3cEX375ZaX44cOH19t3CxfPP/88n376KWA2Jufk5FSbNi0tjfnz53P44YfXl3hVaLJ1GszYsTGEcMxZ/f3vf69iAuSLL74IOd9wkpqaqllZWWqq2uCV1fe4V69eCuicOXN06NChCmhOTk5FuvT0dM3Ly1NAExMTq9zv8XgqlRMOiPCc1Y8//lhh3qVVq1Y1mndJSkrSlStXhu27BUtTqlN3aHCl0VAhVGU1d+7cKo352WefDSnPSDFs2LAqDXvMmDGqqjpmzBi/jdx7np6eXimNauUG7E2jqvrDDz8ooHl5eZqSklLl9xGRiuNjjz1WVVXj4+MrpfElksqqpKREjz/+eL/y+YZmzZrpkiVLAso/kjSVOnWHBlcaDRVCUVa7du2qUmlXXnll0PlFGn8N2427kYpIpfj4+Hht3759pXvS0tIq3TNq1KhK90ycOLGiYY8ZM0bnzJlT7ZMf0HPOOUdVVe++++4qskdSWd13330VsiQmJlarqOLj4/Xrr78OKO9I01Tq1B3snFWAqGqV7TJHHnkkL7zwQgNJFHlM+6o77jegb775JqeccgpAxQuH5OTkSumnT5+OiNQ4qR1usrOzufPOOwEjr++8oxePx8Nnn33GkDD5vosWYrJOg9FwjSEE27M6rlkz9bieuunp6VpaWhpUXvVFQw4Z3Hm1b99eVbXKtTZt2lSUfcMNN1SSnQj0rA4cOKDdu3ev6HUMqqZH5fF4dPbs2XXKs75pKnXqDg2uNBoqBKOszkpM1ELQ10HjnO50LCz69New/TXmnj17KlQ/GZuRkVHnydi6NmzvPIc3TefOnSvJHglldfnll1fIPQa0FPRRP3Mx0aqoVJtOnbpDgyuNhgqBKqtz4+M1H1RBi0AHg+7fvz+gPBoC396CiOgNN9xQKW7SpEl+03sbr6pWzHF4Q15enqqqzpw5s1L8CSecoKp1b9i+8nnzdd0XVmX1zDPPVJQ1HrTEqdMDoN04OJm+YMGCgH7n+qQp1ak7NLjSCFcAemAcFa93PrvXlL6uyuo/y5bpHaBlTqMuBj07RhRVYyCcymratGkaFxenzUEfcepTQfNAhzl/rJYtW+q6desi+p2aOsEqq8Y0wf4c8LSq9gCeBp4PKbMHH+QSj4e4wYO5E7N6NhfjtPCtggJSU1NDlddSTxQVFXHPPffwP+efz+VlZawG/te5tgtja2wecOyxx7J582Z69OjRYLJaqqdR7A0UkbbAQOBUJ+pN4CkRaaOq/n2yV76f44DJQDOgNXAlMMGVZhXwl1atmLt7t3XTHgOICKcDt2Hq9CLgVipv2VjsxG+Ki+ORhx7i+uuvt3UbxTQKZQV0An5V1TIAVS0Tka1OfK3KCowbpRP9xP8OPAwMePVV5l1ySZjEtdQH6cAJfuK3AA9gut7nXXAB37/6aqVtKJboRMwQMrYRkUHAa6ra1xW3BhinqstdcVcDVzunPaHC9JSXTGBnhMW15QRGF1Wt1v2PrdOYLKPGOq2OxqKs2mIm1jOcXlUcZjqie12Gga58lqnq4EjJacupfxrbb1Af5URrfTaKCXZV/R1YibEijPO5IhBFZbFYopvGMmcFZj58iojcDuwBxjewPBaLJYw0GmWlqmuBY0PMpr42+Nly6o/G9hvURzlRWZ+NYs7KYrE0fhrFnJXFYmn8WGUFiEgPEflWRNY7n90jXN4dIqIickSE8j9TRFaIyEoR+U5Ezg1TvpNFZKNbdhHJEJGZIrLOKet9EQn4tXS4sXVa53xjpk4bfE9fNATgS8yaLIBxwJcRLGsg8CmQAxwRgfwF84LhCOf8SIwTYU8Y8h6CWWi7yZV/OjDcleYh4CVbp7ZOwx2afM/KtVXnTSfqTWBgJJ4kIpKI2bf4Z8zG2UhRDrR0jlsB21S1PNRMVXWBqm7xidutqnNdUYuALqGWFQq2TutOrNQpNKK3gSEQ8ladALgbmKqqGyO1B01VVUQuAD4UkTzMTqIzIlKYDyLiAa4BZtRHeTVg6zRMRFGd2p5VfSEixwNHA89EuJx44B/A2araBRgFvC0iKZEs1+FJjHGKp+qhrAbH1mn9YpWV2dfawdmig/OZ5cSHk2FAL2CjiGwCOgKzRGRkmMvpD2Sp6kIA5zMP6B3mciohIpOB7sDocAxPQsTWaRiIsjq1ykrraauOqj6gqlmq2lVVuwK/AH9U1c/DWY6Tb0cR6QkgIr2B9sCGMJdTgYjcBwwCzlFV/54X6hFbp6ETbXUKdlEoACLSC5iCMWW1Bxivqr6798Nd5ibgTFVdFYG8LwJuxkzKAtyhqtPDkO8TGPuD7TG78ncBF2DMfa0HCpykG1X1v0ItLxRsndY539ipU6usLBZLLNDkh4EWiyU2sMrKYrHEBFZZhQERudN5cxKJvPs7a2zccStFJGx2eEXkHBE5ppY014nI38NY5lwROTRc+YUbW6dBlRnROrXKqg4461waiv6YCc8KVLW/qhZUkz4YzgGqbdgikgRci1mpHS4eA+4MY34BYes0Buu0off7RGvAbJ24EZgL3OPE3QQsAZYDHwHtnfg7gcnOcT/gayfNGuB6J94DfA5c65z3wewl64jZSTALWAasBl7BOGXJADYDezGv4p9wyZbiHB+N8ZP4nfN5tBPfFfN25z5gBcY2+RA/3/OPwG7M6/GVmLdmvmkuBl5xnV/qfJd3gLXAHOf7fIJ5g/RvDr68uRr4wcn7O6CXE5+A8ceRauvU1mmd6q+hlUK0BqfxTHKdj8MYJfM459cA//bTsFOBROc4xWncvZ3ztsBGYCjwPXCGEy8Y+/He49eACa5GNM2PbClO498MjHDi/+CcN3MatmJepYPxOrWwmu/6KjCxht/iJa88Lpn2AB2d84+dhtsS8yf9ziXTPqCTc5wIJLny+RI4zdaprdO6BLs3sGamuI7PAgYDy509YPGYSvMlCXhWRI7CrInJAo4CflDV30XkckyFPqGqnzj3eIC/i8jpQBxmbVB+HeTrCRSr6mwAVZ0jIsVO/AEgV1U/dtIuwngVC4aOmMbrZqGq/uIcrwA2qeo+ABH5D9ANmI35rq+IyIfAJ6r6syuP35y86xNbp4aYq1M7Z1Uzua5jAe5VM7fQX1WPUFV/rgbvx1TYAFU9CjPEaO66PgDTlXdX6IUYUx1DVbUfZq+Z+57qEPzv9PfGuVcelxH8xvUCP/IU+uTte+4t61zgFiAZ+Mr583ppzsFFh/WFrVNDzNWpVVZ1ZwbwZxFpDcY0iPOk9aUVsEVVSx1jZkO9F5y3MxMxT+U2IjLBdc9OVT0gIi0xDd3Lfg6aBvFlLZAoIic7+Z+MmTdYH+B3q6kMMMObngHm6Z3EPkxVl6jqA5g5kQGuJL2B/wSabxixdRogDVmnVlnVEVV9HTPJOE9EvgOy8e/E+V7gKhFZivFePh9ARFoBbwCXqtm7dhHwDxHpj5nPSBWR1cC7mMlcL3OAZBH5j7M1wi1TMXAecL8j0/3A+U58ILwOXOi8PvfnFeh9zKRtoMQBr4rI984w4hCMI2REpIvzHcK+NaWu2DqNrTq1220sdUJEZgE3q+qKMOX3f8BPqvpSOPKzBE6s1antWVnqykTMUzRcbMW8zrc0HDFVp7ZnZbFYYgLbs7JYLDGBVVYWiyUmsMrKYrHEBFZZWSyWmMAqK4vFEhNYZWWxWGICq6wsFktMYJVVFCMir4rIvQ0tR1NERC4VkQURzP9TEbkkxDwiKmO0YZVVDCAiw0VEReR9n/ijnPi5rriznf1g+0Vkp4jMEZGuzrU7nfTX+uRzvRN/Z+S/jQVAVU9X1Sm1pwwOEenq1Olyn/hMESkW4zbMGzdERL4RkX0isltEForI0c61S518HvHJ5xwn/tVIfQdfrLKKHXYAJ4hIhivuEly78UWkG2YD7Q2YHfeHYkyTuL3prnfuczOewHf1W2KDZMdShJcLMcYCARCRNIxdqyeBdKADcBeVTdFsAEb7mIKu9zZjlVUUISIDRGS5iBwQkbepbG+oGJgOjHHSxmHseP/blaY/xhnlHDUcUNX3VHWzK81SIElE+jr59AVaOPFNDhHpJCLvi8gOEdklIk9Vk+5xEdni9FizRaSSmRgRWeZc2+7thYhIcxGZ6uS7V0SWikg759pcEbnSlcdVIvKDU/drRGSgE3+ziGxwxQfqaPR1Kj+cxmMeaF56AKjqm6papqoFqvq5qn7nSvMbxqTMHx2Z0oETMCZ26g2rrKIEEWmGUUavY55w72JMhbh5DdPYwDSc1ZjNo16WA71E5FEROVlEUqop7nVXPpdQufE2GRyF/zHGbnpXTK/irWqSL8U8DNIxZmHeFRHvw+Rx4HFVTQMOx9gxB/PbtgQ6YWyvT8CPYToR+W+MGeXxQBrGguku5/IGjP2slpgez1QRCWTz8VRgjIjEiXE7nwosdl1fD5SJyBQROV0c215+cLe9McCHVO59RRyrrKKH4zBG1h5T1RJVnYZPb0dVvwHSRaQnVZ+QOOZlh2P+dO8AO51Jel+lNRUYKyIJmIY3NQLfJxY4BmOi+EZVzVPVQlX1O2GtqlNVdZeqlqrqwxjb417jdSVANxHJTkCaIAAAIABJREFUVNVcVV3kis8Aujm9lmxV3e8n+yuBB1V1qdMj/klVc5xy31XVraparqpvAz9Sg9caP/yCcSwxAj8PJkeeIRhLpP8P2CEiM7w9QBcfAMPFGBKs0vbqA6usoocs4FetbAYjx0+61zGmPU7GNKBKqOoiVb1AVdtgnsgnAbf6pNkM/IQx7Pajqm4Jz1eIOToBOapaWltCEbnBGabtE5G9mJ5OpnP5Csxwaq0z1DvTiX8d4+HmLRHZKiIPOg8If3JsqKbc8c4Lk71OuUe4yq0rr2EcQozFz4NJVX9Q1UtVtaOTfxbGrZY7TQHG081tQKaqLgxQhpCxyip62AZ0EDGeCxw6+0n3OvBnYKaq1uiAQFWXYixCHuHnsncivkkOAR22AJ2lFh+CzvzUJMwcYWtVbYVxLCEAqvqjqo7FeLr5FzBNRJKdHvJdqtoHM8dzJgeHUr5yHO6n3C6Y3s5EjKecVsAqb7kB8B5wBvCzt8dWHaq6FuMZp6Y283qA5YcFq6yih2+BUuBaEYkXkXPx091X1Y3AMHx6S1DxCvoqEWnrnPfCzH8s8k0LvA2M5OD8SlNkCeYh8YCIJDsT4v7MGqdi6mYHEC8it2PmlgAQkXEi0kZVyzH+AMHMA50sIv2cubH9mGFhmZ/8X8R4whkkhm6OokrGDM92OOVchn8lUiOqmgecghluVkJEejm9xo7OeSdMD8xfm5kHnIp5c1jvWGUVJTg2ts/loP+20Zhekb+0C1R1q59LezHK6XsRyQU+wwwVH/STR4GqztbwegGOKVS1DBiFcTG1GTO/M9pP0lnAp5jJ6ByM1xf30Pk0YLXzmz8OjFHVQqA9MA2jqH7A/Nn9DcPexTgufQPjbms6kK6qazCutr4FtmOcrQY1/FLVZarqb6h5ADgWWCwieRgltQrTg/LNQ503zbuDkSFUrKVQi8USE9ielcViiQmssrJYLDGBVVYWiyUmsMrKYrHEBDWuL2nMZGZmateuXRtaDEstZGdn73QWuNaKrdPYIJA6ddNolJWITMbspesK9KvNhXXXrl1ZtmxZfYhmCQERqXERoxtbp7FBIHXqptEoK8zalMeBrxtaEF927txJXl5eQ4sRUzRv3px27Xy3p1maMo1GWXk3oFberVILr75qQiDMnRtQ8mUTJyJTArOx9j89e9aeyMWZO3cyateu2hPGSBkfZWTw22mn8dxzzwV0P1AvdWrLCKKMSy81IQQajbKqCyJyNXA1QOfOnWHTJpg3L2LlqSob5sxhdG5uxMoAyCouZlAjKiM7NZXf6nhPfdcpYMsIpozhw0POqkkpK1V9AXgBYPDgwUrXrjBsWMTKW7duHT8UFJCdkkJmZiapqal1uu++++4LqJyszz9n9xdfBHRPNJdx/Kmn0n+0v10vVanvOgWwZQRRRhhefER0u42IDFPVCKvuKmVuAs6sbYJ98ODBGunJ2Mcee4ypU81WsE8++cTOwQSBiGSr6uC6pK2POrUEzg8//MAHH3zA9ddfT1JSUkB16ibSPatHgYH+LojIn2q6UVVnRkSieqK8vJzPP/8cgAEDBlhFZWmSFBYW8s9//pNNmzaRnZ3Nm2++GXRekVZWNc123+h8NgeOxth4BrOz/FsgIGUlIk9grBa0B2aLyC5V7RuYuOFj5f9v77zDq6jSBv47tyY3vZCEQEIogdAiICgCCihFRUBdXRBUZC279sIqImxQEFSU/RRRd5W1ImIBG3aaKNio0gJEEyAJpPdyb+695/tj7ow3kIQkJCTB+T3PfZI7mTnzzpzJO+e85y07d5KdnQ3A2LFjW0oMHZ0W5fnnnyctLQ2AcePGYbFYGt1Wc3uwv1/bH6SUI6WUI1EyVg6VUvaXUvZHSVK2r6EnklLeI6XsKKU0SSmjWlJRAXz11VcAGAwGLrnkkpYURUenRfjxxx959913AUhMTGTatNMqk9i8Iysp5cJ67NZTSqklsJdS/iyEeLEZxWo2du/ezYcffojT6WTTpk0ADB48mJCQ2nLw6+icXXz++ef8+KOSt0/96evry7x58zAajafVdpMrKyHEEJRkb1087QuUvF0RtRziFEJcL6Vc7jl+KkpWxjaF2+1m5syZ2tRPRZ8C6vxZ2LJlC0lJSSdtnzFjBh07djzt9ptjZPU/YD5KxsGaUrieyHTgLSHEMpRinHs4uQhnq2fHjh2aogoPD8dqtZKQkMDo0aNbWDIdneanqKiIefPmAWA2m4mIiEAIwfDhw5k4cWKTnKM5lFWFlHJFfXeWUu4HBgohAjzfS5pBpmZHtVEZjUbeffddgoKCWlgiHZ0zg5SSBQsWkJubC8CsWbOYMGFCk5+nOQzsnwshLqvvzp4E+TcDs6WUJUKIOM9Uss1QVVXF2rVrAbjgggt0RaXzp+Krr75i/fr1AAwfPpzx48c3y3maQ1n9HfjMU18tWyhlubPr2P/fwCXAlZ7vJZxQs6y189NPP1FcrNSu1G1UOn82Vq1aBUBISAizZ89uWHxuA2iOaWBDPVNHAv1RSp8jpczzKsvdJlCngFarleHNHb6go9OKcDgc7N27F4DRo0cTGhrabOdqcmUlpTzsKRrZA6Xm2cFTVLytlFJKVRsLIQw0vIhji1FZWclGTyT68OHDsdlsLSuQjs4ZJDk5GYfDAcA555zTrOdqDteFgSgVYO0oSsckhPiLlHJ7LYfs9rgrCCFEHDCLVpiTqjbWrl1LRYVSek+fAur82di5c6f2e79+/Zr1XM1hs3oOmC6l7C6ljAf+Rt0VXB8ARgDtgZ88Mj3UDHI1Obm5ufzf//0fAMHBwVxwwQUtLJGOzpll165dAERGRjZ7/Gtz2Kz8pJTr1S9Syg1CCL869hdSyluBW7UNQgTWsX+rQErJvHnzKCoqAmDmzJmnFfeko9PWkFJqyqq5p4DQPCOrciHESPWLEGI4UF7H/hvrua1VsWrVKrZs2QLAZZddpjt/6vzpOHLkCIWFhUDzTwGheUZW9wIfCCHsKAZ2K0ohh2p4jPAWwCCE8OUPo3oQ0Kqt1EeOHOHZZxXvisjISB56qE3MWnV0mhR1VAVtVFlJKX8RQnRDWQ0UQLKUsqqGXWcDcz2/e1dTKAYWN7VcTYXT6WTOnDlUVlYC8Oijj9Y7A6iOztmEaly32Wx069at2c/X5NNAIcQowCal3COl3A34CSEuPnE/KeVjUkoD8KKU0uD1CZZSzm9quZqKV199lX37lAw2U6dOZdCgQS0skY5Oy6Aqq8TERAyG5q+X3BzTwKepnh20GHiGWjKGAv8RQvhJKcsAPMb4OCnl3maQ7bTYs2cPy5YtA6Br167ceeedLSyRjs6ZZevWrXz11Vc4nU6OHDkCnJkpIDSPshLSK7G7lNIthKgrkc0bwGCv71WebQ3O0dycVFRUkJSUhNvtxmw2M3/+fH31T+dPxZEjR7jvvvs0E4jKmVgJhOZRViVCiPPVhHpCiPOpbpM6EaO3TUtK6fAY31sVzz77rPYmuf322+nevXsLS6Sjc+Y40VYbEaGkpxswYADnnnvuGZGhOZTCQ8BHQgh1GtcLJTd6bVQJIbpIKX8HEEJ0pX55sM4Ymzdv1oI1BwwYwPXXX9/CEunonFm8bbXXX38999133xmXoTlWA38QQvQCLkBZDdwipSyo45DHgM1CiM883y/Hy0G0pSkuLuaxxx4DlFWPRx999IwYE3V0WhKn08n8+fP59ddfAcjIyAAUW+0dd9zRIjI1y3TLo5zqVZ1GSrnG4zg6CkW5PSmlTGkOuRrDhg0byM/PB+DBBx8kOjq6hSXS0Wl+li1bxmeffVZtm8lkalFbbWuxDR0Hfqwj2LnFUJdn/fz8GDduXAtLo6PT/OzevZtXX30VgOjoaM01YfTo0S1qq21xZeUpdvpfFDtVnCdrw1wpZfOkG2wgqrI655xz9OmfzlmFlJL8/Hy8q7JXVVVVW/VevHgx8fHxLSjlH7S4skKxWQ0CvgCQUm71GNlbnPz8fI4ePQqcueVZHZ0zQVVVFbfffnu1FC8ncscdd7QaRQXNWORUCJEohNgohPhFCFFnoicp5fETNtmbS66GcKZjn3R0zhQvv/xynYrq3HPPZerUqWdQolPTZCMrIYRRSuntcjAbuAklmHkN8FUth5YIISI9+yGEGAEUNpVcp4PamUajkd69W7TAs45Ok7Fz507eeOMNALp3784111xT7e8Wi4WRI0e2OrNHU04D1wshZkkpt3i+CxQ7lKDuNMUPo0wBOwshNgLxQNPX8WkE6sgqISEBH582lRZeR4fi4mJ++eUXXK4/xhBSSl544QXcbjcWi4XHH3+cLl26tKCU9acpldVVwNNCiGnATGAh8DZKupeZtR3kKRc/EhjCH35ZLT6yqqysZP/+/YBur9JpezgcDv7xj39w8ODBWve5++6724yigiZUVlLKfOBmIcRFwMfAy1LKi+p5bJEQYp0qjxDCJqWsK2Ffs7Nv3z7tjaTbq3TaGi+99FKdimrYsGFMmjTpDEp0+jSlzcoAXAo4gDHAQ0KIT4H7pJS/1XHc1cASlBzsoIyuJFBX8HNN7XRHCYAOA/KAG6WUhxp6HSrexkd9ZKXTlti2bRvLly8HoFevXjz22GPVavkZjUY6dOjQ6mxSp6Ipp4HvAUUo076rpZR3CCHigeeEED/VkaNqEfBXFKdQ92mc/z/AC1LK5UKI61F8t07Ko6WSl5fHk08+WWtjP//8MwAdO3YkLCzsNMTSOVPk5OTU2ad/FjZt2oSUEqvVyvz58+nUqVNLi9QkNKWy6ial7AcghNgJ4BnZXOEptVUb+V5G+UYhhIhAyZelJkJ/B1gqhGgnpcyp6Zji4mI++OCDU7atj6raDkVFRfXq0z8L991331mjqKBplVWqEGIZysjqZ+8/SCnfruO4D4UQtwPvAlqinAbarGKADNV1QkrpEkJkerZrykoIcRtwG4C/vz8hISF1NhoSEtLqfE10quPdpzab7ZR9+mdACMGIESNOcklo6whvV/vTakjJQTUWJXneN7KeDQshvKd+Eo/NSkpZb5uVEOJc4E0pZW+vbfuA62uLNxw4cKDcunVrfU+h00IIIbZJKeuViFHv07ZBQ/rUm6ZcDXQCn51yx5OPawor31Ggg+qY6slMGu3ZrqOjcxbQGmIDTxspZbbHTnYdsNzzc0dt9iqAtLQ0Bg6sn3J3Op2YTGfFrWqL1Ja7/yRO7FO73c7vv/9ebZ+uXbs2KMVJYWEhx44dAyA+Pr7JnoOysjIt82xsbCx+ftXrAJeXl5OVlVUthbDFYsHhcGjfjUYjwcHBhIaGtrXns9596k2busJT8A/gDSFEElAA3FjXznFxcdRnypCZmcmECRPo1asXo0ePZvTo0URFRTWNxLWwfv165syZwwUXXMBTTz3V1h7EJkUIUe+0QSf26c8//3xSorhHHnmEq6+uK3Ftdd59912efvppANatW0dQUFC9j62Ljz76iMcffxyATz75pMY8aW63m7feeosXX3xR8/kbO3YsO3bsIDs7W9svNDSUBQsWtJlKSw3pU2/alqNFHUgpk6WU50spu3t+HmiKdteuXQsoTqILFy4kMTGROXPm1Lq/213d+6K0tJQ33nhDy7h4Kn7//XeSkpJwOBx8++23LFq0iKayK/7ZyM3N1X5XfYq2bdsGwPHjxykvP/Uajt3+R0y9xWLBbrfz/vvvc8899/DBBx9UC2VpCGrmTYPBQGRkZI37GAwGpk2bxuuvv0779oob4g8//MA777zDggULSExMBJTsIHfccQcvv/zySc/f2cSf95VdT0aMGEF5eTkrV64kOTkZKSVLlizhyiuvrDblyMrKYsGCBWzfvp05c+Zw6aWXApCUlMSmTZsAuO6667jrrruwWq01nqu8vJwHH3yw2tB/9erVxMXFMWXKlGa8yrOTnJw/rAADBw7k559/ZuvWrSxevJh33nmH6OhoPvjggxqnhW63G4PBwLp16zh06BBms5kFCxbw888/a5ljt2zZwrvvvsvUqVNxOBwUFhYSERFBv3796NSpE0II3G43O3fu5LPPPuPHH39k8ODB3HnnnWRkZCClxMfHh2PHjtGxY8dar6Nnz57cf//9PPTQQxQXF/PJJ59www03MGbMGD799FOeeuop7HY7L7/8MsePH2fOnDltzuGzPjTZamBboyErR9u3b+fuu++moKCAI0eOIKWkd+/e/PDDDxgMBr744gsWLVpEaWkpAAEBAXz00UccPXqUadOmkZOTo9kVwsPDCQwMxOl0EhgYSHR0NOHh4TgcDvbv38/evUqdjUmTJrFhwways7MRQvDss88ydOjQZrsfrZXTWQ185plnWLlyJb6+vtx5550888wzJx3z6quvaiMUUAJ933nnHV555RUGDx7MypUrtSlXz549tf1MJhNOpxO32619cnJy8PX1JSQkBD8/P4QQOBwOqqqqFyT38/OjsLCQ1NRUfHx8iIuLY/Lkyfz973/Hz88Ph8OBEAKz2awd43a7mTRpEqmpqYSFhfHpp59qSva3337jwQcf1GxgV1xxBUlJSa1WYbX4auDZyq+//sq9996L3W7H39+fzp078/vvv7N//37mzZtHenp6tbxXACUlJbzyyiukpqaSkZFBaWkpRqORoKAgcnNzq01P9uzZQ1VVFUIIjEYjQgiGDBnCjBkzmDBhAjfffDOVlZXMmTOH5cuX06FDhzN9C9osOTk52mfVqlXaaMmb3bt3a8rK4XDw5JNP8sknnwDwxRdfkJubi5QSX19fgoODiYyMZOrUqYwcOZIVK1Ywc+ZMysvLMRgMSCkpLi7G6XQCnBTi0rVrVw4ePEhZWRmpqam4XC7MZjNut5sVK1awZs0aTCYT+fn5BAYGsmTJEvr06QMoU8Lp06eTlJREXl4eH3/8Mddeey2gLBosW7aM22+/nd9++401a9bgcDiYO3duraP4tkibGVkJIcKAt4CuKMn5UoC/SylzGhMXWN+R1X//+19eeeUVhBDMnz+fgIAArrjiChwOB8HBwfj7+1NRUUFwcDALFy5kzZo1bN26FYPBgMvl4tixYxgMBqKiovDx8aFfv374+vpqKWU3bNhATk4OBoMBo9FIly5deO+99+jaVUmW+swzz7Bo0SKcTidWq5UhQ4YwfPhwhg4dSlRUFDabjcDAQIzGBoVSthkaO7KqrKxk+PDhmsHd19eXiIgIAgICiI2NZcOGDRiNRm655RYWL15McXExDzzwgBYTKoSgoqKCgwcPIqWkb9++bN9e3S78yy+/cMMNN3D06FF8fX3p1KkTDoeD8vJyEhMTufDCCzGZTERHRzNq1CiCg4PZvHkzTz31FOvWrcPPz49rr70Wu91eYyK8wMBAli1bpmVGcLlcXHXVVWRmZuLn50fPnj1p164dw4cPZ+TIkRQXF3PHHXdw6JDy6Hfv3p2nn3661b3gGjuyakvKKhRIlFJu9Hx/GgiVUt4shFgPvOoVF/g3KWWtcYFQf2UlpeTll18mMjKSK6+8EoBrrrmGL774QpWL0NBQwsLCMJvN9OjRg19//VUbwpvNZoYMGcK3334LKFOJkSNHsmTJEo4fP05oaCg2m63aOc877zyWLl3K0aNHmTx5Munp6ZqdRJ06qm/t8vJyhBAkJSW1uSj6+tAYZbV3715mz57NunXrKCsrQ0qJwWAgMDAQi8VCaWkpDocDt9tNVFQUP/zwAzNnztT+yfv27ctjjz3GP/7xD9auXYvL5aJbt24cOFB9zWbJkiXMnj1bM2o/+OCD/PLLL2RlZQHQo0cPLr74YgYOHEi7du0ICwvDarWSnJzMlClTMBgMLFy4kNGjR/Pll1+yYcMG/P39MZvNWp3Kdu3a8b///U9bLVy9ejULFy486drbt2/Pddddx6WXXsrDDz+sKdbAwED+9a9/MXLkyMbc/mbhrFdWJyKE+AtwOzAFOAiEeTmE5gHxdflZNcbb2eFw8NZbb9GpUyfuuOMO7HY7ISEhJ7kWOBwOzGYzQgiuu+467rvvPubMmcM333wDKP5ABQUF2tQvPj6eyy+/nIyMDJKTk3G73Vx22WWkpaVx8OBBhBBUVlaSn5+P3W5HSklUVBT5+fnk5ORo/4xxcXEsWrSIq666qkb5XS4XO3fuZMOGDRw5coTg4GDCw8M555xzuPDCC1uljaOhyuree+9l6dKluFwudu/erdmLjEaj9gHlXlRVVWE2m+nZs6e2fdSoUdx///0sX76c0tJSnnjiCaqqqoiKitJW8FTGjh3L999/D0BUVBRTp07l1ltv5bbbbqtm3Pe6Fi666CIGDx7MU089BcAbb7xRYxbaN998kyVLlgCKIlq6dCmdOnXC7Xbz2muvsX//fvLy8khOTtZsXACjR49m7ty5vPTSS7z99h9RbpdddhkPPvgggYGB9bmVzUqrUVZCiEnAF1LKYiHEPOA8YLaUclsTnsMAfA18AmymnqE23nFksbGx5x4+fLje56ysrGTKlCkcOXKEwYMHM378eJKTk+nfvz8DBw4kJSWF5cuXs379epxOJxaLhY4dO7Js2TJCQ0NxOp08+eST/O9//9OM9F5yYbPZCA4OJi8vT3tTCyEIDAxk8uTJzJ07lzvuuIOUlBSKiorIzMzUDLwul0trTwiB1WolOjqaiIgIzGYzZrOZkJAQUlNTKSkpweVy4Xa7MRqNmEwmLBYL3bp1Y9q0aYwYMaLJfImagoY82JGRkTImJkb7vm/fPhwOB1JK/Pz8sFgsWCwWhg8fTnJyMnv27EEIQVxcHIGBgUycOJHZs2fzz3/+k02bNlFYWMiRI0eoqqoiICCAlJQUrWx6VVUVMTExlJSUYDKZ6NatGyaTiU8//RQhBCtXruT777/nt99Ozo5UXFyM3W4nPDyclStXav3066+/0rdvX81BdMmSJbz55psABAUF8dxzz2k2LIBFixaxcuVKAK3/jUYjI0aMYOHChWzcuJFHH32UyspKDAYDBoOBoKAgAgICuPrqq5k6dWo1u9qZojUpq91Syr5CiPOA54HngLuklEOa8BwvAB1QytL3p4FxgdC4kdXChQtZvXo1oDwoF19cfaZpt9v561//SkZGBiaTiRdeeIFzzz0XUHxhnnzySVavXs3Ro0dxOp2YzWacTqf2oKmrOyd6KQcFBREaGoqvry9Hjx6lsLCwmkLz8/OjvLz8JB8b7wdRSokQAiEEBoMBq9WKEAKXy6Vtt1qtWCwWIiIi6N27N+PHj+fKK6/E19e3Wrvqec7ESKwhD3ZQUJCMj4+nQ4cOjBo1iocffpjKykp8fHwICQnhnHPO4fnnn6dLly78/e9/Z/ny5dpUcMaMGfTu3Zvt27ezYsUKhBAcOnSI0tJSLQXw22+/rQUHb9iwgSuuuAKASy65hMzMTAA6d+5McnIyBQUFXHnllUybNo2UlBTy8vJ46623yMvLo7KykoKCgmp2rsrKSsxmM+Hh4UyYMIEhQ4aQk5PD6tWr+eGHH/Dz88PX15fIyEgcDgclJSXs2bOHyspKpJRIKXG5XISEhNChQwf69OlDSUkJv//+O0eOHKGiogK3213tpRYWFsbIkSMZNGgQffr0QUrJ4cOHSU9Pp7KyUrPLXn755QwcOLDJ+rs1KavtUsoBQojZQLaU8hV1WxO1/wyQCIyXUto96WHOyDSwqKiIq6++mqKiIiIiIli1atVJ/8hbtmzh1ltv5fDhwwQEBHDJJZcwZMgQ3nvvPY4dO6a9qX18fOjevTsGg4HU1FTN+dBisSClpLy8XFNkJ7791IdTVTw2m43ExERyc3NJTU2tpuxquYdamwaDAZPJpD3IRqOx2kNpMpmIiYmhb9++dOvWjZSUFA4cOIDL5SIiIoL27dszfPhwJkyYQLt27U46lyqr0+msNg07EZfLhcFgOOlaG/Jg9+rVS1511VXMmjWLGTNm8NZbb2G322nXrh1Wq5V27drx2WefERkZyX/+8x/mzJlDRUUF48eP5/bbb+eBBx7g8OHDSCkJCgrSXiqgvDSuueYabSRz2223adOsL7/8kldffZX169dTUFCgOYqqSj0yMpLS0lLt/lZUVFBVVaWNwNW+dLvd2Gw22rVrR0BAgHZdBQUFZGVlYTAYcDqdmr3Nuz+977X6XBgMBs12qo6+VaXmfWxAQAAhISEEBgbWGi0RHR3NNddcwzXXXHOSjRUU52cfH596RVu0JmW1Dfg3MAtFoaQKIfZIKfuc4tD6tL0AJVf7OO8UMp5CE8u8DOw3SynrtCg2NkLfO0xi4sSJBAQEsGnTJs455xxmz54NKAZa1efFarUSEBBAXl6e9nYzGo1ERUWRlJTE9OnTefvtt1m6dClSSqZPn87mzZv54IMPMJlMCCGoqqrSlt3VhzE4OBij0Uh+fr6maPz8/JBSUlFRUatntWonU9tS+99qtWIymaiqqsJoNOJwOLTpovfDX5NCUQkICKB9+/Z07dpVWwnNy8vTjnG5XJrC8vZBslgsWK1WwsLC6NmzJ71796ZPnz4MHz68UQZ2u93OoEGDSElJoaqqis6dO2ujjnnz5nH55ZezZcsWJk2aRH5+PjExMQQEBPDbb79RUlKC2WzGYDBo/QWKsoqIiCAlJQWbzUaPHj1IT08nMDCQPXv2MH36dNavX69N6R0OB2VlZdrIFf4YOdvt9mojY4PBoBUVBeUF4uPjg7+/P35+fthsNsrLyyksLKSoqEizW6r7qvezqqrqpGgH9WVkMBjw9/enY8eOWK1Wdu3aVa0d7xF3aGgoPj4+uN1uCgsLMRqNmnz+/v4MGjSIAQMG0LdvX3766Sfef/99Dh8+jNlsJiAggNDQUCIiIggPD9emnX5+fiQmJjJx4sTmVVZCiEuAnlLKpUIpmxUkpawxwbMQ4gKUAhEbpJTPebKF3iOlvPsU5/BFSW1cIaU8VsPfewN7UEZRFZ7NqVLKq4R5bFGUAAAgAElEQVQQCSiuCyF44gJPFW7TGGX1448/smzZMo4cOaKtznkzZswY4uPjefrpp7UAWu9pHvzhGT1ixAg+/vhj/P39qaqqYsqUKZpNKSUlpdo/SUBAAEII7HY7vr6+2lL2/v37NUfU00V96M1ms/YWVm1iJ+KtvLzf1EKIatd6Oi9C1fWjocpq7dq1/POf/2Tt2rVkZWUhpaRXr17aKGbChAkkJSVRUlJCv379SE9P114gTqdTu4bo6GitwK2qbHx9fXn11VcZMGAAAwYMwO12M3jwYPz8/EhPTycnJweTyYTNZtOmkd4Ooep987Yxel+vr6+vJoP64rHb7dWeIXWk5z0VV5+HiooKsrKytH288R5Nq22d6cU1m81GWVlZ8ykrIcTDwOVAeyllvBCiI7BSSjmscSJXa9uAEnB8C9APJS2yD0pOrA+B/6tNKZ4uAwcOlFvvugtef71e+5dXVHAgOZlb4+Ox2+24XC58fHzIz8+nvLwcl8tFQEAABoOBnJwcLSL+0qwspnimeFJKLWGXn58ffn5+dOvWDbPZTEFhIXv37sXpdKIOCb0fMBXvob1qr5iGUqCxITR0IftMnuN1lDePZxRZ7wc7MTFR3mIyMfTQIRwOB06XC6TEZrNpoxOL1UqPHj04duyYNs0bCZrytVgsBAUF4XQ6NYdQVYFZrVYGDBjADS4XvT0vOh8fH4wGA46qKqTbjcnj5On0jIYlgJRc7HkZqMpdbVcdtahtqQsm1zud3Fi/gYTyE7jUx0cbpVVUVNR9IG2jT72pjwf7dcBAPNk/pZTpQoha1z89I6nXgA5Sys5CiAHABCnlozXsvhn4AXgA2CY9mT49dqixwMtCiP9IKVc24JrqT1oaePyfToUNxZIPiu+Uuqrm7++Pr68vx44dq7ZcHRMTw8iRI/nrvn1ctHHjyQ2WlSkfTyhHCKBqfyEE/v7+9O3bl71792qxglJK7WFWbT9Go5FEm40RJSUNunTVtcLlctUrGDcOGNGgMzQc9RwbG3l8WloaPv7+nHviaLO0VPmov2/ZQifAO+GvyWSiR48edOrUiczMTKqqqjT7k6pUnE4nu3btYpzDwUXqiPPEYOgTQmu821en1BaLBYPBQGhoKKGhoeTn5+N0OqmsrNSmiJ2krN/9PmFVWZ3ae5sMvKf73nQRghHNPLqKo3qfns5orj7KqkJKWXXCG76uM74EPA6omft3onieP1rDvuOllLknbpRSZnuOeUsIEV4PGRtHXBwMH17nLlJKDh06RIlHGVx88cV8/PHH1RSTGtKg2gAiIyNZtGiRsjL08MPsDAqiuKREeXgAg9GoKAjPSAs8oyjP+fz9/QkODiYsLIzXXnsNPz8/vv/+ez766COysrIICQmhsFAprWiz2XCEh7PN4wMUEhoKQHFRkfIGBwxCEBYejs3Xl8LCQoqLi0ns0QNQphPHjx8nOzu72pRPHWUYjUbsdjtHnE42NvBBMxqNmud4586d6dmzJ+PGjSMkJESzwRQXF3Po0CF++uknTKmp7MjJwRIVxbQLLmjQuUCZpmVaLGTGx1NUXKytmvpYrXTv0YODHqdOh8NB1QnTqdDQUD7++GM6d+7Mhx9+yBNPPEFAQAAFBQXa9Nhut1NZWclvUvKt0mmYTCbtHqnKwWwy4VRHT579VDuhwWAgPDycCy+8kPvvv5/IyEjCw8O5//77+eWXXygqKlIWYior2eR2g0e5mUwmJfDZaqV9+/ZKnq3jx3GfUMDUYrHgcrk05+HJkydjNBp5/fXXSUtLo7KyUlOah4FvvUb7qP/jNfWzZ5Qf6O9f7aWp2jW9nx1vhZSGoqiOGgzEd+16Wq4S9ZkGvo/ifrAEZYT1CNBLSlljGgAhxC9SykFCiB1Syv6ebdrvpxRIGVV1kVL+2IDraDD1tVl5+7qMGTOGBQsWsGLFCu655x7cbre2cuN2uzly5AglJSW0b9++2spYbm6uZv8ICgrC7XZrq0HecWTqNMBgMGh+O6B4IRcXFyOl1Jah1elEREQE7dq1QwjBX/7yF2bNmgX8UenlW6+RY+fOnTn33HNJTU0lJSWFnJwcCgoKNDtLVVUVNputxhUd9R+htLRUU5SdO3emS5cuWK1WqqqqMBgM9O3bl9jYWHr16kV8fPxpL3c3ZMowYMAA+e9//5sRI0Ywbtw4du3ahdvtJiEhga+//pqJEyeSnJxMXl6epnxUD3cfHx+WLFnCLbfcAsCHH37I/fffz+HDh7XRiaqY1BGplFIzxpeXl2urfeHh4URFRWEymQgKCiItLY3jx48DaEbn0tJSLbA9ISGB888/n6NHj/LJJ59QXFysGdItFgtms5nIyEj69+/P3XffzZ49e5g1axalpaXatM8bo9GIv78/1157LZ07d+btt9+uZvPKz8+nuLi42oKP0WjU/NDUBRbVzGC1WrXr9PQJZrMZk8mkLQKpLzf1HvXo0YPOnTtz/vnnM3DgwGoJAptzGng38CbQBygHvgPqqqLgEkKY8Yy+hBAdgDqT7AghvgOuQFHwO4BCIcTnUsoH6yFfs3LFFVewfv16/Pz8SEpKQgjB1KlTsdvtpKenExcXR2xsLCkpKSxatAiXy0Vubi5Wq5XAwEDtQVYf+A4dOhAcHExRURFBQUGEh4eTkpJCdnY2paWl2krfFVdcwZYtWyguLqa4uBhQRg5BQUGUlJRowc+qMTcxMZF//vOfmtzt2rXjmWeeYcOGDTz77LNkZmaSmppKamqqto/FYqmWS0ltv2/fvkRHR1NZWUllZSVxcXEMGzaMhIQEnE4nhYWF+Pr6Vltebw2oCxdut5vc3FzNqB4eHo7JZKJTp07aip3b7SY2Npbk5GRt8SIpKYlBgwbx008/sXjxYtLT0zW7oWpbUlfFQkNDyc3N1UZU2qjKbKa0tJTg4GAef/xx3n//fdLS0oiNjcVsNhMTE8OmTZs0BZORkcGuXbv48MMP8fX1pbS0FJfLRVBQEMHBwdrIJT8/n7Vr17J+/XqqqqooKyvDaDRy2WWXMXbsWN555x3279+vKVF/f382b97M5s2bAWX0P378eEpLS8nIyODAgQPY7XZKSkq0azoxW6nVaqVjx46EhoZqspSXl2ujMyklISEhJCQkkJCQQEREhOYT2BwB1PV2XRBC2ACDlLLO5SchxI0odQATgVdRDOizpZTv1HHMDillf4/bQSLwMLBLStm3fpfRcBqyGlhcXExlZaXmvVwbhw4d4uabb9aM7jfeeCOJiYlMmzYNt9tNv379eOmll+jVq5f2lnG73cybN4+lS5dSVFSE0Whk8uTJjBkzhj59+vDNN9+wc+dOIiIi2L59O9nZ2Rw7dozy8nKCg4Pp0qULQ4YMYdasWbXWN3Q4HKxcuZLVq1dTVlaG2+3GarXStWtX4uPjiYuLo0OHDnTo0IHIyMgW8Wqujca4LuTn5zNmzBjS0tIIDAxk7NixTJw4kWuvvVabAnXq1ImEhAR+/PFH8vPzNd+00NBQ/P39yc3N1V4yQ4cOJS8vj4yMDOx2u+Z2oSoc9X/I39+fwMBAAgMDueqqqzh06BA7duzQHEaHDRumRSgUFxdTVlamjbDVVUOr1aoFqIMyEjp+/LhmhlBRPdXfeecdTckUFxezceNG1q1bx+7du7WXXM+ePXn88cerleVyOp0cOnQIp9NJ7969ycjIICcnR3OTiYmJITo6ulkcf5vbdaErSrYDbSQmpfy8jv2HAeNRRkqfSim/O0X7+6SUvYQQ/wPelVJ+3ZCpY2Nojkoomzdv5vXXX2fr1q2aATw9PZ2SkhIMBgP/+te/MBgMJCQkkJiYyNatW3n22We17JVms5nY2NhqyeC6d++OxWJhz5492raAgACGDh3K+PHj6d+/f4Nyirc1GqOsDh48yJQpU7TwmISEBFJSUlDDq26++WaSkpJYu3YtTz/9NIWFhWRkZGhTWdU+pI5QLrzwQjp16sT27dspLCykoKBAi8lUbTcGg4E777yT/Px8duzYgdvtxm63k5+fr/mudevWTfN3mj9/PsOGDSMjI4P169fzxRdfcODAAc477zwmTJhAYGAgX3/9Nd9//73mwqC6qUgpSUxM5PnnnyfUY6M8ESklR48epbS0lO7du7eq1NjNNg0UQixCWYE8AKjWPAnUqqyklN8D3zdAjo1CiAMoaZZvF0IEe52rTbBmzRrmzZuH2+0mODhYexMGBQXhcDgIDw9n9erV1Xxg1Bg/UKYwI0eOpEuXLuzevVvz4zp48A/PjYiICO666y5GjRp1Viuo00VVJE6nk6KiInbs2KHZVqKjo5k7dy5hYWHs27cPIQSdO3fm6quv5qWXXtJiCdV+Kioq0gLQAwICqgUCe/s6GY1GvvvuO+x2O94xp2pBh9DQUMrKyujcuTOLFy8mNjYWUCp+33jjjdx448klA4YNa7x3kBBCO8fZQn3U7VVAZ3mKoqMeQ3ytwzQp5V/rOPxO4Bzgdymlw+MacWs9ZGsVrFq1iieeeEL7XlJSwnXXXUdwcDAZGRkUFBQQExPDihUrtH3y8/PJysrCbDbTrl07nn32WcaNGwegGdK//vprvvrqKzIyMpg0aRK33XZbjaEOOtXJysrC7XZrU7XIyEjMZjNxcXFER0dr0+Xk5GRAmSbNmTOHQ4cOkZ6ejsvl0qZIqluHlFJLnOjt/6auvkVHR2M0GrHZbPj7+2t2q27dunH++efz2GOPcezYMWJiYqplANWpP/VRVkeAuoPNFNY0VggppRRCHAQ6epxOQUmw1+pZsWIF//73vwG0BzU7O5s1a9awevVqQkJCcLvdmjIzmUwMGTKEjz/+mI4dOxIREcFzzz3HgAF/hE6qNpVbb72VW2+9tZpflc6pOXbsGFlZWVomzqioKAICAjhw4AClpaXk5uZis9lITU3VRsL+/v6sWrWKzMxMMjIyuPHGGykuLkYIQWRkJAaD4aR8VmFhYVr20FdffVWzJYKSZ0zNY6UqNtXTXKdx1EdZzQA+FUJ8Q/Xy7i967ySlfKOxQggh7gEWoITKeE81W33vqkvSgYGBLF26lJycHGbMmEFJSQkPPPAADoeDAwcOkJ6ejsViwd/fn02bNhESEkJwcDBLly4lISGhznPoiqphfPvttxQUFABocYlbtmzRbEw33ngjt956K2VlZRw7doz33nuPw4cP88QTTxATE0N6ejoVFRVUVlZqcW75+fl07dqVmJgYpJQcO3aMfv36cc011zBq1CisVive6Wl0mp76KKuHUWL2+lFdkdSIUMrI/82zv4+6XUr5tzrOcR/QQ0qZWQ95WhX3338/ZrOZsWPH0r17d6SUDB06lM2bN7N7924AbYkY0HyaevbsyVNPPVVjvTidxrNv3z6+/fZbzefJYrFw/PhxbcXOx8eH7OxsZs6cqSXT8/HxYd++fUyZMoW+ffuyceNG8vPzqaysxGazaX5l/v7+vPbaa0RERNSYz12neamPsjoX6C7r7yf/X0+7I1G82acAm05xTHpbVFSgTNnuvrt6jPaDDz7I3r17KSoqomfPnloGgoqKCvz9/Zk4cSIzZ87UjeRNjMvlYubMmZqR3GazUVxcrBXs6Nq1K9OnT2flypXayprZbCYhIYHvv/+ekpISdu3aRVhYmGY8LygoICAggICAACZPnqy5r+iK6sxTH2V1EPAD6hvef55Uku/9KqV8QgjxIvDuKY6ZK4RYhrLC6D3VrHXFsTWybds2lixZwjPPPMMnn3yiBTffddddxMXFERoayksvvaQVg9BpWjIzMzEYDJonvpo4T3XILSoq4u6772bw4MH85S9/wWQyaalh/Pz8NB8om82mpfVRHW/9/f256aabWvoS/9TUR1kVA9uEEF9RXZE8VMv+ari3Swhhk1IWebzY62K859OderpHtDYyMjKYMWMGpaWl3HTTTVpGSrvdriXvHzFihK6omhE1W6rFYtECzYUnfk8IQWlpKUuWLNHCbJxOpxbwW1paqoUNqc65FotFyxt2yy23tIr85X9m6qOskj2f+pIvhAgBvgS+EELkAsdPccxVQJyU8tR5LVopUVFRjBkzhtWrV5OVlcXf/vY3HnroIXJzczXv6AsaEZyrU386d+5M+/bttfJZoNijQkNDyc7OpqqqiqSkJM2tARTv/oiICM34npqaqrksBAQEcO+99zJu3LgaizronFmaI1OoUSrphQ0o9qpglBzpxXUcsw4YK6U8OWtYM9EcHuxSSl5//XVeeOGFk/5mNBpZt24d/v7+TXrOs52GeDuHhYVJNX5N9UZX83+pyshgMGjVYNSAW9V2mJiYyHnnncd///tfHA4HVquV559//k9ZCbs5aXIPdiHEtVLK94UQd9T09xNdF7wIFUIUSSkdwHIhhAU4VbmUQ8B6IcRH1OEe0doRQjB9+nSioqKYN29etSyRY8aM0RVVM6OWN/OO2VPzdqmZLLKzs/Hx8cFisdC1a1eKioqoqqqif//+PPfcc5hMJj7//HPN/tW/f7NFfOk0kLqmgX2A94HRQOEJfwsCalMka1BWAlVHUjPwKTC4jnNZgd8A78DltlnQEKVG2+DBg0lNTcXHxwdfX1/i4uJaWqyzHqMnT5hqUAclk4TqyW40GjUP8lmzZjFhwgQKCgo4dOgQ/fv31zzLH330UR577DEuvfRSPWKgFVGffFYnVaapaZvX33ZKKfudaltL0xzTwBNxOp289957vPTSS8TGxhIREUFUVBQ9evQgISGBTp06nVQdBxRD8datWwkPD6dnz56tKgvCmaYhUwZfX18JaJkD1IR3at6l+Ph4unTpwpw5c0hMTGxewXVqpTmmgSbAAhiEUsxB/Y8JQsnyW5cw7aSnDJYnmV6NTilCiHNlHcVPhRBWlLjEhhj4m42ysjK2bdtGUFAQ8fHx2Gw2Ld3Hzp07+fnnnzl8+DCRkZG0b9+eL7/8km3btpGTk8O+ffvo0KHDSb5V/v7+RERE0KdPHwYOHEh6ejorV67U0ntERUUxbNgwxo4dyznnnKP799SBWqFHTVSoKiyz2cy0adOYMWPGWRfc21w4HA6OHTtGUVERRUVFWkJAdZW1V69e9OvXj+DgYO655x4tc4jq+tG1a1cOHTqExWKhS5cuJCQkcMEFF3Deeec1WqZaR1ZCiLnA3Br+VAwsllLOr+W4v6F4vb+JouBuAJ6QUr5Ww76rUBTfCuAnIAvF670HcClKoYoZUsp1DbusU9PQkdXWrVuZO3cuWVlZwB/11lRfnpqQUpKWlobL5cJqteJ2uwkJCdH+piYwU/NLqUvsFosFIQQlJSVagjdQpjRquXo1p5Kai8lqteLj40NwcDBRUVHExsZqoza1Nl1BQYFWvFJNuGaxWAgNDSUmJoaYmBiCgoIIDAwkKCiIdu3aERYWVi2Rmtvt1hK2+fr6asVSm4uGjqy8y5aZzWa6d+/Oiy++qBvJG8i+fftqzARxIj4+Pmzfvl17RtVnQX22vZ+NgIAAioqKmrW6zVIp5V0NalSI4cA4z9c681kJpXLzbSh55TuiZCPdjVLd5n9SyoZVQqgn9VVWDoeDF198sVpq2Nowm8106tSJ/Px88vPzCQgI4KabbiIvL48VK1aQnZ1NWVmZljIW/uhU77a9O7ipV2ubGjUDgfcHTr6uBiR5RHiqQ3uS39X7wTaZTFqghdlspmPHjrz55psMGdJkxcDPGqSU2nRZ7StfX19t5J6ZmcmECRNqPPbEVMqHDh3SArhrep5VfHx8qKioaN7kew1qUIiVwN9RDOy7gHBgoZTymSY90WlS31JcR44c0YpDSClZcdttWhGJgoICrbM7duyIw+Fg3759+Pv7a2/28Xl5XJaVpRQerWUEdiKtuUzWmTjH6zSubFNAQICcXFnJTaCNLuuVjqWm6kN18frr9S7h1pLnyMzMPKmGoPoykVLyytChTPAEfAMYPdW9bTYbNj8/zRRhMhqV3OomEz8/9RRCCHJycsjIyMDpdJKTk8PWrVtJS0vTXsRqyNPfDAZucLt5UwjWxcSQlpbWrDnYG0oPj9f6NcA6lKwNPwKtSlkB9SrFFev5qPxj0yacTieZmZlakraysjJ27NhRraSV+hIYJyUN7pUGEseZK5N1Js6xsZHHx8bGMqiggKHHjkF6uvJpDhpQwq0lz3HnnXdqhUpq4qrjxxl0YkFUr/TJNZWVeuyxx6p9V3O6qyOqE0vbd5KS4VLyncFw2jn7m0NZqa+y4cAXUspyIUT9hhRnmnqU4gKlwKmvjw+paWnaNrvdjslk0kI81M46kVQa/s/n4+NDQEAAFRUVdZaCB8VAX+nnxz5PCl4fHx8EUOV0UllRQUVlJS6nE7eUOD3lxSUQ36FDNWdJtTiEyWTSqpSYzWbtY8zPZ1tODtLtxqXaI7zkcHuu3eDxaTKaTFw6ZAjh4eG0b9+eiIgIQkNDiYyM1AoKGAwGSkpKOH78OJmZmcRv3szhrVvp368fL4wa1cC7pqy+ZlosbPP3JyEh4aQCCE1GPZ+blj7HqWZNBUFB/OhVY9EtJW6XS+tLAZgtFgxeZgl18UJVTHl5edp3Pz8/rZAJKFPxdCnZ5HRy1FNZ/HRojmngeyhe6wlAT5TKNj+cDa4LmZmZCCHw9fXlgQceIC0tjW3btlVTWOBVJdfLQVFFtcn4+fnh66njpxqE1VJIqrJwu91a6IfRaMTPz48OHTpoBvDJkyczaNAgAgICzlqH04ZMGbp27SrVBYwvv/yS8PDmKznZFti0aZNWIgyqKy/1eVWfLdVxVkrJd999x5o1a7SyWqoy8v64XC6trH1eXh4+Pj5YrVYtk6rNZiMxMbFa+a4BAwawYMGCVmWz8kWpprxLSpnqCWLuK6X8sklPdJqcrp9VZmYm1157Lenp6eTk5CCE0EYmaueD4vOjjoxMJpNWiddisZCTk4PNZiMkJASbzVZjkr2oqCiuvvpqLrroIrqeZpHItkhDHuzY2FgZERGByWRiy5YtupvHabBlyxYeeeQRLZVObcTGxhIYGFitoAnA4sWLGV7LyLDV2Kw8wcgfeX3PADKa+jwnIoTojmKXDQPygBullIdOt93XXnuN+fPna28l1WXBaDRSWVmJEIKIiAhtrq7GowFa8ciqqiqsVqs2alKnVl26dMHPz49HHnmESy+9FFCKHezfv5+0tDS6devGsGHD9Eyh9UQNb4qIiNAV1WkyZMgQPv/8c1JSUgC0bBZWq1V7ls1mM0FBQSQnJ2sFTkCpWdnDU/G7KWk99XlOn/8AL0gplwul/uB/gYtPt9GcnBzNt0rF+3twcDDt27cHaHDWz/DwcBYtWlTNmzo6OlrPHtpI1JWvU9V31Kkf6lTOm8zMzGoVxe12O+Hh4URERFRzXSksLKz2XV0dP53ip2eFsvJ4yQ9AiWMEeAdY6u1J31g6depEfHy8dtNdLpdWC85kMmkVjWNjY7n77rs1T9/CwkLtTS+EIDw8nOjoaNq1a6eV3o6Li9OzhTYh6v32rjKt07Rcf/31mktDQxk9enS1KlAN5axQVkAMkCGldAF4UtRkeraflrKaNGkSkyZNqrbN7XZz+PBhsrOzq6UaOVuN3G0Fg8GAr6+vrqxaKadrb21yA3tLIIQ4FyVnVm+vbfuA66WU27223YbiLQ9KSE/12kotRziQ29JC0Drl6CSlbFefg4QQOcDhU+54ZmiN97IlaVSfenO2KKsIlFzxYZ5RlRHFyB5/utPAM4EQYmtjVkd0OVovreUaziY5zoolEyllNrATuM6z6TpgR1tQVDo6OvXjbLFZAfwDeEMIkYRSLPXUIeM6OjpthrNGWXlyXp3f0nI0kpdbWgAPuhxNR2u5hrNGjrPCZqWjo3P2c1bYrHR0dM5+dGWlo6PTJtCVlY6OTpvgrDGwn20IIYYB01HKlBU2NLV0E8rRC7gXMKI8L9NlCxo6hRBvAFVSyltaSobGovdprfLUr09rylWjf5r2g5IlNRWlFmIfr+3dgR9QHFp/QHFiren4jwH/ViDH+y0pB3AXSpXvZXqf/vn6VJ8Gnhk+Ai7i5FAQNVNEd+AFlEwR1RBCjAP2SynrTizUjHIIIUYKId5GCZcobwk5PCFVNmBLE5y/KdD79DTlaHCftvQb6s/0AdLwvHWACJRK10bPd6Pnezuv/W8CHm1pObyOex4Y0BJyAEnAK8DbwH7gwpbuT71Pz2yf6jarlqPOTBFCiCuAx4E1Qoj/AP+SzRM+dCo5RgDX4EnJDeypraHmlENKOQ9ACBEHzJF1lHdrQfQ+bYAcDe1TXVm1UqSUa1DqKLa0HBtpfMGZJkdKmQa0OeM66H1aG/XtU91m1XIcBTp4MkTg+Rnt2a7L0XJynA6t5RrOSjl0ZdVCyFaSKUKXo+loLddw1srR0gbKP8MHWAKkA07gOLDXsz0B+AllWfcnlAKxuhxnSI6z4Rr+THLogcw6OjptAn0aqKOj0ybQlZWOjk6bQFdWOjo6bQJdWbVihBAbPY6ECCHmCSEmneqYBrb/qBCi1sKFQgiDEOJ7IUST+AYJIfoIIT5virbaInp/nh66smojSCmTpJTvNnGzc4G6qqxei7Kqk94UJ5NS7gEcQoiRTdFeW0bvz4ajK6tTIISwCSHeF0LsE0LsEkK859keJYTYIITYJoTYK4RY5HXMo0KIlUKIz4UQKUKId4UQ/YUQ64UQvwkhnvbad6MQ4lnPzxQhxMJa5HhdCHGXV/vveNpPFkJ8JoSwef4WJIRY5dm+TgjxphDimRrae8Hz6xYhxE4hRHANp70NWHGCrIuFEN8JIY4KIf4phLhOCLFFCJEmhLi2rnvm4R1a0ANd78823J8t7a/S2j/AVcBar+8hnp8+eFJroMRXrQcu9Xx/FDgEBKEEb+4CvkLJY+QHZONJlYES9vA1SuiTP7AbuMLrb+rvrwN3ndB+MEp819fArZ6/LcaTbgMIRZ9rkDwAAAJRSURBVEnb8Uwt1yapJT2I55oqAF+vbRuBd1FectEo0foLPH87D0iv6555fo8Fjuv9qfdnQz/6yOrU7AIShBAveN40ds92I/C0EGIXsA3oA/TzOu4rKWWRVII4fwW+kVLapZRlKJWgu3rt+4aU0imVlCErgYvrIddXUspCqTwxP3m1NxJ4DUBKmY+SuqMxhAMOKWXFCdvfl1K6pZSZKIVkP/Rs34YSWuFD7fcMFIfBSCGEuZFynS56f1anzfSnrqxOgZTyd6An8A0wCtjl6cAHgBDgfCllIspD5ON1aKXX764avtcWRC5Q3pCnorb26nv8qaig+vXUdt5KUCLqPdtMddwzPG1WSSmrmkDGBqP35ynP22r7U1dWp0AoKycuKeVHwP1AO5TheDBwTEpZKYToAEw8jdPcIIQwCSH8UIygG06jrQ3ANI/sIaeQqwRlanMSUspCIEso6TsaRB33DJSH/teGttlU6P0p4hoqQGvpT11ZnZq+wA+e6cHPwBOeIfMSYKgQYgfwErDuNM6xHViLEvT5mVRSiTSWeUCEEGIvsBzYDBTVsu9iYH0dBtkPgbGNkKG2e4anvVWNaLOp0Puz4bSO/jwThjH9U/sHL6NrE7VnBnw8vweivPVGNbKtzsCPeIrhNoFsFpR/4PCWvu96f7a9/tST7519hABfCCV3kA+wQkq5tjENSSlThRCLgfZA5qn2rwedgEeklLlN0NafBb0/PehZF3R0dNoEus1KR0enTaArKx0dnTaBrqx0dHTaBLqy0tHRaRPoykpHR6dN8P+bFUkjkd36AQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x367.2 with 11 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "# set up axes\n",
    "fig = plt.figure(figsize=(4, 5.1))\n",
    "gs = fig.add_gridspec(6, 2, height_ratios=(1, 1, 1, 1.15, 1, 1), wspace=.1, hspace=.0)\n",
    "ax_top = np.array([[fig.add_subplot(gs[row, col]) for col in range(2)] for row in range(3)])\n",
    "ax_placeholder = fig.add_subplot(gs[3, :])  # space holder between upper and other axes\n",
    "ax_placeholder.axis('off')\n",
    "ax_mid = [fig.add_subplot(gs[4, 0]), fig.add_subplot(gs[4, 1])]\n",
    "ax_bottom = [fig.add_subplot(gs[5, n]) for n in range(2)]\n",
    "\n",
    "\n",
    "for trial in range(ntrials):\n",
    "    # relaxation plot: subsample to make plot more clear\n",
    "    for n, index in enumerate(example_relaxation_curve_indices):\n",
    "        label = f'{int(sampling_time_ms[index])}ms'\n",
    "        ax_top[n, 0].plot(times, \n",
    "                       relaxation_arr_dec[:, trial, index]*normalization_factor, \n",
    "                       lw=lw, c=cols[index] if plot_colorscaled else 'k')\n",
    "        \n",
    "        ax_top[n, 1].plot(times, \n",
    "                       relaxation_arr_full[:, trial, index]*normalization_factor, \n",
    "                       lw=lw, c=cols[index] if plot_colorscaled else 'k')\n",
    "        for _ax in ax_top[n]: _ax.text(8, 2, label)\n",
    "    \n",
    "    # implied timescales\n",
    "    for its_n in [0, 3]:\n",
    "        if plot_colorscaled:\n",
    "            plot_segmented(sampling_arr_dec[0, trial], sampling_arr_dec[its_n + 2, trial], cmap, norm=norm, ax=ax_bottom[0], lw=lw, alpha=alpha)\n",
    "            plot_segmented(sampling_arr_full[0, trial], sampling_arr_full[its_n + 2, trial],cmap, norm=norm, ax=ax_bottom[1], lw=lw, alpha=alpha)\n",
    "        else:\n",
    "            ax_bottom[0].plot(sampling_arr_dec[0, trial], sampling_arr_dec[its_n + 2, trial], lw=lw, c='k', alpha=alpha)\n",
    "            ax_bottom[1].plot(sampling_arr_full[0, trial], sampling_arr_full[its_n + 2, trial], lw=lw, c='k', alpha=alpha)\n",
    "        \n",
    "# connectivity\n",
    "if plot_colorscaled:\n",
    "    plot_segmented(sampling_arr_full[0, 0], (sampling_arr_full[1, :] == 1).astype(float).mean(axis=0), cmap, ax=ax_mid[1], norm=norm, lw=lw, alpha=alpha)\n",
    "    plot_segmented(sampling_arr_dec[0, 0], (sampling_arr_dec[1, :] == 1).astype(float).mean(axis=0), cmap, ax=ax_mid[0], norm=norm, lw=lw, alpha=alpha)\n",
    "else:\n",
    "    ax_mid[1].plot(sampling_arr_full[0, 0], (sampling_arr_full[1, :] == 1).astype(float).mean(axis=0), lw=lw, c='k', alpha=alpha)\n",
    "    ax_mid[0].plot(sampling_arr_dec[0, 0], (sampling_arr_dec[1, :] == 1).astype(float).mean(axis=0), lw=lw, c='k', alpha=alpha)\n",
    "    \n",
    "# set labels\n",
    "for _ax in [ax_top[0, 0], ax_mid[0]]:\n",
    "    _ax.set_title('dMSM')\n",
    "for _ax in [ax_top[0, 1], ax_mid[1]]:\n",
    "    _ax.set_title('classical MSM')\n",
    "ax_top[1, 0].set_ylabel('conductance\\n' + r'(m S/cm${}^2$)', va='bottom')\n",
    "ax_mid[0].set_ylabel('% con-\\nnected')\n",
    "ax_bottom[0].set_ylabel('time\\nscales\\n(ms)')\n",
    "#ax_bottom[0].set_ylabel('its/ms')\n",
    "\n",
    "# sharex = True for connectivity and its plots\n",
    "for _ax in np.concatenate([ax_bottom, ax_mid]):\n",
    "    _ax.semilogx()\n",
    "    _ax.set_xlim(0.9*sampling_arr_full[0, 0, 0], \n",
    "                 1.2*sampling_arr_full[0, 0, -1])\n",
    "\n",
    "# relaxation plot: reference, labels, sharey\n",
    "relax_ylabels = [0, 15]\n",
    "ymax_relax = reduce_max_conductance * max(relaxation_arr_dec[-1, :, :][np.isfinite(relaxation_arr_dec[-1, :, :])].max(), relaxation_arr_full[-1, :, :][np.isfinite(relaxation_arr_full[-1, :, :])].max()) * normalization_factor\n",
    "for _nax, _ax in enumerate(ax_top[2]):\n",
    "    _ax.plot(time, g_k*(ode_solution**4), 'r', linestyle='--', dashes=dashes, #label='reference' if _nax == 0 else None, \n",
    "             lw=lw, zorder=100)\n",
    "    _ax.set_xticks(np.arange(0, 16, 4))\n",
    "    _ax.set_xlabel('relaxation t (ms)')\n",
    "    _ax.set_ylim(0, \n",
    "                 ymax_relax)\n",
    "    if _nax == 1:\n",
    "        _ax.set_yticklabels([])\n",
    "    else:\n",
    "        _ax.set_yticks(relax_ylabels)\n",
    "\n",
    "for _nax, _ax in enumerate(ax_top[0].tolist() + ax_top[1].tolist()):\n",
    "    _ax.plot(time, g_k*(ode_solution**4), 'r', linestyle='--', dashes=dashes, #label='reference' if _nax == 0 else None, \n",
    "             lw=lw, zorder=100)\n",
    "    _ax.set_xticks(np.arange(0, 16, 4))\n",
    "    _ax.set_xlabel('relaxation t (ms)')\n",
    "    _ax.set_ylim(0, ymax_relax)\n",
    "    if _nax % 2 == 1:             \n",
    "        _ax.set_yticklabels([])\n",
    "    else:\n",
    "        _ax.set_yticks(relax_ylabels)\n",
    "\n",
    "\n",
    "# connectivity plot: reference, xlabels, ylim\n",
    "for _nax, _ax in enumerate(ax_mid):\n",
    "    _ax.plot(sampling_arr_dec[0, trial], [1 for _ in sampling_arr_dec[0, trial]], 'r', linestyle='--', dashes=dashes,\n",
    "             lw=lw, zorder=100)\n",
    "    _ax.set_xlabel('')\n",
    "    _ax.set_xticks([])\n",
    "    _ax.set_ylim(-0.1, 1.1)\n",
    "    if _nax == 1:\n",
    "        _ax.set_yticklabels([])\n",
    "    \n",
    "# its plots: reference, xlabels, ylim\n",
    "for _nax, _ax in enumerate(ax_bottom):\n",
    "    _ax.set_xlabel(f'sampling t (ms)')\n",
    "    _ax.set_ylim(0, 25)\n",
    "    _ax.set_xticks(sampling_time_ms[example_relaxation_curve_indices])\n",
    "    _ax.plot([sampling_arr_full[0, 0, 0], sampling_arr_full[0, 0, -1]], [its_ref[1], its_ref[1]],\n",
    "                    lw=lw, linestyle='--', zorder=100, color='r', dashes=dashes)\n",
    "    \n",
    "    _ax.plot([sampling_arr_full[0, 0, 0], sampling_arr_full[0, 0, -1]], [its_ref[4], its_ref[4]],\n",
    "                        lw=lw, linestyle='--', zorder=100, color='r', dashes=dashes)\n",
    "    if _nax == 1:\n",
    "        _ax.set_yticklabels([])\n",
    "        _ax.set_yticks([])\n",
    "\n",
    "fig.align_ylabels([ax_top[1, 0], ax_mid[0], ax_bottom[0]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## decomposed MSM example\n",
    "Here we show at example how to compute a full MSM from submodels and to access full system experimental observables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "using 750.0 ms sampling time\n"
     ]
    }
   ],
   "source": [
    "end = 7500\n",
    "print(f'using {end * msmlag * 1000} ms sampling time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100.0 % connected\n"
     ]
    }
   ],
   "source": [
    "# subsystem count matrices\n",
    "singlet_ccmats = [msmtools.estimation.connected_cmatrix(\n",
    "    msmtools.estimation.count_matrix(singlet_traj[:int(end)], 1, sparse_return=False))\n",
    "                   for singlet_traj in singlet_trajs]\n",
    "print(f'{100*np.mean([ccmat.shape[0] / 2 for ccmat in singlet_ccmats])} % connected')\n",
    "\n",
    "# subsystem transition matrices\n",
    "singlet_tmats = [msmtools.estimation.tmatrix(ccmat) for ccmat in singlet_ccmats]       \n",
    "\n",
    "# full system transition matrix on all 16 combinatorial states\n",
    "tmat_full_uncg = np.kron(np.kron(np.kron(singlet_tmats[0], singlet_tmats[1]), singlet_tmats[2]), singlet_tmats[3])\n",
    "\n",
    "# full system transition matrix in empirical state formulation\n",
    "tmat_full = cg_transition_matrix(tmat_full_uncg, chi) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# tmat_full is fully connected\n",
    "msm = pyemma.msm.markov_model(tmat_full, dt_model='0.1 ms')\n",
    "\n",
    "# relaxation experiment (to reproduce electropysiology)\n",
    "times, rel = msm.relaxation(p0, observable, maxtime=15)\n",
    "assert np.all(times == time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f76cffa9358>"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAEgCAYAAAAUr2zFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hU1fbw8e9Kg9AhhN4EE6p0ECyAoCI2VPypCGK9eO2+XK961ateu2C5FlSwYgHkShVRUIpKUQQBRREU6YiEUAIESJn1/nFOdIghmWRmcmaS9XmePDPnzJl9VuCZWdlnr7O3qCrGGGNMJIjxOgBjjDEmjyUlY4wxEcOSkjHGmIhhSckYY0zEiPM6AGOMKY+WL19eJy4u7jWgHeWvg+ADVufk5FzXpUuXnf4vWFIyxhgPxMXFvVavXr3WycnJe2JiYspVGbTP55O0tLQ2O3bseA043/+18padjTEmUrRLTk7OKG8JCSAmJkaTk5P34fQSj37Ng3iMMcZATHlMSHnc3/0vOciSkjHGmIhhSckYYwwNGzY8ITk5uX1OTs4f+5577rkkEeny2GOPJaelpcWed955x6WkpLRNTU1t06pVqzYzZsyoCvD8888niUiXxx9/PDnvvT6fj0aNGp1Qs2bNDsWJw5KSMcYYAJKTk7MnT55cPW/7vffeq922bdtMgBEjRjRs0KBB9tq1a39Yt27djwsWLFjbunXrI3nHtm7dOnPChAlJedsfffRR1Ro1auRQTJaUjDHGAHD55Zenjxs3Lgngp59+Sjh06FBMy5YtDwFs3749vkGDBlkxMU7aqFevXm5KSkpW3nubNGlypEKFCrp8+fKKAG+++WbS5Zdfvqu4MVhSMsYYA8CZZ565f82aNYlpaWmxY8aMqT148OA/ksqtt9668+mnn25wwgkntL7mmmsa51268zdkyJBdr732WtK+fftili1bVuX888/PKG4Mdp+SMcZEABHpEq62VXV5gDHowIED97zxxhu1pk2bVnPJkiU/LV++vDLAwIED92/atOm7WbNmVf3yyy+rDhs2rMWNN96447HHHtuR9/6rr756T4cOHdqkpKQcOeuss/bGx8cXu7rQkpIxxkSAQBNHuA0fPnzXySef3PrEE0/cX69evVz/12rWrOkbMmTIviFDhuzr0qXLwVGjRtX3T0rVq1f3derU6cDDDz/ccM6cOWtLcn5LSsYYY/7Qpk2brHvvvXfbySeffNB//9SpU6v17t37QK1atXw+n48VK1ZUaty48ZH877/vvvt2dO/e/WC3bt0Or127NqG457ekZIwx5ih33HHHXwoUVq5cmXjXXXc1UlUBaNas2eExY8Zszn9cly5dDnfp0uVwSc8ttvKsMcaUvlWrVm3s0KFDsavTypJVq1bV7tChQzP/fVZ9Z4wxJmJYUjLGGBMxLCkZY4yJGJaUjDHGRAxLSsYYYyKGJSVjjDERw5KSMcaYiGFJyRhjTMSwpGSMMYZ58+ZVPvHEE1Pbtm3bum3btq0nTpxYfebMmVWbNm3aLj09PRZg0KBBzW688caGec8vu+yypieddFJqs2bN2l122WVNDx8+LMHGYUnJGGMiQffuLXn+eWeRvCNHhO7dW/LSS7UA2L8/hu7dW/LqqzUBSE+PpXv3lowbVwOA336Lo3v3lowf7yzQt3mzs/3BB9UCOfWuXbtib7rppqaTJk3a8MMPP6yZNWvWL7fddlvTHj16ZF5yySXpQ4YMafriiy8m/frrrxWfe+65bXnv+/bbbyt//PHHv/z888+rt2zZkvDMM8/UDvafwZKSMcaUc3Pnzq2ydevWhP79+6e0atWqTf/+/VNEhB9//LHCE0888dvevXvj/v3vfzeaNGnSr/Hx8X+878ILL9xdvXp1X3x8PEOHDk1fsGBBQEmwMDYhqzHGRIKlS/9c6qFCBT1qu2pV31HbSUm5R23Xr59z1HaTJkdvF0FVadmy5aFly5b95T1paWmx27dvT0hISNCdO3fGpqSkHLMNkaCv3llPyRhjyrt+/fod2LRpU4UPP/zwj9VkP//880o+n4/Bgwc3Gzp06K6xY8duGDp0aPM9e/b8kTemT59eKyMjIyY7O5sJEyYk9e7du9grzeZnPSVjjCnnkpOTcz/44INf/vnPfza+4447YrOzs6VJkyZH+vTps//w4cMxjzzyyI6YmBjmzZu3Z9iwYU0//PDDDQA9evTYf9ZZZ7XYvn17hR49euwfMWJE0LOe29IVxhjjgWhfumLQoEHNunTpcvCee+5JK2kbtnSFMcaYiGaX74wxxhTb5MmTN4aj3ahJSjExMZqYmOh1GMYYExJTpkwhNze3abjP4/P56Nq16/JwnydUoiYpJSYmcvDgQa/DMMaYkFizZg2tWrUKSRl1YZYvX+4L6wlKyOfzCfCX2GxMyRhjPFCxYkXS09Mpj8VmPp9P0tLSqgOr878WNT0lY4wpSxo1asTWrVtJSytx8VpAdu3aJatWrQp6+p8Q8wGrc3Jyrsv/QshLwkUkFRgHJAHpwDBV/TnfMW8D7f12tQcuUNUZx2q3cuXKapfvjDGmeEQkU1Urex1HoMKRlOYBb6jquyIyFLhGVfsWcnwHYB7QQFWPHOs4S0rGGFN80ZaUQjqmJCJ1gM7ABHfXBKCziCQX8rZrgfcKS0jGGGPCQ0SeEpENIqIi0s5vf0UReVlEfhaR70VkbGnEE+oxpcbANlXNBVDVXBHZ7u7/y4VTEUkALgdOL6gxERkODAdISEgIcajGGGOAacBzwJf59o8EDgOpqqoiUrc0gvG60OECYLOqrizoRVUdC4wF5/JdaQZmjDGe8h9aCWPZuKoudE7x5zlEpAowDGik7hiPqv4etiD8hDopbQEaikis20uKBRq4+wtyDfBGiGMwxkQqnw+ys6FCBWc7LQ2ysqBhQ2d7zRpnu0MHVJXc+fPJOXKEwz17kpWVReyUKWSrsve008jKyqLapElkxcezvXdvcnJyaPT++xypUoWNffqQm5tLq3ff5UCtWqzt3Zvc3Fy6v/02u+vU4btevcjJyeHMt95ie8OGLD35ZHJycrh03Dh+adaML7t1Iycnh+vHj+f75s2Z36kTPp+POydO5KvUVOa1b4/P5+OhiROZ17o1n7ZpQ0xODk9Onsystm2ZnZpKQlYWz86YwQdt2jD7+OOpdPgwoz/+mHfatuWTZs2onpnJa3Pn8lLbtnzSsCF1MjN544sveKZNGz6pV48GBw7wztKlPNSyJVfNmUOTJk1K83+qBU6h2gMichpwALgvL4GFU0iTkqruFJGVwGDgXfdxhaoWdOmuEXAqzuU7Y0wk2rMHfv8dWrUCQJcuJXvVKvacfz4ZGRnETJtG7HffsfrSS8nMzKTZtGnU+OUXPnS3T/vwQ+pv28aT55xDZmYmN86fT4OMDK7s3JnMzExGrl1L3awszqlfn6ysLN5JS6Nabi6nxMeTlZXFZ6pUiInhnKpVSUhIYOa+fWTHxHBt06YkJCQwYf169iUk8MC8ecTHx/PiN9+wp3JlXl+/nri4OB758kv2Vq/O7MOHiY2NpdvatWTs2cPqWrWIjY2l7759HKpcmb179xIbG0tWfDwkJFClShViY2PJSUqicv36NG/eHBHhUJMm1Dj+eNq3b09MTAwZK1eS3KkTJ3ftSqwq6Rs30qBrV/p36kRcbi479+yhebduXNyxI3FZWWz3+Wh94onUO+EE4g8fZmO1anQ96SRSWrUiITOTjQ0bcuopp9C5ZUviMzLYMGMGZ/fqRe3aQVV0x4nIMr/tse5VqELfAzTH+f7+p4icCHwoIseratDLUxQmHNV3rXBKwmsCe3BKwteKyCzgflVd5h53L3CCql4WSLtWfWdMCahCZiYkJkJMDPz8M7poEfvOOYf0jAxyp0yh+uTJzLn+enbt3UvLjz6i18KFXH/hhezbv5/LVq3i0m3bSG3ShIz9+7l3zx6uV+W45GSqVavGAxkZnJGRwXWnn05iYiKXrF9Pq/R03rr4YhITEzlp9WrqZGTw7WWXkZiYyPHffUfVw4fZNWgQlSpVovbKlVTIySFrwAASEhKotGYN8SLEnnQSCQkJxG7a5MTdrJnz++zb52xXdZf9ycmB2NiwXt6KdoFW34nIRuBcVV0tIrWB34CEvMt3IvIjzvf5skKaCT7eaLmb2JKSMQXYswdWrEC7dWO/Kns/+ojEl1/my8GD2ZCZSePZs7nk00+55KST+GHvXs7dvJknDxwgtUoVcpKTGabKsPR0Hh8wgMr169M9PZ12mzfzwzXXUKVWLer//ju109PxDRpEtRo1qJqVRQWABg28/s1NgEqSlNztOcBTqjrHvf90MXC8qu4Na7yWlIyJYAcOwNdfo+3asTsujp2zZlH3P//ho/PP55ucHGp//TX3L11K30qV+Bq4pEYN7t+7l+d69YKWLTlBhI7bt7P/0ktJSkmhTlwcNYGElBSI87rOyZSGopKSiDwPXATUA3YB6araVkSa44z5JwHZwL2q+nHY47WkZIzHDh6E3FyoVo2sDRs4ct11rOzdm4Xx8WQuWcLDH37IlZUqMT0+nl716/Of9HTm9etHTPfupCYlkXr4MElnnknNpk3DPrmniT7RdvOsJSVjSlNuLsyejdarx5batfnhyy8ZMHQob7ZrxxNZWezeuJElwP/at2d3nz60bdGCrtnZ1B8wgKTjj/c6ehOFLCmFiSUlE7VGjmRvYiLzGzViyeLF3P/f/zIpNpZ7a9akY8eO3HrkCHFnnknD88+nRYsWVMgrlzYmBCwphYklJRMtfI89xs7165ncuTOLFy/mn5Mn86Mq751+Oj179uS0WrU4/owzqJuS4nWophyItqRkI53GBOuNNzjw4YeMHzCAzz77jPNmzKB2QgIrROjXrx8V77qLy9q14/IYW77MmKJYUjKmuBYsIPf555lz1VVMnzWL4/73P/rs38+SxETOOeccTnvmGRo1asQAr+M0JgpZUjKmKOnpMH48Gf37M+vbb9n6wgtcvGQJb23dSvdLL6X/55/Ttm1bTrTKN2OCZknJmIKkpUFuLvsrV2b+mDGcf++9/KNiRX7r148Lr7qKSlOm8H7dUpk02ZhyxZKSMfnk7t+PNm3Kp02aMHjHDnqdeiqxzz7LM9deS9W86W2MMWFhSckYgGee4eCKFTydksKrr77KtfXqkTJoED/ffjvJyYWtUWmMCSVLSqZ8UoVvvkG7dWP+/PnsHjOGxA0b+P2aa5g5cyYdOnTwOkJjyiVLSqZcynnrLeKuuYahqamsiI3l9hEjuHzIEM6pUsXr0Iwp1+zmWVM+5OTA229zpEEDXv7pJ14eNYrhVavS5rHH6H/BBcTYPUSmjLKbZ42JQIcOHCD7jjuYfuQIn595JuNnzKBLly5eh2WMycf+PDRl16ef4vu//2PMSy9xfNu23N61K+2+/JKpU6daQjImQllPyZRJqsqKuXOpOWMGC7ZvZ4b1jIyJCpaUTNlx+DD84x9sqV+faz7/nK1btvDUpEmMP/98W2fImChhl+9MmXEwJ4dfZ8xg4uOPM3DgQL77/nvOGTjQEpIxhRCRp0Rkg4ioiLQr4PUHjvVaOIQ8KYlIqogsEZF17mOB8/OLyCUi8r2IrHYfbc4WU3wbNsCwYXw8aRJt27XjwVNPZdivv3LzzTcTHx/vdXTGRINpQC9gU/4XRKQz0APYXFrBhOPy3SvAaFV9V0SGAmOAvv4HiEhX4EGgr6ruEJHqwJEwxGLKuH3r1hE3cSLvzJvHa2+9xemnn+51SMZEFVVdCPzlioKIVABGA5cD80srnpAmJRGpA3QGznB3TQBeFJFkVU3zO/T/AU+p6g4AVd0XyjhMGbdtGyxezMdVqnDddddx5fDhvD5qFImJiV5HZkwkihORZX7bY1V1bADvewh4V1U3lOYl8FD3lBoD21Q1F0BVc0Vku7vfPym1ATaIyBdAFWAK8KhGy528xlPZ99xD9sSJ3Fm3Lu+99x59+vTxOiRjIlmOqnYtzhtEpCfQDbg7PCEdm1fVd3FAe5weVQLwCc41y7f9DxKR4cBwgISEhFIO0USUQ4fgyBG+37KFaxctos9ZZ7Hw7bepXr2615EZUxb1BlrhdB4AGgGzReRqVZ0TzhOHdJoh9/LdOiDJ7SXFAulAiv/lOxGZCUxS1bfd7TuBJqp687HatmmGyjGfD049le0ZGXT47TeeefZZrrjiCq+jMiYqBDrNkIhsBM5V1dXFeS3UQlp9p6o7gZXAYHfXYGBFvvEkgPHAmeKIB/oBq0IZiyk7Dmdl8UZ8PA/t3s38BQssIRkTQiLyvIhsxekNfSYiP3gaT6iHcUSkFTAOqAnsAYap6loRmQXcr6rLRCQGeAoYAPiA2cAdquo7VrvWUypnVOHpp9nVqBH9R42iRYsWvP7667bInjHFFG0Tstos4SYyHTzIoQ4dmLxjB1vvu4+77rrLboI1pgSiLSnZjA4mshw4ALm5TP/sM9rv2UPFN9/k7rvvtoRkTDlhc9+ZyJGZCb168W2VKty4fj1TP/6Y7t27ex2VMaYUWVIyEUMTE5lXsSKTNm5kyZIlNGnSxOuQjDGlzJKS8d4PP+CrVImbRo5kWXY2n6xYQVJSktdRGWM8YEnJeCsnB73gAtbt28ea1q2ZO3cu1apV8zoqY4xHLCkZTx3JzeWe+vXZ36gRH8+aZfPXGVPOWVIy3vjmG3JWr+bS6dOJTU5m4sSJttSEMcZKwo03fCNHknbrrUhWFhMmTLCEZIwBrKdkPODz+RheoQKZnToxYcoUm2zXGPMHS0qm9Pz2G/rQQ9yalcXPW7bw8SefULFiRa+jMsZEEEtKpvR8/jlH3nyTPSkpzFy8mEqVKnkdkTEmwgQ8piQiiSLSXETqhzMgU3a9duAAp9avz3/nzbOJVY0xBSq0p+TO5j0MuA7oCOwDKopINjAVeFZV14U9ShO9VOHOO1ncqBH/fuIJvvjiC5KTk72OyhgToYq6fLcIWAKMAJbnLXPuLubXHxgrIq+o6sTwhmmi1q5dHBk/ni/27mXKZ5+RkpLidUTGmAhW1OW781R1hKouzUtI4Czmp6rvqGof4LOwRmii2pbDh+nk85H69tv07NnT63CMMfmIyFMiskFEVETaufuSRGSWiKwVke9EZIqIlMoljkKTkqruKqqBQI4x5dCyZWTffz8XXnABV48YwUWDBnkdkTGmYNOAXsAmv30KjFTVlqraHlgPPFEawRQ1prS0sNdV1dYVMAXS999n78sv07F/f+644w6vwzHGHIOqLgSOWrNMVXcDC/wO+wq4oTTiKWpMqQ3wI/AW4Om67Sa6vNi4MVOaNmXm22/bAn3GRDG34O0GYEZpnK+opFQX+D/gSpwqvLeBCaq6J9yBmSj11lssqVqVRx59lCVLllC5ctSswmxMWRUnIsv8tseq6thivP8F4ADwov9OEekH9AMaAYeAVcB0Vd0WTLCiqoEdKNIMuAqnPPxOVR0fzImLq3Llynrw4MHSPKUprp078bVowes+H82mTeOMM87wOiJjyj0RyVTVIv86FJGNwLmqutpv31NAe5yityPuvsuAB3HGmb4CdgAVgdY4Seor4C5V3VGSeAOa0UFEqgNnAWcCv7jBHOvYVGAckASkA8NU9ed8xzwI3Ahsd3ctUtWbihu8iSw5tWoxLDWVjuedx98sIRkT1UTkUaALcE5eQnJ1Ak45VpGbiJwJnAxMLtF5C+spicgAnEt37YFJwFuqurHQBkXmAW+o6rsiMhS4RlX75jvmQaCKqgY8Am49pQiWkwOLF/PA3LksWrSIOXPmEBNjE9AbEwmK6imJyPPARUA9YBdOZ+ISYDWwDufSHMAGVb0wzOEWmZR8wLc4VRi+/K+r6p35jq+D80skqWquiMTi/IIpqprmd9yDWFIqO158EW65hTOTkhj3/ffUr28zURkTKQK9fFeCduOA4UBfnBLy+TjjVTnBtFvUn7MPAR8C+4GDBfzk1xjYlnejrfu43d2f32XuTVlzRKTAuypFZLiILBORZTk5Qf2eJox2XXABt9SqxYh337WEZEz5MRoYCEzBudfpPPIVQ5REwIUOATUm0gV4W1Xb+u37ERiqqt/67asHpKtqtoicAbwHtFbV9GO1bT2lCHTwIBoby8BLLqFly5aMGjXK64iMMfmEsae0Bmirqj53Ow74XlVbB9NuoT0lEWkjIrX9ti8VkWki8oyIFDTN8xagoXvZDvexgbv/D6q6Q1Wz3eefuq+3C+YXMR64+WbSWrZk57ZtPProo15HY4wpXVuBCn7b8eT7ri+Joi7fveWeCBHpCrwKLAbq4NSuH0VVdwIrgcHursHACv/xJLethn7POwLNgLUl+QWMd37r0YMXd+1i3PjxtnqsMeXPBmCJiPxLRP6FM4H3jyJyo4jcWNJGiyoJr6iqv7nPLwbGqepItwe06hjv+TswTkTuB/bg3HSLiMwC7lfVZcBj7qW+XCALuKKkNe3GGz6fj8vGj2fgww/TsmVLr8MxxpS+eGAFkOpurwKqA91wCh9KpKik5F9xdyLOwBZuZV1uQW9Q1Z/cY/PvP9vv+ZXFD9VEjBtuYO6+feTm5nLbbbd5HY0xxgOqenU42i3q8t0WEblJRAbi3EQ1F0BEKuBe1jPlzKFDZK5YwfJp03jjjTeIjY31OiJjjAdE5DHxm9hSRCqLyHvBtltUUroJGAD8B/i735x3/YCZwZ7cRJ/chAT6x8ZS6ZFHSE1NLfoNxpiyqgHwpYg0dGsDluHM+BOUkJaEh5OVhEeA997j5S1bmDBrFgsWLLBZG4yJAuEqCXfbvgJ4CjgCXKmq84Nts6j1lM4u7HVVnRVsACZKbNiAXnkl+xMSeHXFCktIxpRz7m1B5+NMR9QE6Iozq0NQiip0mAl8jzNVUP5FcRSwpFReHHccd/bqRY2ePa3azhgDsBxn8u3HcKru3hSRj1V1QDCNFpWUHsaZmC8deBP4JO/uXVOOZGQw84svmL51K9/9+99eR2OMiQx/U9XP3ed7gQtFJOhy3IDGlESkD85aSicBU4HRqro52JMXh40peWTzZrR9e26Pi+O8iRM5/fTTvY7IGFMMoR5TEpG4Y026KiL1gr3nNKCBAVVdAFwN/Au4BhgUzElNFKlUiaXNmsEpp1hCMsYALM17IiLT8r0W9JBOkYv8iUgrnIR0AbAEuCQUFRYmOnz/22+cu20bq2fP9joUY0xk8K8vaFrIayVSVPXdV0AszmDWqTjrtCMilQBUNTPYAEyEOnIEve02Hluxgv/85z/UrVvX64iMMZFBj/G8oO1iK+ryXXecmRyeB37DWVdpP05y2h/syU0EW7mSnHHjqLFzJ8OHD/c6GmNMmIjIUyKyQURURNr57U8VkSUiss59THFfihGRRLdz8sfzvO1g4ym0AVWN8fuJ9fuJUVWbX6YMO9yhAyfWqcP/vf46cXFFXuU1xkSvaUAvYFO+/a/gFLWl4sx7Osbd3x6nY3LAfX6QPzsrJwQbjH3bmL/6+Wf+O3kyTTt3pm/fvl5HY4wJI1VdCOA3jR0iUgfoDJzh7poAvCgiyaoa1jvnixpTSqOQa4SqWifkERlvrV6NdujAjsRERq1c6XU0xhhvNAa2qWou/LEyxHZ3f1qh7wxSUT2lru7jtUASMBanuuIanFUHTVnTogUfdOpEjZ49Of74472OxhgTvDgRWea3PVZVx5a0MbcM/DFVXVrAa1Vx8sOhkp4j0JtnP1fV3kXtCye7ebZ0LF++nHPPPZe1a9dSrVo1r8MxxgQp0JtnRWQjcK6qrnYv360DktxeUizOzD4pOLODP46zuN9S4HegItAKZxXxV4BnVTWrJPEGOqbUQERqq+ouN/jaQP2SnNBEqNxc9KqreHP1ah566CFLSMaUY6q6U0RWAoOBd93HFaqahnP57mwRaQz0BhoBh4D/AQtLmozyBNpTuglnNoeP3F1n43TfXg7m5MVhPaUw+/VXDvXowT3x8YzatMkq7owpI4rqKYnI88BFQD1gF5Cuqm3diRPGATWBPcAwVV0b9ngDXU9JRE4A+uCMKc1X1e/DGNdfWFIKL5/PR4+OHfnXgw9y4UUXeR2OMSZEwrmeUjgE/Oewm4SKTEQikoqTXZNwrkEOU9Wfj3FsS2AF8JKq3hFoLCbE1q5l4jffEFOpEhdceKHX0RhjyrGiSsJLUmWRd8PVuyIyFOeGq7/c7OIOnI3BuXHLeGX/fvSUU+DIEZ6YMeOoexWMMaa0FXr5TkQ6UIwqi8IqNtwBMv+278VZQrcKUKWonpJdvgsTn49Prr2WmevW8eKiRV5HY4wJsTAvh54KtFbV6SJSBUhQ1d3BtFloT0lVV1G8KouAbrgSkfZAf+A0wFaN89CBzEyu/uQTPv74Y69DMcZEERG5EqcALgGYDjTEmY4oqDVuAhpTUtUtOGWBQROReOBV4Go3aRV27HBgOEBCQkIoTm/8/fe/zPv8c/r27UvHjh29jsYYE11ux5lg4UsAVV0rIvWCbTTUdb9bgIYiEut3+a6Buz9PfaAFMMtNSDUAEZFqqnrUdNTuWNVYcC7fhTjW8k2V7KlTyfnqKx768UevozHGRJ8sVT2Qr2NR4Iq0xRHSpFTEDVd5x2wGaudti8iDBDCmZEJMhEd692ZnkyZc1KKF19EYY6JPujumpABuYVvQ08+F4w7JvwPjROR+3BuuAERkFnC/qi4r7M2mFOzdy56DBxn90kt8/fXXXkdjjIlOtwPjgZbuFEWZwHnBNlqcm2dDXmVRHFZ9F0L/+AcZr77KXRdeyMvjxnkdjTEmjMJcfReLU50twNq8IrdgBLQuhltlMQN41t3VEJgU7MmNNzJOP51Rubn884EHvA7FGBOlROR0nKGXNar6I1BVRIJegC3QxZryqiz2gVNlgTNPkolCTy1ZwvbLLqN58+Zeh2KMiV6jgAy/7QzgqWAbDXRMKSxVFqaU7d5N5siRvDd2LJ8us6E9Y0xQRP3Gf1TV517OC0qgPaWwVFmYUjZ7NhVHjmRonz7WSzLGBGu/iJyYt+E+D3rgP9ClK1Jxqixa48zMkAmcp6rrgw0gUFboELx9+/bRu1kzpk6qvuwAABp3SURBVCxfbknJmHIiXIUOItITmAL8gFPo0Bq4SFW/CqbdQGd0WOdmwZBWWZhSlJXFyy+/TLtzzrGEZIwJmqouEZE2QE+cvLBYVfcE225AScmtsvhGVde42zVEpLOqzgs2AFMKDh7E16oVuzIyuNsmXTXGhIiq7hGRz3BziYhUUtXMYNoMdEwpLFUWppRkZvJTkybEtm9Pu3btvI7GGBNBRORcEVkhIitF5DsRCWiVTxG5SES24kzSvR844D4GJdDqu7BUWZjSkV2jBmdv28bEiRO9DsUYE0HEKal+BzhVVVe7KzgsEpFpquor4u0jgUuArwI4NmCBJqX9InKiqn4NoauyMKVg6VJmfvEFxx13HD169PA6GmNM5PEB1d3nNYDfAkwyu1V1caiDCTQp3QlME5GjqixCHYwJPb31VtqtXMm/pk/3OhRjjDfiRMT/xsSxeauFq6qKyCXAdBE5CFQFzgmw3akicgPwPnA4b2ewY0rFmfuuJiGusigOKwkvmU/efJPxI0cy7scfbalzY8qhwkrCRSQO+AR4QFUXicjJwASgjaoeKKJd/96U4uQGVdWghnYCLXTATUKfAfOBIyJSKZgTm/BTVR4cM4aBjzxiCckYU5COQANVXQTgPh7EuRpWKFWN8fuJzXsMNqBAJ2QNS5WFCaNvvyW9b1/id+zgggsu8DoaY0xk2go0EpGWACLSGmde01KbGCG/QHtKeVUW8aHMiCaMfvkF/fprht18M7Gx9l9ljPkrVd0B3AB8ICKrgInA1YEsSyQiHURkiYhkikhu3k+wMQU6zdBSVe0e7MmCYWNKxbN+/XpOOvFEft20icqVw7KUijEmCoRxmqFFwP3AM8BZwE3AflV9Mph2A+0pTRWRG0SklohUyvsJ5sQmjLZv57nnnuPa4cMtIRljwqWiqs4FYlT1N1W9DxgQbKOBloQ/6j6Oxq/KArDrQpEmIwNt2ZIGOTkMW+/ZZWFjTNmXt3zRbhHpgDM+1TTYRgOdkDXgKj3jsbg45vbpw0FVGjRo4HU0xpiy630RSQIeBxbidFKCXs464PuUvGZjSoHJysqiefPmfPTRR3To0MHrcIwxHgvXmFK+c8TjXM4Luio70JLwgKssRCTVPXad+5hSwDFXuxP/rRSR70Xk1mB/EQN8+SWL776b1i1bWkIyxoSViCzMe66q2aq6339fSQU6pvQScB/5qiyOcewrwGhVfdddoXYM0DffMZOBt9wpLqoCq0Vkgap+V+zfwPxBX3iBlOnT+X8ffOB1KMaYsu+oYjd3ku5awTYa6FhRQFUWIlIH6IwzTQXuY2cRSfY/TlUz/GYdrwTE4y61bkru87/9jWsbNeKscwKdusoYY4pHRP4pImlAOxHZmfcD7AO+DLb9QJPSUVUW7uBWQVUWjYFteavSuo/b3f1HEZHz3QleNwGjVPX7YkdvjvL0888z6O67iYmxuhRjTNiMBboBc9zHvJ/Gqnp9sI0Hevku5FUWqjoDmCEiTXBmIJ+lqmv9jxGR4cBwgISEhGBOV7Zt387hfv3w/f47QydN8joaY0wZpqr7gH0icjWwT1WzwCl2EJFkVU0Lpv2A/qRW1WdUNV1VP8G5ZlhXVUcVcOgWoGHeAoDuYwN3/7Ha3gwsBc4t4LWxqtpVVbvGxQWaP8uh339nR3o6ZwwZQmJiotfRGGPKh5kc3bFJAD4MttFAq+8CqrJQ1Z3ASmCwu2swsCJ/5hSRVn7PawOnAXb5roTSmzShU3Y2g++7z+tQjDHlRwX/tZNU9SBQMdhGAx18KE6Vxd+BW0RkHXCLu42IzBKRru4x14vIDyKyEpgLvKiqc4odvYFNm3hrzBjOP/986tat63U0xphyxL+IzS10C3pAu9BrYiLyT5xVZ6u71RV5KgHvFfQeVf0JOLGA/Wf7Pf9/JYrW/IUOGcJp33xD70WLvA7FGFO+PA8sEpG33e1hOHUHQSlqoGYs8D/gRZx7k/JklPbKs6ZgS844g6lpaYzq2rXog40xJkRU9Q0R+RU4G2c+1OtU9Ytg2y00KYW7ysIE74GFC7nq/vu9DsMYUw6p6gJgQSjbDLSkbSZOMUKWu51XZdEjlMGYYtizh50PP8zWVau4+OKLvY7GGFPOiMg3FDDpQbBr7wWalP5SZSEiQVdZmCDMmUOdZ5/llr/9jQoVKngdjTEmSrnf5c8CpwOHgSWqOjyAt97h97wiTrX19mDjCfjmH//LdaGqsjAllzFgAL2rVWPmA0HPFG+MKd9G4iSjVHc+0oDKeFX1c/9tEZmDM8tDUAJNSmGpsjAlN27cOFLPOouGDRt6HYoxJkqJSBWc7/NGefORqurvJWyuGtA82JgCXeQvLFUWpmR08GCYM4ebp03zOhRjTOSLE5FlfttjVXWs+7wFkA48ICKnAQeA+1S1yCUo8o0pxeAkpKeDDjbQA8NRZWFKwOdje3o6mpDAKaec4nU0xpjIl6Oqx7pnJA4nmaxQ1X+KyInAhyJyvKpmFNGu/5hSDrBBVUtnTClcVRamBGJiuKFiRc5/+GFExOtojDHRbRNOQpkAoKpfi8guIBVYVtgb848phUqgPaWwVFmYYjp0iM3LlrF48WImTpzodTTGmCinqrtEZD5wBjBHRFKBOsAvx3qPu5bSMde/U9U6wcQU6JhSWKosTDFNmECj667j7mHDqFSpUtHHG2NM0f4OvCEiTwPZwBWqureQ4/MuBV4LJOHM/CPANcDWYIORPxeALcabRKoD36pqi2ADCFTlypX14MGDpXW6iHRo3Toe79SJq777juYtSu2f3hgTxUQkU1Urh6Hdz1W1d1H7iqskY0ohq7IwxfPeF1+wsl8/S0jGmEjQQERqq+ou+GMZovrBNlqSMaWQVVmYwOnLLzP/2We5+cUXvQ7FGGMA/gusFJGP3O2zgceCbbREl++8UK4v3x06RHb9+kwUYUh6OjExNpmGMSYw4bp857bdHuiNM6Y0X1WDXqy1qPWUwlplYQKUmMh1fftySrdulpCMMRFDVb8Dvgtlm4X2lESkqfu0wCoLVR0ZymAKU557Stu2beOEE05g48aNVKtWzetwjDFRJIyFDifhzJvXHKeDI4CGtSRcVTe5J++dr6LiVhH53A3IhNOUKWSNGMH1AwdaQjLGRJLXgYeBr4DcUDUaaKFDWKosTNGyDx1ix/btDPvHP7wOxRhj/B1S1fGhbjTQpBSWKgtTtEnAW71782m7dl6HYowx/maJyABV/TiUjQY0aq6qo3ES0WrgB+BsVX25oGNFJFVElojIOvcxpYBj/i0iP4jIKhFZLiL9g/klyqy1axn9wgvcfPPNXkdijDH5XQ98JCL7RGSniKSJyM5gGw15SbiIzAPeUNV3RWQocI2q9s13TH/gS1XNFJEOwOdAfVU9dKx2y12hQ0YGufXr81pcHNft3k1sbKzXERljolAYCx2aFrQ/rxahpAKd0SGgKgt3RdrOOJP7gTPz7Iv+q9a6Qc/2e9t3bntJhGDepDKjYkVe69KF+C5dLCEZYyKOqm4SkTigJc6tQ+tUNSfYdgMdUwq0yqIxsE1VcwFUNVdEtrv7047xnmHAelX9S0ISkeHAcICEhIQAQy0b0vbt4+7vv+eXqVO9DsUYY/5CRLoCk4EjOB2LOBEZpKrfBtNuoEkpLFUWItIbJ9mdUdDr7uqIY8G5fBfq80espUtZ+uSTXDJwIElJSV5HY4wxBXkOuFpV5wG4K9e+AJwcTKOBJqVAqyy2AA1FJNbtJcUCDdz9RxGRnsC7wEBVXVusqMs43zvv0HPaNBouWuR1KMYYcyyV8xISgKrOF5Ggx64CnbMmoCoLVd0JrMRZBBD3cYX/eBKAiHQD3gcuDrarVxZN7d2b6zp2pGOPHl6HYowxx5Lp9o6AP658ZQbbaEDVd8WpshCRVsA4oCawBximqmtFZBZwv6ouc5fCaAZs83vrFYVN5lduqu9U6d2nDzfddBOXXHKJ19EYY6JcoNV3IvIA8CBwgqquDuD4bsAHOGNKClQABqnq8qDiDbQkPBxVFsVRLpLSoUMc6tCB29LTGb1jB/Hx8V5HZIyJcoEkJRHpDDwKtAHOCSQpue+Lx8kLAvykqtnBxhvQ5Tu3ymI9MBWYDvzs/hImlNLS+OXgQXpddJElJGNMqRCRCsBo4EYKWRWigPedDlRS1dXuVa7KItK3qPcVJdAxpbwqi1RVTcGZJfyFYE9ujrazYkV6ZWZy1uOPex2KMab8eAh4V1U3FPN9o4AMv+0M4Klggwm0+i4sVRbGz+bNjBszhkGDBlG7dm2vozHGlB1xIrLMb3use7tNXhV0N+DuErQr6jf+o6o+t+I6KIEmpUwROU1V50PoqizMn3wjRjBk6lTSlwc1RmiMMfnlqGrXY7zWG2gFbBARgEbAbBG5WlXnFNHufhE5UVW/BhCRE4GgB/4DTUq3AR+IyFFVFsGe3PxpdpcuLP3pJx7o2NHrUIwx5YSqPgE8kbctIhuBcwMsdLgTmCYiP7jbbYCLgo0poKSkqt+IyPGEuMrC/Ok/06dz9yOPeB2GMcYERFWXiEgboCdOXlisqnuCbTfQCVlPB77Jy54iUkNEOvuPM5kSOnyYbX//O3Hbt3Peeed5HY0xphxT1WbFPH4PMCuUMQRafReWKgsDfP01dd9+m9vPPttmAzfGlHuBJqW/VFkA9g0aAttTUjihWjVOtzJwY4wJOCntdysrgNBVWZR7qrzyyiv0HTKEGjVreh2NMcZ4LtDqu7BUWZR3uUOHcvyUKVy+YoXXoRhjTEQItPouLFUW5ZoqP+zdizRsSKtWrbyOxhhjIkKgPaWwVFmUZz5VBm/cyHOvvOJ1KMYYEzECTkomhPbtY+Fbb1GhQgX69evndTTGGBMxAi10MKE0diy9br+dh664AndqD2OMMVhS8sTSDh24tXZtzrrlFq9DMcaYiGJJyQNPvPIKqQ88QFycXT01xhh/9q1YmlTZe9llZM+fz9XvvON1NMYYE3Gsp1SatmxBZ87kmt69qVzZlqMyxpj8Qp6URCRVRJaIyDr3MaWAY84UkWUickREys0cejsSEmiVkMCpL7/sdSjGGBORwtFTegUYraqpOOu+jyngmF+Bv+FM9Fo+7N7Ns888wyVDh1K7fn2vozHGmIgU0jElEakDdAbOcHdNAF4UkWRVTcs7TlV/cY8fGMrzR7IjgwbRf+FCUjds8DoUY4yJWKHuKTUGtqlqLoD7uN3dX36pMjkxkW29e9OoUSOvozHGmIgV0dV3IjIcGA6QkJDgcTQlt+P337n5q69YvTqQFYaNMaZ0iEgS8A7QAjgC/AJc739lq7SFuqe0BWgoIrEA7mMDd3+xqepYVe2qql2j9p6eFSuYN2wYVw4ZQoMGDbyOxhhj/CkwUlVbqmp7YD3whJcBhTQpqepOYCUw2N01GFjhZdb12sGXXuKszz7jrltv9ToUY4w5iqruVtUFfru+App6FA4Qnuq7vwO3iMg64BZ3GxGZJSJd3eeniMhWYARwvYhsFZH+YYjFc/ckJvLKVVdRL+UvlfHGGBMxRCQGuAGY4WkcfqucR7TKlSvrwYPRtdjttl9/5YSuXVmzZg1169b1OhxjTDkkIlnA9367xqrq2AKOGw00BC5SVV9pxZdflA7URIGFC6ncvz//vugiS0jGGC/lqGrXwg5wJzFIAc7zMiGBJaWw+fX33/kxJ4dhjz7qdSjGGHNMIvIo0AU4R1WPeB2PzX0XJrePG8faxx4jqUkTr0MxxpgCiUhb4B6cKunFIrJSRKZ6GZP1lEItK4tfbrmF9d99x//+9z+vozHGmGNS1R+AiFpp1HpKIZb78cccP3YsL19+ORUqVPA6HGOMiSrWUwqx0Zs2sbJ7d163sSRjjCk2S0oh9PuaNTz88MN88cUXiERUj9gYY6KCXb4LlblzqXHCCTzarx+tW7f2OhpjjIlK1lMKkdmbN7OzShWGjh7tdSjGGBO1rKcUArt37+aa++6j6YwZVEpK8jocY4yJWpaUgqTz5vFLx44MPfdcevXq5XU4xhgT1SwpBWnRO+9QfedO7n/4Ya9DMcaYqGdJKQg//vgjF86cSc5XX1G5Th2vwzHGmKhnSamEDj32GI8OGMCTTz5J244dvQ7HGGPKBFu6ogSy9+3jtwYN2NasGT1/+MHrcIwx5phEJFNVK3sdR6CsJLyYfD4ffx8xggOnnML4KVO8DscYY8oUS0rF4Fu4kE9vuYWfEhP5ZPZsYitHzR8fxhgTFezyXYBycnL4pHt32q9ZQ821a6lqS1IYY6JAtF2+s0KHoqiye/16zjrrLEbXqUPNdessIRljTJhYUirC5osvZn+rVpzUqhUzP/qIqo0bex2SMcaUWSFPSiKSKiJLRGSd+5hSwDGxIjJaRNaLyC8icl2o4wjKwYNs/uknrr76av7fwoVkDhvGQ889R2xsrNeRGWNMSAXynV2awtFTegUYraqpwGhgTAHHDAGOB1KAnsCDItIsDLEUS25uLl9++CEZSUm82bkz9erV482ff6b166+DJSRjTNkUyHd2qQlpoYOI1AHWAUmqmisisUA6kKKqaX7HfQS8qaofuNsvAptUddSx2i5poUN2djb7P/oIPXSII7174/P5iJ87l8P79rGlc2e2bt1Kmyef5Lc9e7giM5NGjRoxukYN2vzrX1Q/44xin88YYyJJYYUOgX5nl6ZQl4Q3Brapai6A+0tud/f7/4JNgE1+25vdY44iIsOB4QAJCQklCmj16tWkX3optX0+zk5OJiYmhom7dlE5JoY7O3akcePGNIyLo0HPnix78kmaWBGDMaZsiRORZX7bY1V1rPs80O/sUhPR9ym5/3BjwekplaSNTp06wc8/Q24u2487ztn5+++QkMDimjVDFqsxxkSoHFXt6nUQgQp1UtoCNBSRWL+uYAN3v7/NQFPgG3c7f88ptPL3furWDdupjDEmigT6nV1qQlrooKo7gZXAYHfXYGBFAdcm/wf8TURiRCQZuACYHMpYjDHGFK4Y39mlJuQzOohIK2AcUBPYAwxT1bUiMgu4X1WXudn4ReBM921P+l3jLJDXMzoYY0w0KmpGh2N9Z5dWfH+Jx6YZMsaYssumGTLGGGNKyJKSMcaYiGFJyRhjTMSwpGSMMSZiRE2hg4j4gEMlfHsckBPCcMLBYgwNizE0Ij3GSI8PIifGRFWNmg5I1CSlYIjIski/o9liDA2LMTQiPcZIjw+iI8ZIFDXZ0xhjTNlnSckYY0zEKC9JqdDZIiKExRgaFmNoRHqMkR4fREeMEadcjCkZY4yJDuWlp2SMMSYKlPmkFGnrz/sTkSQRmSUia0XkOxGZ4s6aHpFE5AERURFp53Us+YlIRRF5WUR+FpHvRSTiLp2IyLkiskJEVrr/3xdFQExPiciG/P+vkfS5KSjGSPrsHOvf0O/1iP3cRKIyn5SIsPXn81FgpKq2VNX2wHrgCY9jKpCIdAZ64KyFFYlGAoeBVFU9Afi3x/EcRUQEeAe4QlU7AkOBcSLi9WdwGtCLv65nFkmfm4JijKTPzrH+DaPhcxNxvP5AhJW7/nxnYIK7awLQOVJ6I6q6W1UX+O36Cmfxw4giIhVwvphuxPkyiCgiUgUYBvxb3UFSVf3d26gK5AOqu89rAL+pqs/DeFDVhap61IJukfa5KSjGSPrsFBQfRP7nJlKV6aREAevPA3nrz0cU9y/mG4AZXsdSgIeAd1V1g9eBHEMLIB14QESWicgCETnF66D8ucnyEmC6iGzC+ev6Sm+jOqao+dxARH92Iv1zE5HKelKKJi8AB3AWP4wYItIT6Aa85HUshYgDmuOsmNkVuAuYIiLVvA3rTyISB/wLGKiqTYHzgPfdXp4JTsR9dqLkcxORynpS+mP9eYBIWH++ICLyFJACXOr15ZwC9AZaARtEZCPQCJgtImcW+q7StQlnjrEJAKr6NbALSPUyqHw6Ag1UdRGA+3gQaO1pVAWLis8NRPRnJxo+NxGpTCelSFx/Pj8ReRToAlygqke8jic/VX1CVRuoajNVbQZsBfqr6hyPQ/uDqu4C5gNngFM5BtQBfvEyrny2Ao1EpCWAiLQG6uEM0EeUaPjcQGR/dqLhcxOpyvzNs5G2/rw/EWkLrAbW8ecM6BtU9ULvoiqc+1ffuaq62utY/IlIc+ANIAnIBu5V1Y+9jepoIjIEuBun4AHgAVWd5mFIiMjzwEU4CXIXkK6qbSPpc1NQjDjjcxHx2TnWv2G+YzYSgZ+bSFTmk5IxxpjoUaYv3xljjIkulpSMMcZEDEtKpkwTkQdFJMFv+yERubQUzvuZO2ZojCkGG1MyZZqIKFBVVQ+U4jlrAEtUNRLLvY2JaNZTMmWWiIx2ny52J0GtISJvicjN7usPishEd2LPX0TkfRHpJCLzRGS9iIzya6u+iHwgIkvdCV/vKeTU5wCzCoinmYjsEpHH3YlZfxKRLiLyqjup6NciUs899iQR+daN+wcRGfyXsxhTBllSMmWWqt7kPj1JVTuq6t4CDuuCcx9OS5ybHZ8ABgDtgSv9Zsd+G3heVbu77xkgImcc49QXANOP8VoSsFBVOwGvA3NxJj5tDywHbnaPuwt41p28tR0QUeXtxoRLnNcBGOOx2aq6D0BEvgNWuTdiHhGRtUALEdkO9AGSncm+AaiKMxvDp/6NuZNwdgUWHeN8B1T1I/f5t8BWVV3pbi/HvQEY52bgf4lIU+BTd5YKY8o8S0qmvDvs9zy3gO04nCsKCnRT1ewi2usHfJ43mWkB/GceONb5UNX/isiHwOnACyIyR1XvK+qXMSba2eU7U9bt58/lIkpEVfcDX+LMxgCAiDTOG//JZyDHvnQXMBFJVdX1qjoGeA7oHmybxkQD6ymZsu5pYJ6IHMK5BFdSQ4BnReR7d3s/cA2wI+8AdyG/M4ARQZwnz60ichqQhdO7uiUEbRoT8awk3JgQEZEeOHPuned1LMZEK0tKxhhjIoaNKRljjIkYlpSMMcZEDEtKxhhjIoYlJWOMMRHDkpIxxpiIYUnJGGNMxLCkZIwxJmL8f7nqv8wqcGJMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(times, rel, 'k', linewidth=1, label='MSM')\n",
    "tax = ax.twinx()\n",
    "tax.plot(time, g_k*(ode_solution**4), 'r:', label='exp')\n",
    "ax.set_ylabel('conductance (dMSM)')\n",
    "tax.set_ylabel('conductance (Exp)')\n",
    "ax.set_xlabel('time / ms')\n",
    "fig.legend()\n",
    "\n",
    "# this plot does not include scaling of relaxation observable to conductance as above"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}