tautomer_energies.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c4d56e8c-a6b5-4d58-8933-a927dc1cace2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "16273999-6477-4663-a427-93ca622ccd3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "temp = 298\n",
    "RT = 0.59  #kcal/mol\n",
    "x = np.linspace(-3, 3, 61)\n",
    "p = np.vstack([np.zeros_like(x), x]).T\n",
    "\n",
    "weights = np.exp(-1*p/(RT))\n",
    "weights = weights/weights.sum(1)[:,None]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cceb2b38-7838-4c59-a96e-dca26b854330",
   "metadata": {},
   "source": [
    "# plot \"error bars\" of 1kcal/mol\n",
    "This is the range of values you would observe if you had accuracy of +/- 1kcal/mol. \n",
    "\n",
    "Explanation: The blue line is the **true** occupancy of the minor tautomer with relative energy in kcal/mol compared to the major tautomer. The occupancy of the major tautomer is 1-occupancy(minor).\n",
    "\n",
    "If we had error bars of +/- 1 kcal/mol, then you would be fooled into thinking the minor tautomer has a population of within the gray bars. \n",
    "\n",
    "As an example, let's say the 'major' tautomer is equivalent in energy to the 'minor' tautomer. This means the difference is 0 kcal/mol. At 0kcal/mol, the blue line is at 0.5 occupancy for the minor tautomer. But, if we had error of 1kcal/mol either side, you would measure an occupancy for the minor tautomer of between 15% and 85%. \n",
    "\n",
    "Let's say we have a minor tautomer with a true energy difference of +1kcal/mol. It's true population would be 15%. With an error of +/- 1kcal/mol, you would predict a population of between ~3% to ~50%\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1c6208bd-6a55-4a98-ad84-5a2e6b2ddd62",
   "metadata": {},
   "outputs": [],
   "source": [
    "p = np.vstack([np.zeros_like(x), x+1]).T\n",
    "weights_up = np.exp(-1*p/(RT))\n",
    "weights_up = weights_up/weights_up.sum(1)[:,None]\n",
    "p = np.vstack([np.zeros_like(x), x-1]).T\n",
    "weights_down = np.exp(-1*p/(RT))\n",
    "weights_down = weights_down/weights_down.sum(1)[:,None]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4363f487-1f11-4998-9a5d-f693e0038b8c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "low: 0.8448694704908302 high: 0.15513052950916986\n"
     ]
    }
   ],
   "source": [
    "print('low:', weights_down[30][1], 'high:',weights_up[30][1]) #equienergetic tautomers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b5e18a21-f4e9-4676-ae37-13b914dbe763",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "true: 0.15513052950916986 low: 0.5 high: 0.03261480252335784\n"
     ]
    }
   ],
   "source": [
    "print('true:', weights[40][1], 'low:', weights_down[40][1], 'high:',weights_up[40][1])  #1.0 kcal/mol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "7cf93a34-8114-4b83-841a-a62abd62d16e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAACQ1ElEQVR4nOzdeXxU5dUH8N+9d/Y9k0km+wYIuAAarIKiooJiXd9qaWnFBfuK1AVTtVLUght9W0VsKygqUtfiWmuLSqqyKCKCQVyQLQnZZjLJ7Etmv+8fcabEJJMZmHAnM+f7+eQT5ubOnZOHzMyZZzkPw/M8D0IIIYSQLMEKHQAhhBBCSDpRckMIIYSQrELJDSGEEEKyCiU3hBBCCMkqlNwQQgghJKtQckMIIYSQrELJDSGEEEKyikjoAI61aDSKjo4OqNVqMAwjdDiEEEIISQLP83C73SgpKQHLJu6bybnkpqOjA+Xl5UKHQQghhJAj0NrairKysoTn5Fxyo1arAfQ2jkajSeu1Q6EQNmzYgJkzZ0IsFqf12tmG2ip51FbJo7ZKDbVX8qitkjdcbeVyuVBeXh5/H08k55Kb2FCURqMZluRGoVBAo9HQH/8QqK2SR22VPGqr1FB7JY/aKnnD3VbJTCmhCcWEEEIIySqU3BBCCCEkq1ByQwghhJCsQskNIYQQQrIKJTeEEEIIySqU3BBCCCEkq1ByQwghhJCsQskNIYQQQrIKJTeEEEIIySqU3BBCCCEkqwia3GzevBmXXHIJSkpKwDAM/vGPfwx5n02bNqG2thYymQw1NTV48sknhz9QQgghhIwYgiY3Xq8XEydOxF//+tekzm9qasJFF12EadOmoaGhAb/73e9w66234o033hjmSAkhhBAyUgi6ceasWbMwa9aspM9/8sknUVFRgRUrVgAAxo8fjx07duCRRx7BT37yk2GKkhBCCCEjyYjaFfzTTz/FzJkz+xy74IIL8OyzzyIUCg24+2ggEEAgEIjfdrlcAHp3LQ2FQmmNLxgMxh8jFsvhu5fG/j3QjqY//NlQ30e6WNun+/8gG1FbJY/aKjXUXsmjtkrecLVVKtcbUcmN2WyG0Wjsc8xoNCIcDqO7uxvFxcX97rNs2TIsXbq03/ENGzZAoVCkLbZwFNhkYiDlGOx46xNIOUDGAVKOP+zfgJQFsiQ/SYv6+nqhQxgxqK2SR22VGmqv5FFbJS/dbeXz+ZI+d0QlN0D/ngue5wc8HrNo0SLU1dXFb7tcLpSXl2PmzJnQaDRpi6vLHcBv/rhpyPNYBpCLGCglLBRiFkpx7N8MlGIWagkLlaT3mFrCQCXpPUctYaEQAyzDxH9n4L+/f6Jjh38//LyYwdru8HMZhunTe8SybJ/jh/+cZVkwDAOO4+LfRSIRWJYFx3HgOA7RaBS7du1CbW0txGIxWJYd8Cv2WLF/56JQKIT6+nrMmDFjwN5J8l/UVqmh9koetVXyhqutYiMvyRhRyU1RURHMZnOfYxaLBSKRCPn5+QPeRyqVQiqV9jsuFovT2uhSSRSXnFSIQ+1msBIF/BEevhCPnlAUPaEofKEoojwQ5QFviIc3FAEQSekxWAbQSDloZRx0st7vsX/r5Bz0cg75ChH0ChHy5BxEbPqSgVhidHiC9MNjh39Fo1GEw+GEPweA5ubmeKIEoE9CE0uYYslNLEk6/CuWLMW+Dk+gDj+WDdL9N5vNqK1SQ+2VPGqr5KW7rVK51ohKbqZMmYJ33nmnz7ENGzZg8uTJgv+x5aukeOiSsfjoow7U1BRCIpH0+TnP8whEeHiDUXiDvcmOJxiJ3/YGo/AEo3AHInAHonAHv//+/e1YcuTwR+DwR3BoiHgYABrZfxOeAqUIhcre7wWq778rRJCIknvjP7xnJh0ikQisViv0ej04jgOAfsnP4d8jkQhCoVC/n8ViiyVch/fyxBIbkUgEsVgMiUQCsVjcJzH6YbKUq71DhBCSTQRNbjweDw4cOBC/3dTUhF27dkGv16OiogKLFi1Ce3s7nn/+eQDA/Pnz8de//hV1dXX41a9+hU8//RTPPvssXnnlFaF+haQxDAOZiIFMxCL/CKb6hKM8nP4InN8nN/F/9/R+t/VEYOsJw+aLwN4TRoRH/Jwme3DQ6+bJORQqRSjWiFGiFqP4+68SjRh6OXdM3+wPT6BiCU+qotFov69AIACfzxe/fbhYT0/su1gsjvf2SSSSeGJ0+HdKgAghJLMJmtzs2LED06dPj9+OzY255pprsHbtWphMJrS0tMR/Xl1djfXr1+P222/HE088gZKSEvz5z3/OiWXgIpZBvkKEfMXQ/2VRnofLH4G1JwKbrzfh6fKG0eUNw+INwfL9vwNhHvaeCOw9EeztDvS7jpRjUKQWo1wrRoVOgnKtBBU6Ccq0YsiS7PE51mI9N8mKRCKIRqOIRCKIRCLo6emBx+NBJNI7ZHh4shXr6ZFKpZDL5ZBKpX16hMRi8REnZYQQQtJH0OTmnHPOGXCCa8zatWv7HTv77LPxxRdfDGNUIx/LMNDJRdDJRRil7z/fCOgdAnIHorB4wzC7QzC5Q+j4/rvJ1ZsABSI8DjmCOOQIAoe88fsyAApVIlRoJajKk6BGL8WofCnKNGJwaZzncyzEkpZEw5qxHp9wOIxIJAK32w273d5nKCyW+EgkEiiVSshksnjvT6wHiHp8CCHk2BhRc25I+jAMA42Mg0bGYXR+/wQoHOVh8YTR7gqi1RlCqyOIFmcQLY4gXIEoOj1hdHrC+Lz9v0vzJByD6jwJRumlqMmXYpS+99/SDO3lSVasN0gkGvjpcnji09PTA5fLFU98YomTRCKBQqGAQqGID3tJpVLq6SGEkGFAyQ0ZkIhlUKLpnXtzalnfnzn8kd5kxxFEkz2Ig7YAGm0B+MM89nYH+gxxiVigRi/F8QUyjC/s/SpUZtefHcuy/SaQA4hPhA6Hw/D7/fB4PPE5P7H5O0qlsk9Pj1wup4SHEEKOUna9y5BjQifjoCuS46QiefxYlOdhcoVwwBbAQVsQjdYA9lsDcPgj2NcdwL7uAP6xxwkA0Ms5jCuQooBhwOoDGF0gB5uFQzYMw8RXYclksvhxnucRDocRCoXgcDjQ3d0NoLeXJzaspVarIZfLIZfLaTiLEEJSRMnNMPD5fGAYRvDl6ccSyzAo1UpQqpXg7OreYzzPo9MTxp4uP/ZY/NjT5cdBawC2ngi2tvgAcHj7UAfUUhYTi+SYVKLApGI5yjTirH5Dj/1t/PDvIxwOIxgMwm63o7u7O35erFfIarVCo9FQwkMIIUOg5CaNYm9WUqkUPp8PoVAILMvG51cMNmcjWzFM72qrIrUY02vUAIBAOIr91gC+MfuwrdGKJg8HdyCKjw958fH3k5YNCg4TixU4tUyBU0sVUElzY5gm1ssT2xYkGo0iFArB7/cDAPbv3w+ZTAaZTAatVguVShWfw0PJDiGE/FduvdsOs1jyMn78eITDYfh8Png8HjidzvjyYo7j4vMrcvENSSpicaJRjvEGCSZKu1A9qgYH7SHsMvVgV0cPvu3yo9sXwQcH3fjgoBscA5xolOP0CiWmVChRrM6h3rDvE+PY35Ver4/X7WlrawMASCQSyOVyaLVaqNVqKJXKnEuiCSHkh+hVcBiwLBufKFpQUBBfRePz+eBwOOB2u+H1esEwDKRSKWQyWc5OIhWxDI4vlOP4QjnmTOzt2dlj8WNnhw+ftfpwyBHEl+YefGnuwVPbu1Gpk+D0cgWmVqow1pBbCSLDMPGl5UDvsF8wGITP54PT6ewtFCmTQafTQavVQqlUDrj1CCGEZDtKbo4BjuOgUqmgUqlQWFiIQCAAr9cLp9MJh8MBl8uFaDQa/xSey5+8pSK2d+5NiQLzJgMdrhC2tXqxrcWLrzp74nV31n3lQJFKhOk1apw7So0KXf/VStkulhzHEphIJAK/3w+z2QyTyQSpVAq1Wg2dTge1Wg2ZTJZTySAhJHfl7ruogGJvSHq9HuFwGF6vN14Yzu12IxqN0rLg75VoxPifE3T4nxN0cAci2NHuw6ctXnzW6oXZE8Yru+14Zbcdo/QSTK9RY3qNGoYsW2qeLI7j4j2GPM8jEAjEJydLJBKo1Wrk5+dDrVZTjw4hJKvl5rtABhGJRNBqtdBqtSgpKYnP0bHZbHA6neB5HnK5HDKZLGt2tz5SaikXT2D8oSg+bfXio0Y3drT5cNAWxEGbFc/usOKkIjkuGKPGWVWqpDcGzTaxISqZTBYfvootO48NXeXl5UGtVufUqj5CSG6g5CaDsCwLjUYDjUaDkpISuN1uOBwO2O122O32+FyegQrG5RqZmI0nOk5/BJubPfio0Y1vOv3Ybe7B7u/n6FwwRoOLxmpRosndN/DDh694noff70dXVxc6Ozshl8uh1+vjiQ4NWxFCsgElNxmK4zjodDrodDqUlpbC5XLBarXC6XTC5XLFC7zlem8OAGhlHC4Zp8Ul47Qwu0P44KAb7+5zocsbxmtfO/D61w7UlipwyTgtTi1TjLj9r9KJYZj43040GkVPTw/a29thNpuh1WphMBig1WqpN4cQMqJRcjMCiMVi5OfnQ6/Xw+fzwW63o6urC3a7PT7Pgt6MehWpxfjFJD1+NiEPn7V68a+9Luxs92HH91+FShF+PE6Di8dqc6Z+zmAOX9UXCoXiw6EKhQIGgwF5eXlQKBTUm0MIGXEouRlBGIaJvxkZjUY4nU50dXXB5XIhEonE9ygiAMcymFqpwtRKFdpdQazf68L7+12weMN4bqcN63bbcck4La44QYc8OT0NxGIxdDpdvDfn0KFDMJlM0Ol08d4c6iUkhIwU9Ko+QonFYhgMBuTn58PtdsNqtcJqtcLj8cSTHPrE3atUI8GvTjVg7sl6bG7y4LWvHfHl5G9968Ss4zS48kQdClXU+3V4b04wGER3dze6u7uh0+lQWFgInU6X8yv4CCGZj5KbEY5hmPgkZKPRiO7ubnR1dcFqtUIul9OwwmGkIhYzxmhw3mg1Pmv14ZUvbdjbHcDbe5z413dOnDdajdkn5aFMSxO2gd7qx7FyBS6XC3a7HVqtFkajEXl5eZTkEEIyFiU3WUShUKCiogIFBQWwWq2wWCywWq2QyWRQKBQ0rPA9lmEwpUKJ08sV2GXqwStf2vGluQcb9rtRv9+Nc0epce0peurJ+Z5IJIJOp0MkEoHH48G+ffugVqvjSQ7N9yKEZBpKbrKQXC5HWVlZnyTHZrNRT84PMAyDk0sUOLlEgW8tPVi3245trT58cNCNzc0e/M/xWsyekAelhHoogN4VfFqtFpFIBF6vF/v374darUZxcTH0ej315BBCMgZ9lM9iUqkUJSUlOP7441FVVQWe52G1WuO7TJP/Or5QjqXnl+Avl5RhQpEMoQiPdV85cN3rh/D2tw6Eo7zQIWYMjuOg0Wig1+sRCASwf/9+7N27F3a7HTxP7UQIER713OQAiUSC0tJS6PV6dHZ2oqurC16vFxqNhoYUfuA4gwx/vLAUn7X68MyObrQ6Q1j5WTfe3uPEvMn5mFqhpJ6v78WKTkYiEbjdbrhcLuj1ehQXF0OtVgsdHiEkh1Fyk0PkcjmqqqqQn58Ps9kMq9UKANBoNDSkcBiGYXB6hRKnlinw7j4XXmiwod0Vwv0fmnGCUYZbpxSgKo/2ZoqJFZwMhULxQpMGgwFFRUWQy+VCh0cIyUE0LJWD1Go1Ro8ejbFjx0Kj0cDhcMDtdtOQwg9wLIOLx2nx3JWVmDMxD1KOwTedfix4uxVrd1oRDEeFDjGjiMVi6PV6SKVSdHR04Ntvv0VHRwfC4bDQoRFCcgwlNzmKYRjk5eVh7NixGD16NEQiEaxWK4LBoNChZRyFmMU1p+Tj2Z9UYmqFEhEeeGW3HTe93Yrd5h6hw8s4UqkU+fn5YBgGTU1N2LdvH5xOp9BhEUJyiODJzcqVK1FdXQ2ZTIba2lps2bIl4flPPPEExo8fD7lcjrFjx+L5558/RpFmJ47jUFhYiHHjxqG4uBherxdOpxPRKPVK/FCBUoTfn1eMe6cXQS/n0OYK4c532/HYxxa4AxGhw8sosWraeXl5cLlc2Lt3L1paWih5JoQcE4ImN+vWrcPChQuxePFiNDQ0YNq0aZg1axZaWloGPH/VqlVYtGgRlixZgm+++QZLly7Fr3/9a7zzzjvHOPLsI5PJUF1djeOOOw4ymQw2m41WVQ3izCoVnvmfCvx4rAYA8N5+F371Vgs2NtLQ3g9xHIe8vDxIpVK0trbiu+++g9VqpXYihAwrQScUL1++HPPmzcMNN9wAAFixYgXef/99rFq1CsuWLet3/gsvvIAbb7wRs2fPBgDU1NRg27Zt+L//+z9ccsklAz5GIBBAIBCI33a5XACAUCiEUCiU1t8ndr10X/dYUqvVGDVqFCwWCywWC7xeL9RqddonHEcikT7fRxoZB/z6tHycU6XEn7f1rqpatqkTm5rcuG2KAeo0bso50tsK6J2Po9Vq4fF4sHfvXhgMBhQXF6d9L7RseA4eS9ReyaO2St5wtVUq1xMsuQkGg9i5cyfuvvvuPsdnzpyJrVu3DnifQCDQ78VQLpdj+/btCIVCAy5rXrZsGZYuXdrv+IYNG6BQKI7iNxhcfX39sFxXKF1dXcN27cbGxmG79rEgAbBwHFDfzqC+ncXWFh++NR/C1WMiGK1J72ON9LY6XFdXF/bs2TNs18+25+Bwo/ZKHrVV8tLdVj6fL+lzBUtuuru7EYlEYDQa+xw3Go0wm80D3ueCCy7AM888g8svvxynnHIKdu7ciTVr1iAUCqG7uxvFxcX97rNo0SLU1dXFb7tcLpSXl2PmzJnQaNL77hMKhVBfX48ZM2ZkTf2YcDiMzs5OmEwmMAwDtVqdljovkUgEjY2NqKmpyYpl6OPHAhdZA/jDZgs63GE88a0IPztJh59P0IFjj669sq2tAIDneXg8HkQiERQVFaGoqCgtz5lsfA4OJ2qv5FFbJW+42io28pIMwevc/PCNkuf5Qd887733XpjNZpx++ungeR5GoxHXXnst/vjHPw76oi+VSiGV9q9JIhaLh+0PdDivfayJxWJUVlZCq9WitbUVDocDOp0OIlF6/nQ4jsuaN+yxhQqsvKwCT2zrQv0BN17e7cDuTj9+e5YxLftUZVNbAYBOp0MgEIDJZILf70dFRQVUKlVarp1Nz8FjgdoredRWyUt3W6VyLcEmFBsMBnAc16+XxmKx9OvNiZHL5VizZg18Ph+am5vR0tKCqqoqqNVqGAyGYxF2TootGz/uuONQWFgIp9OZUvdgLpGLWdwxzYjfnmWEQszg604/bnq7FZ8c8ggdWkaSSqXQ6/VwOp3Yu3cvzGYzrdQjhBw1wZIbiUSC2trafmNy9fX1mDp1asL7isVilJWVgeM4/P3vf8fFF19MO14fAzKZDKNGjUJ1dTUikQjsdju9EQ3i3FFqPHFpBY4zSOEJRnH/h2b85VMLghFaJfRDLMtCr9eDYRg0NjaiqampzyIAQghJlaDDUnV1dbj66qsxefJkTJkyBatXr0ZLSwvmz58PoHe+THt7e7yWzb59+7B9+3acdtppsNvtWL58Ob7++mv87W9/E/LXyCksy6K4uBhKpRKtra2w2WzQarXUTTuAEo0Yyy8qw9++sOK1rx3413cuNNqC+P25RdDJBR8RzjhKpRISiQRmsxlerzc+HEoIIakS9BV29uzZsFqtuP/++2EymXDiiSdi/fr1qKysBACYTKY+NW8ikQgeffRR7N27F2KxGNOnT8fWrVtRVVUl0G+QuzQaDcaMGYO2tjaYzWYoFAraR2gAYo7BDacaMLFYjmWbOvGtxY9b3mnDkvOLMUpP+1P9UGwLB6fTif3796OiogIFBQW0WSkhJCWCf3xcsGABFixYMODP1q5d2+f2+PHj0dDQcAyiIsmQSCSoqqqKF2gLh8O0G/QgTi1T4vGLy/D7/5jQ7gqh7t9t+O1ZRkytTM8E2mzCsizy8vLg8XjQ2NiIQCCAkpKSrJpMTQgZXjRRhRwVlmVRUlKC0aNHg2EY2O12qj47iHKtBI9fXIaTS+Twh3ks/dCMV760UXsNQqVSQS6Xo6WlBU1NTbR1AyEkaZTckKPGMAwMBgOOO+44KJVK2Gy2EV1NdzippRwemlGCy8b3ziVZ+4UNf9jciQDtMD4gmUwGnU6Hzs5OHDhwgFbpEUKSQskNSRu1Wo0xY8YgPz8fdrudPmkPgmMZLDi9ALdOKQDHABsbPbjj3XZYfWGhQ8tIIpEIer0eDocD+/fvh8PhEDokQkiGo+SGpFVsuXhJSQk8Hg96enqEDilj/XicFg9fUAK1lMW+7gBu/3cb2l2UEA4ktlzc7/dj//79MJvNNJxHCBkUJTck7UQiESorK1FZWQm/3w+v1yt0SBlrUrECf764HKUaMTo9YfxmfTuabFTjZSAMw0Cn04HjODQ3N6Ojo4MSHELIgCi5IcMiVg+nuroa4XAYbrdb6JAyVolGjEcvKkWNXgJ7TwR3vNuOPRa/0GFlLKVSCZlMhpaWFrS1tVEhSUJIP5TckGHDMAyMRiNqamoApLbpWa7Jk4vwxwtLMb5ABk8wirvfb8euDpo8Oxi5XB4vJNna2koT2AkhfVByQ4adwWBATU0NWJaFw+GgoYRBqKUcll1QEl8qfs9/TNhKe1INSiqVQq1Wo62tDYcOHUI4TBOyCSG9KLkhx4Rer8fo0aMhkUgowUlALmZx//klOKNCiVCExwMfmfFhIyU4g5FIJNBqtTCZTGhubkYoFBI6JEJIBqDkhhwzWq0Wo0ePhlwup+W8CUg4BounF2HGaDWiPPDIx13YYqbtBwYjFovjtXCam5uFDocQkgEouSHHlFqtxujRo6FUKgGAJoMOgmMZ1J1ZGC/293oTh3e+ozlLgxGJRMjLy4PVagUAqrFESI6j5IYcc0qlMj7JmIaoBscyDG46zYDZJ/UmOKu2W7FhPyU4g+E4DjqdDgBouwZCchwlN0QQsR3EVSoV7UeVAMMwmDspD2cX9/ZwPfaJBZubaA7OYGKba9rtdjQ1NdEcHEJyFCU3RFDV1dVQKpXUg5MAwzC4ojKKC0arEOWB/9tsxvZWKoyYiE6nQ3d3N5qammgVFSE5iJIbIiiFQoHq6mrIZDKaZJwAwwA3n27AOTUqhKPAAx+Z8aWJ6uAMJjZE1dXVhebmZqqDQ0iOoeSGCE6tVqOmpgYSiQROp1PocDIWxzK4c5oRp5crEIzw+P1/TPiuiyoZD0YkEsVXUR06dIgSHEJyCCU3JCNoNBrU1NSA4ziqZJyAiGWw+JwiTCqWoyfMY/GGDjTSXlSDEolE8To4ra2ttDqPkBxByQ3JGDqdDtXV1QBAe1ElIBGxWHJeMY4v7N2qYdH7HWh10sqgwYjFYmg0GnR0dKCtrY3mdhGSAyi5IRlFr9ejuroa0WiUdhNPQC5m8cD5xRitl8Lhj2Dxhg7Ye2ji7GAkEglUKhXa2trQ3t5OCQ4hWY6SG5JxDAYDqqqqEAwG4ffTnJLBqKQcHrqgBCVqMTo9YSz9wIRAmIZdBiOVSuObbXZ2dgodDiFkGFFyQzJSQUEBysrK4PV6qVZJAjoZhwdmFEMlYbGnK4BHP7YgSr0Sg5LJZJBKpWhtbYXdbhc6HELIMKHkhmQkhmFQUlKCoqIiOJ1OWumSQJlWgvvOLQLHAJuaPHihwSZ0SBlNqVSC53k0NzfD46GCiIRkI8GTm5UrV8brnNTW1mLLli0Jz3/ppZcwceJEKBQKFBcX47rrrovvJ0OyC8uyqKiogMFggMPhoJUuCUwsVuC2MwoBAC9/acd/DtCKs0Q0Gg38fj+am5sRCNBqM0KyjaDJzbp167Bw4UIsXrwYDQ0NmDZtGmbNmoWWlpYBz//4448xd+5czJs3D9988w1ee+01fP7557jhhhuOceTkWBGJRKisrIRGo6EqxkO4YIwGs0/SAQBWfGLBV+YeYQPKYAzDQKfTwel0orm5maoYE5JlBE1uli9fjnnz5uGGG27A+PHjsWLFCpSXl2PVqlUDnr9t2zZUVVXh1ltvRXV1Nc4880zceOON2LFjxzGOnBxLMpkMVVVVkMlkVANnCNfW5uPMSiVCUeD+D01od9ES8cGwLIu8vDx0d3dTDRxCsoxIqAcOBoPYuXMn7r777j7HZ86cia1btw54n6lTp2Lx4sVYv349Zs2aBYvFgtdffx0//vGPB32cQCDQp9s59uYYCoXSPlE1dj2aADu0VNtKKpWitLQUTU1NcLvdUCgUwxleRonNN0p23lHdGQZYPCHsswZxb70Jy2cVQy3lhjPEjJFqWwG9m7e2t7eD4zgUFRWBYZjhCi/j0GtW8qitkjdcbZXK9RheoH7+jo4OlJaW4pNPPsHUqVPjxx9++GH87W9/w969ewe83+uvv47rrrsOfr8f4XAYl156KV5//XWIxeIBz1+yZAmWLl3a7/jLL7+cU2+QJLe4gsCjX3FwBBmM0UQxf3wUIsFn2BFCyJHz+XyYM2cOnE4nNBpNwnMF67mJ+eGnJJ7nB/3k9O233+LWW2/FfffdhwsuuAAmkwl33nkn5s+fj2effXbA+yxatAh1dXXx2y6XC+Xl5Zg5c+aQjZOqUCiE+vp6zJgxY9Bki/Q60rbieR7t7e1oa2uDTqcDx2V/j0QkEkFjY2N8e4pkPVQSxB3vdmC/i8Vmlw43npo/jFFmhiNtKwDxopE1NTVpf23IVPSalTxqq+QNV1ulMi1BsOTGYDCA4ziYzeY+xy0WC4xG44D3WbZsGc444wzceeedAIAJEyZAqVRi2rRpePDBB1FcXNzvPlKpFFKptN9xsVg8bH+gw3ntbHMkbVVRUYFIJILOzk7o9XqwbG50SXAcl9Ib9miDHL89uwhLPjDh7T0unGCU4+xq9TBGmDlSbSugdwWV3W5He3s7lEolZDLZMEWXeeg1K3nUVslLd1ulci3B3hUkEglqa2tRX1/f53h9fX2fYarD+Xy+fm9ksRcwWkWTOziOQ3l5OTQaDe0iPoQpFUr89PsVVI99bEGLgyYYJ6LT6eB2u9HS0kK1lQgZwQT9yFtXV4dnnnkGa9aswZ49e3D77bejpaUF8+fPB9A7pDR37tz4+ZdccgnefPNNrFq1Co2Njfjkk09w66234kc/+hFKSkqE+jWIAKRSKSorKyESiWgPqiFce0o+JhT17iL+4Edm+EO0KmgwsSXiXV1dMJlM9KGJkBFK0ORm9uzZWLFiBe6//35MmjQJmzdvxvr161FZWQkAMJlMfWreXHvttVi+fDn++te/4sQTT8RVV12FsWPH4s033xTqVyAC0mg0KCsrg9/vRzBIPRKD4VgGi842Qi/ncMgRxONbLfSmnYBIJIqvoKItGggZmQSfULxgwQIsWLBgwJ+tXbu237FbbrkFt9xyyzBHRUYKo9GInp4edHR05NT8m1TpFSL87pwi3PVeOz5s9OB4oxyXjNMKHVbGkslkCAaDaGlpgUwmo5WVhIww9E5ARjSGYVBWVoa8vDyqYDyEk4rkuL62d8XUU591YW8X7bieiFqths/nQ0tLC1UwJmSEoeSGjHhisRiVlZWQSqW0EeIQrjxRhzMqeisYP/iRGS4/TZodTGz+jdVqRXt7OyXOhIwglNyQrKBUKlFRUYFwOEwbISbAMAx+M60QJWoxLN4w/ri5E1F60x4Ux3FQq9Xo6OhAV1eX0OEQQpJEyQ3JGvn5+SgpKYHb7aZlvAkoJRzuObcIEo7B5+0+vPIlTZpNRCqVQiKRoLW1FW63W+hwCCFJoOSGZA2GYVBaWgqDwUDzb4YwSi/FLVMKAAAv7rLhWwvtIJ6ISqVCKBRCS0sLrcwjZASg5IZkFY7jUFFRAblcTvNvhjBzjAbn1qgQ5YE/brbAR/VvEtJqtXA4HDT/hpARgJIbknXkcjnKy8sRCoXoU/YQfn16AQqVIpjcITz1WbfQ4WQ0lmWhVqvR2dkJm80mdDiEkAQouSFZKT8/H0VFRXC5XIhGqUdiMCophzvPMoIB8N5+Fz45RL1diUilUnAch9bWVvj9tJSekExFyQ3JSrH5NxqNJqWdZHPRhCI5rvp+/6kVn1hg9VFNl0Ri9W9aW1spcSYkQ1FyQ7KWRCJBeXk5GIahT9lDmHtyPkbpJXAFonj0Y9qeIRGGYaDVatHV1QWLxSJ0OISQAVByQ7KaTqdDSUkJPB4PLQ9PQMwxuPvs3uXhO9t9+Oce2m09EZFIBJlMhvb2dpq4TkgGouSGZL2ioiLo9Xo4nfSGnUiFToJfndq7PcMzO6w45KDJ2IkolUoEg0G0trbS9gyEZBhKbkjWE4lEKC8vh0QigdfrFTqcjHbJOC0mlyoQjPD4v01mBCM0PJWIVquFzWaDyWQSOhRCyGEouSE5QaVSobS0FH6/H6FQSOhwMhbDMPjNmYXQSlkctAXx/BdWoUPKaBzHQaVSwWQyweFwCB0OIeR7lNyQnFFYWIjCwkK4XC6aMJuAXiHCwjMKAQCvf+3AV2aqXpyITCYDz/NobW2lukqEZAhKbkjOYFkWZWVlUCgUtEfQEKZWqnDhGA14AI99YkEgTEueE4mVHKDqxYRkBkpuSE6RyWQoKytDOBym4akh/O+P8pGv4NDuCuHFXVSRN5HDqxfT8BQhwqPkhuSc/Px8FBQU0PDUEJQSDrdO+e/w1P5uqhWUiFQqBcMwaGtro+EpQgRGyQ3JObHqxXK5nFZPDeH0CiXOqe7dXPPRjy0IRykZTCQ2PEWrpwgRFiU3JCfJ5XKUlpYiEAhQjZIh3HR6ATRSFk32IF7dbRc6nIzGsixUKhU6OzuprhIhAqLkhuQsg8EAg8FAb0JD0Mk43HRaAQDg5S9tVNxvCDKZDNFoFG1tbZQ4EyIQSm5IzmJZFqWlpZBKpTQ8NYTpNSr8qEyBUBR47GMLIjQ8lZBWq4XD4YDZbBY6FEJyEiU3JKcplUqUlJTA7/fT3lMJMAyDW6YUQCFmsKfLT3tPDYFlWSiVSphMJtqVnhABCJ7crFy5EtXV1ZDJZKitrcWWLVsGPffaa68FwzD9vk444YRjGDHJNoWFhbT3VBIKVWLMm2wAADz3hRVmNy2lT0QulyMSiaC9vZ0SZ0KOMUGTm3Xr1mHhwoVYvHgxGhoaMG3aNMyaNQstLS0Dnv/444/DZDLFv1pbW6HX63HVVVcd48hJNuE4DmVlZRCLxejpoWq8iVw0VoMJRTIEwjxWbLXQUvohaDQa2Gw2dHZ2Ch0KITlF0ORm+fLlmDdvHm644QaMHz8eK1asQHl5OVatWjXg+VqtFkVFRfGvHTt2wG6347rrrjvGkZNso1KpUFxcDK/XS5+yE2AZBgunFkLCMWjo6MGG/VTpORGO46BQKGAymeDxeIQOh5CcIRLqgYPBIHbu3Im77767z/GZM2di69atSV3j2Wefxfnnn4/KyspBzwkEAggEAvHbsfHvUCiU9gq1setR5duhZWJb6fV62Gw2OBwO6HQ6ocOJiyVbmZJ0Fak4/HKiDmu+sGP15904tVQGrYwTOiwAmddWQG9xP5vNhpaWFowaNQosK/hsgLhMfB5mKmqr5A1XW6VyPcGSm+7ubkQiERiNxj7HjUZjUisMTCYT3n33Xbz88ssJz1u2bBmWLl3a7/iGDRugUChSCzpJ9fX1w3LdbJSpbdXV1SV0CP00NjYKHULciRKgVMGh3RfFio+aMWd0Zu09lUltFWO1WrF//36hwxhQpj4PMxG1VfLS3VY+ny/pcwVLbmIYhulzm+f5fscGsnbtWuh0Olx++eUJz1u0aBHq6urit10uF8rLyzFz5kxoNJojinkwoVAI9fX1mDFjBsRicVqvnW0yua1aW1vR3t4OvV6f1N/icItEImhsbERNTQ04LjN6SACgLs+P37xrwmddLK6sLcUJhTKhQ8rYtgIAj8cDkUiEsWPHQiKRCB0OgMx+HmYaaqvkDVdbpbLyULDkxmAwgOO4fr00FoulX2/OD/E8jzVr1uDqq68e8kVCKpVCKpX2Oy4Wi4ftD3Q4r51tMrGtSktL4XK50NPTA7VaLXQ4cRzHZdQb9olFSlw4RoP39ruw8jMrnri0HBwrfDIIZF5bAb2Ti61WK7q6ulBZWZkRiXNMJj4PMxW1VfLS3VapXEuwwV+JRILa2tp+3Vb19fWYOnVqwvtu2rQJBw4cwLx584YzRJKjpFIpSkpKEAwGqcLsEK6fnA/191szvE21bxJiGAZqtRoWi4Vq3xAyzASd2VZXV4dnnnkGa9aswZ49e3D77bejpaUF8+fPB9A7pDR37tx+93v22Wdx2mmn4cQTTzzWIZMcYTAYkJ+fT29CQ9DKOFxfmw8AeKHBCquPksFEpFIpotEo1b4hZJgJmtzMnj0bK1aswP33349JkyZh8+bNWL9+fXz1k8lk6lfzxul04o033qBeGzKsWJZFSUkJRCIR/H6/0OFktAuP02BcgRS+EI/V27uFDifjaTQaOByOjJy0Tki2EHxC8YIFC7BgwYIBf7Z27dp+x7RabUozpgk5Umq1GkajEa2trZBIJBm1hDeTsAyDW6YU4pZ3WrGxyYMLj/Ph5JLhWYmYDTiOg1Qqhclkgk6ng0wm/ERsQrINvVoTkkBRURHUajUVYBvC6HwpLh6nBQD8dVsXghGqXJyIUqmEz+dDR0cHVXkmZBhQckNIAhKJBCUlJQiHwzS5eAjXnKxHnpxDmzOEN762Cx1ORmMYBhqNBhaLhfY0I2QYUHJDyBD0ej3y8/PpTWgIKimHX53au7HmK1/aaWPNIcTKWLS3t1PiTEiaUXJDyBBik4slEgltrDmEc2tUvRtrRnis+owmFw9Fo9HA6XTCYrEIHQohWYWSG0KSoFKpUFRUBK/Xi2g0s7YayCQMw+DmKYXgGGBbqxeftXqFDimjcRwHmUwGk8lECyUISaOUV0vt3r17wOMMw0Amk6GiomLAisCEjHRGoxF2ux1utxtarVbocDJWpU6CK07Q4fWvHVi9vRunlCgg5jKnGm+mUSgUsFqtMJlMqKmpyajKxYSMVCknN5MmTUr45BOLxZg9ezaeeuopWuJIsopYLEZxcTH279+PcDgMkUjwSgoZa85EPf5zwI02Vwj/3OPAT07MEzqkjBWbXNzV1YX8/PyM2pGekJEq5WGpt956C2PGjMHq1auxa9cuNDQ0YPXq1Rg7dixefvllPPvss/jwww9xzz33DEe8hAhKr9dDr9dT5eIhKCUsrq3VAwBe+tIOh5+q8SYSm1zc0dFBlYsJSYOUP3o+9NBDePzxx3HBBRfEj02YMAFlZWW49957sX37diiVSvzmN7/BI488ktZgCREay7IoLi6G0+lEIBCgIdgEZo7W4F97XDhgC+BvX1hx29RCoUPKaLHKxVarFYWF1FaEHI2Ue26++uqr+PYIh6usrMRXX30FoHfoymQyHX10hGQgjUaDgoICuN1uKsCWAMcymH9a79Lw9/a50GgLCBxRZuM4DhKJBCaTCcFgUOhwCBnRUk5uxo0bhz/84Q99nnyhUAh/+MMfMG7cOAC9dRuMRmP6oiQkwxQXF0Mul9MKlyGcVCTHtCoVojzw5GfdlAwOQaVSwePx0NJwQo5SysNSTzzxBC699FKUlZVhwoQJYBgGu3fvRiQSwb/+9S8AQGNj46D7RRGSDWQyGYqLi9HY2Ai5XE77TiVww+R8bGv14ktzD7a2eHFGpUrokDIWwzBQKpUwm83Iy8uDUqkUOiRCRqSUk5upU6eiubkZL774Ivbt2wee53HllVdizpw5UKvVAICrr7467YESkmkKCgpgtVppafgQitRiXHmCDq/stuPpz7txapkSEloaPii5XB5fGj5q1ChaGk7IETiitawqlQrz589PdyyEjCgikQjFxcXYt28fQqEQxGKx0CFlrNkT8vD+fhdM7jDe+saB2RNoaXgiGo0G3d3dyM/PR14etRUhqTqi5Gbfvn3YuHEjLBZLv2qt9913X1oCI2QkiO07ZbVaodfrhQ4nY8nFLK6fnI9Htljwypc2nD9ajXwF1QkaTCxRNplM0Gg04DhO4IgIGVlSfnV5+umncdNNN8FgMKCoqKhPlynDMJTckJzCMEx8abjf76fClQmcN0qNd/Y4sbc7gLVfWPGbM2nRQSKxpeHd3d20QIOQFKU8C/LBBx/EQw89BLPZHC/iF/v64osvhiNGQjKaWq1GYWEhPB4PrQZKgGX+uzS8fr8b+7v9AkeU2Q5fGh4I0DJ6QlKRcnJjt9tx1VVXDUcshIxYRqMRCoUCXi9tFJnI8YVyTK9RgQewajstDR+KSqWC1+tFZ2en0KEQMqKknNxcddVV2LBhw3DEQsiIJZPJUFRUBL/fT7uGD2HeZAOkHINvOv3Y2kLJYCKxpeGdnZ3weDxCh0PIiJHynJvRo0fj3nvvxbZt23DSSSf1WyFy6623pi04QkYSWhqenAKlCP/z/dLwZ3dYcVq5EiKWljsPJlYssrOzE0qlkpaGE5KElJOb1atXQ6VSYdOmTdi0aVOfnzEMQ8kNyVmHLw2nXcMTu+qkPKzf50K7K4T1e524dLxO6JAymlqtRldXFwwGAyXOhCQh5Vffpqam4YiDkKyQl5eHvLw8OBwOqk+SgFLC4upJevx1Wxde3GXHeaM0UEqoyvNgJBIJvF4vTCYT1Go1VcQmZAj0DCEkjWK7hjMMQ5sfDmHWWA3KNGI4/RG8+pVd6HAynlqths1mg91ObUXIUJJKburq6uKrQOrq6hJ+pWrlypWorq6GTCZDbW0ttmzZkvD8QCCAxYsXo7KyElKpFKNGjcKaNWtSflxChotGo4HBYIDb7RY6lIwmYhnMm5wPAHjzGwe6vGGBI8psIpEIIpEIJpMJ4TC1FSGJJDUs1dDQgFAoFP93uqxbtw4LFy7EypUrccYZZ+Cpp57CrFmz8O2336KiomLA+/z0pz9FZ2cnnn32WYwePRoWi4We6CSjMAwDo9EIu91Ohf2GMKVCiROMMnzT6cffvrDijmlUrC6RWO9NV1cXiouLhQ6HkIyVVHLz0UcfDfjvo7V8+XLMmzcPN9xwAwBgxYoVeP/997Fq1SosW7as3/nvvfceNm3ahMbGxnip+6qqqoSPEQgE+hTAcrlcAIBQKBRP2NIldr10XzcbZXtbSaVS5Ofno62tDSKR6KhWuEQikT7fs828U/JQ964J/zngxmXj1KjRS4/4WtneVgDihf00Gg0kEslRXSvbn4fpRG2VvOFqq1Sux/ApVtG6/vrr8fjjj8d3AI/xer245ZZbkh4iCgaDUCgUeO2113DFFVfEj992223YtWtXv5VYALBgwQLs27cPkydPxgsvvAClUolLL70UDzzwAORy+YCPs2TJEixdurTf8ZdffhkKhSKpWAkhw2vtPhYNVhbjtFHcdDzVCSKE9Ofz+TBnzhw4nU5oNJqE56ac3HAcB5PJhMLCwj7Hu7u7UVRUlPQQUUdHB0pLS/HJJ59g6tSp8eMPP/ww/va3v2Hv3r397nPhhRdi48aNOP/883Hfffehu7sbCxYswLnnnjtoUjVQz015eTm6u7uHbJxUhUIh1NfXY8aMGbRD9BBypa3MZjOam5uRl5d3xCtcIpEIGhsbUVNTk7UbKJrcIdz4dhvCUeDB8404peTIPnjkQlsBgN/vRyQSwdixYwf9YJeMXHkepgO1VfKGq61cLhcMBkNSyU3SS8FdLhd4ngfP83C73X3mEUQiEaxfv75fwpOMH3bX8zw/aBd+NBoFwzB46aWX4rUeli9fjiuvvBJPPPHEgE9yqVQKqbR/N7dYLB62P9DhvHa2yfa2Kioqgt1uR09Pz1En0xzHZe0bdpmOwyXjtHjrWyee3WnHKaUqcEdR2C+b2woAlEolrFYruru7UV1dfdSF/bL9eZhO1FbJS3dbpXKtpD9K6nQ66PV6MAyD4447Ll7PIy8vDwaDAddffz1+/etfJ/3ABoMBHMfBbDb3OW6xWAbdAbe4uBilpaV9iliNHz8ePM+jra0t6ccm5FgRi8UoLi5GKBTK6nkg6TBnoh5KCYsmexAfHKSVZkNRqVTo7u6mVXmEDCDpnpuPPvoIPM/j3HPPxRtvvBGf0Av0TnCrrKxESUlJ0g8skUhQW1uL+vr6PnNu6uvrcdlllw14nzPOOAOvvfYaPB4PVCoVAGDfvn1gWRZlZWVJPzYhx5Jer4dOp4PL5aLCfgloZBx+PiEPz+yw4m9fWHFWtQoyEZXiGoxUKu1T2I+2ZSDkv5JObs4++2wAvRWKKyoq0vJEqqurw9VXX43JkydjypQpWL16NVpaWjB//nwAwKJFi9De3o7nn38eADBnzhw88MADuO6667B06VJ0d3fjzjvvxPXXX39U486EDCeO41BcXAyXy4VQKERd2glcNl6Ld75zotMTxj++deBnE/RD3ymHHV7Y7/APnITkupQ/Fu3ZsweffPJJ/PYTTzyBSZMmYc6cOSlXzpw9ezZWrFiB+++/H5MmTcLmzZuxfv16VFZWAgBMJhNaWlri56tUKtTX18PhcGDy5Mn4xS9+gUsuuQR//vOfU/01CDmmYsO6NISQmETE4ppTegv7vfqVA64ADeUlIhaLwbIsTCYTDXsScpiUk5s777wzXivmq6++Ql1dHS666CI0NjYeUYXiBQsWoLm5GYFAADt37sRZZ50V/9natWuxcePGPuePGzcO9fX18Pl8aG1txaOPPkq9NiTjMQyDoqIicBzXZ/Ue6W96jQrVeRJ4g1Gs201bDQxFrVbD6XTCZrMJHQohGSPl5KapqQnHH388AOCNN97AJZdcgocffhgrV67Eu+++m/YACckWarUaBoMBHo9H6FAyGsswuL62t/fm7T1OWDxUNC0RjuMgFotpWwZCDpNyciORSODz+QAA//nPfzBz5kwAvZMmYz06hJD+YtsySCQS+P1+ocPJaKeWKTChSIZQhMeLu6hHYigqlQputxtdXV1Ch0JIRkg5uTnzzDNRV1eHBx54ANu3b8ePf/xjAL2rlmjFEiGJKZVKFBYWwuPxIMX6mTmFYRhcX2sAANQfcOOQg3ZYT4RlWchkMnR2dtJu9ITgCJKbv/71rxCJRHj99dexatUqlJaWAgDeffddXHjhhWkPkJBsU1hYCLlcHu8BJQMbXyjDGRVKRHnguZ1WocPJeEqlEl6vl3pvCEEKS8FjKioq8K9//avf8cceeywtARGS7WQyGYxGI5qbm6FQKKg+SQLX1ubj01YvPm3x4pvOHpxgpMUDg2EYBgqFAmazGXq9nhZakJyWVM/N4XNpXC5Xwi9CyNAKCgrin7TJ4Cp0Eswc3bttxZqdVhrKG4JcLoff74fFYhE6FEIElVRyk5eXF3+y6HS6PlsvxL5ixwkhQ5NIJDAajQgEAohGaRfsRH55sh4SjsHXnX5sb6OhvEQYhoFKpUJXVxclziSnJTUs9eGHH8arX3744YfUjU5IGhQUFKCrqwsejyftO9RnkwKlCJeN1+K1rx14bqcVk0sVR7WpZraTyWTwer3o7OxMy6aahIxESSU3sa0XAOCcc84ZrlgIySkikQjFxcXYv38/IpFIVu9ifbR+OiEP7+5zockexEeNbpw/mpLBRGKbahoMBkqcSU5KebXUvffeO2CZb6fTiZ///OdpCYqQXBHbVJO2ZUhMI+Xw05N6h72fb7AhGKG5N4lIpVJEIhGYzWaap0RyUsrJzfPPP48zzjgDBw8ejB/buHEjTjrpJDQ3N6czNkKyHsdxKCoqQjQapeqyQ7jseC3yFRw6PWH8+zun0OFkvNimmg6HQ+hQCDnmUk5udu/ejaqqKkyaNAlPP/007rzzTsycORPXXnstPv744+GIkZCsRptqJkcmYvHLSb1z/1750gZfiCZiJxLbfd5sNtOkdZJzUq5zo9Vq8fe//x2LFy/GjTfeCJFIhHfffRfnnXfecMRHSNZjWRZFRUWw2+0IBoOQSCRCh5SxLhijwRtfO9DmCuHNrx345cl6oUPKaGq1Gg6HA3a7Hfn5+UKHQ8gxk3LPDQD85S9/wWOPPYaf//znqKmpwa233oovv/wy3bERkjM0Gg0MBgP13gyBYxnMPaU3oXnjGzuc/v7z/8h/iUQicBwHs9k84FxJQrJVysnNrFmzsHTpUjz//PN46aWX0NDQgLPOOgunn346/vjHPw5HjIRkvdimmmKxGIFAQOhwMtq0KhVG66XwhXj8fbdd6HAynkqlgtPphM1GG5CS3JFychMOh7F7925ceeWVAHorYq5atQqvv/46bcFAyFFQqVTIz8+Hx+MROpSMxjIMrq3t7b155zsnLJ6QwBFlNo7jIBaLYTKZaNI6yRkpJzf19fUoKSnpd/zHP/4xvvrqq7QERUguivXeSCQS+P1+ocPJaJNLFTjJKEMowuOlL6n3ZigqlQput5s21SQ544jm3GzZsgW//OUvMWXKFLS3twMAXnjhBXz33XdpDY6QXKNUKlFQUACPx0P1SRJgGAbX1fZOkN2w34U2Z1DgiDIby7KQyWTo7OxEMEhtRbJfysnNG2+8gQsuuAByuRwNDQ3x+QFutxsPP/xw2gMkJNcUFhZCJpOhp6dH6FAy2glGOU4rVyDKA3/7guaTDCW2USv13pBckHJy8+CDD+LJJ5/E008/Ha+jAABTp07FF198kdbgCMlFcrkcRqMRPp+Pem+GcO0p+WAAbG72YH83DeUlwjAM5HI5Ojs7adiTZL2Uk5u9e/firLPO6ndco9FQJUxC0qSgoAAKhYJ2dh5CjV6K6TUqAMBz1HszJIVCgZ6eHuq9IVkv5eSmuLgYBw4c6Hf8448/Rk1NTVqCIiTXSaVSGI1G+P1+qi47hKtPzgfHADvbffjKTEN5iTAMA6VSCYvFAp/PJ3Q4hAyblJObG2+8Ebfddhs+++wzMAyDjo4OvPTSS7jjjjuwYMGClANYuXIlqqurIZPJUFtbiy1btgx67saNG8EwTL8vmshMslFBQQFUKhW9CQ2hRCPGrLG9O1+vbbCDRvISk8lk8Pv9sFgsQodCyLBJefuFu+66C06nE9OnT4ff78dZZ50FqVSKO+64AzfffHNK11q3bh0WLlyIlStX4owzzsBTTz2FWbNm4dtvv0VFRcWg99u7dy80Gk38dkFBQaq/BiEZTywWo6ioCPv27RM6lIw3Z6Ie9fvd2NMVwDd6BscJHVAGYxgGKpUKXV1d0Ol0QodDyLA4oqXgDz30ELq7u7F9+3Zs27YNXV1deOCBB1K+zvLlyzFv3jzccMMNGD9+PFasWIHy8nKsWrUq4f0KCwtRVFQU/+I47kh+DUIyXn5+fp9EngwsXyHCZcdrAQD/amURpe6bhGQyGUKhEPXekKyVcs9NjEKhwOTJk4/4gYPBIHbu3Im77767z/GZM2di69atCe978sknw+/34/jjj8c999yD6dOnD3puIBDoU87e5XIBAEKhEEKh9FY2jV0v3dfNRtRWycvPz0dzczO11RB+crwG//7OCZMP+PCgG+ePpqQwEaVSie7ubgD0PEwGvWYlb7jaKpXrHXFyc7S6u7sRiURgNBr7HDcajTCbzQPep7i4GKtXr0ZtbS0CgQBeeOEFnHfeedi4ceOAK7gAYNmyZVi6dGm/4xs2bIBCoTj6X2QA9fX1w3LdbERtlbzm5mahQ8h45xQx+Hcrh7/t7EZppBOiI+qbzj30PEwetVXy0t1Wqcw/FCy5iWEYps9tnuf7HYsZO3Ysxo4dG789ZcoUtLa24pFHHhk0uVm0aBHq6urit10uF8rLyzFz5sy0d/eHQiHU19djxowZfWoAkf6orZIXayu9Xg+dTkfDsAmUVISw+c1WWAMMDvJGXDyGem8SCQaDOHToEE477TTk5+cLHU5Go9es5A1XW8VGXpIhWHJjMBjAcVy/XhqLxdKvNyeR008/HS+++OKgP5dKpZBKpf2Oi8XiYfsDHc5rZxtqq+Tp9Xq43W7k5eUJHUrGUkqBC8qieL2Jw9+/cuCC47SQian7ZjASiQQAYLVaUVhYCJalthoKvWYlL91tlcq1BPtLlkgkqK2t7ddtVV9fj6lTpyZ9nYaGBhQXF6c7PEIyTizppzH/xKYU8ihSiWDrieDtPU6hwxkRHA4HFWElWeWIem727duHjRs3wmKx9Cswdt999yV9nbq6Olx99dWYPHkypkyZgtWrV6OlpQXz588H0Duk1N7ejueffx4AsGLFClRVVeGEE05AMBjEiy++iDfeeANvvPHGkfwahIwoGo0Ger0eVqsVer1e6HAylogFfjFRh0c/6carX9lx0VgN1FIaykuEZVmYzWZotVoa9iRZIeXk5umnn8ZNN90Eg8GAoqKiPvNjGIZJKbmZPXs2rFYr7r//fphMJpx44olYv349KisrAQAmkwktLS3x84PBIO644w60t7dDLpfjhBNOwL///W9cdNFFqf4ahIw4DMOgqKgIDocDwWAwPqRA+junWoXXv3HhkCOI17524Ppamk+SiEqlgsPhgN1uh8FgEDocQo5aysnNgw8+iIceegi//e1v0xLAggULBq1svHbt2j6377rrLtx1111peVxCRiK1Wo38/Hx0dnbSBNAEOJbBtafosfRDM/7xrQOXj9dCrxB8/UTG4jgOYrEYJpMJOp0OIhG1FRnZUp5zY7fbcdVVVw1HLISQITAMA6PRCLFYTDs7D2FKhRLjC6QIhHm8/KVd6HAynkqlgsvlgs1GG5CSkS/l5Oaqq67Chg0bhiMWQkgSVCoVDAYDPB6P0KFkNIZhcN33w1Hr9zphctNE7ERYloVUKoXJZKJJ62TES7nvcfTo0bj33nuxbds2nHTSSf2WZt16661pC44QMjCj0QibzYaenh7I5XKhw8lYE4sVqC2RY2dHD15osOGus5IvM5GLlEolbDYburq6UFJSInQ4hByxlJOb1atXQ6VSYdOmTdi0aVOfnzEMQ8kNIceAQqFAQUEBWltbIZPJBi18SYDravOxs6MNHx5046oTdajW9697RXqxLAuZTBaf0zVQjTBCRoKUk5umpqbhiIMQkqLCwkJ0d3fD5/NBqVQKHU7GGmOQYVqVEluavVj7hQ1Lz6e6WIkolUpYrVZ0dXWhrKxM6HAIOSJHVcSP53nwtPsuIYKQyWQwGo3o6emh5+EQrjklHywDbGv14ltLj9DhZDSGYaBQKNDZ2YmeHmorMjIdUXLz/PPP46STToJcLodcLseECRPwwgsvpDs2QsgQCgoKoFQq4fV6hQ4lo5VrJZj5/S7hz+20UjI4BLlcDr/fD4vFInQohByRlJOb5cuX46abbsJFF12EV199FevWrcOFF16I+fPn47HHHhuOGAkhg5BIJDAajfD7/f2qhZO+fjEpD2KOwW6zHzvak99dOBcxDAOVSoWuri5KnMmIlPKcm7/85S9YtWoV5s6dGz922WWX4YQTTsCSJUtw++23pzVAQkhiBQUF6OrqgsfjSftO99mkUCXGpeO0eOMbB9bstKK2VAGWJmIPSiaTwev1orOzE9XV1TRpnYwoKffcmEymATe2nDp1KkwmU1qCIoQkTyQSobi4GKFQCJFIROhwMtrPJuRBIWbRaAtiYyPVCRqKSqVCd3c33G630KEQkpKUk5vRo0fj1Vdf7Xd83bp1GDNmTFqCIoSkRq/XQ6vV0pvQEDQyDj89SQcAeL7BilCE5t4kIpVKEYlE0NnZSfOUyIiS8rDU0qVLMXv2bGzevBlnnHEGGIbBxx9/jA8++GDApIcQMvw4jkNxcTH27t2LcDhMewMlcMXxOvxzjxMmdxjr9zpx2fE6oUPKaGq1GlarFQUFBdDpdEKHQ0hSUu65+clPfoLPPvsMBoMB//jHP/Dmm2/CYDBg+/btuOKKK4YjRkJIEnQ6HfR6PfXeDEEmZvGLSXoAwMtf2uEL0UTsRGJV6M1mM01aJyPGEX28q62txYsvvpjuWAghR4FlWRiNRtjtdgSDQUgkEqFDylgXHqfBm9840O4K4c2vHfjlyXqhQ8poarUadrsddruddqMnI0JSPTcul6vPvxN9EUKEo9VqYTAYqPdmCCKWwbWn9L5Jv/61HY6esMARZTaRSASO42A2m2nSOhkRkkpu8vLy4sWcdDod8vLy+n3FjhNChMMwDIxGI8RiMQKBgNDhZLRpVUocZ5CiJ8zjlS/tQoeT8VQqFZxOJ6xWq9ChEDKkpIalPvzwQ+j1vd22H3300bAGRAg5Omq1GgaDASaTiTY+TIBhGFxfm4+73+/Av/Y6ccUJOhSpxUKHlbE4joNYLIbZbEZeXl58Lg4hmSip5Obss8+O/7u6uhrl5eX9CjrxPI/W1tb0RkcIOSJGoxE2mw09PT2Qy+VCh5OxTi5R4JQSOb7o6MHzDTbcdZZR6JAymkqlgs1mQ1dXF0pKSoQOh5BBpbxaqrq6Gl1dXf2O22w2VFdXpyUoQsjRUSgUKCgogNfrpfokQ7i+tnfuzYcH3Wi00VBeIizLQiaTobOzk4Y9SUZLObnheX7AMtwejwcymSwtQRFCjl5hYSEUCgV8PtpHKZExBhnOrlaBR++mmiQxpVIJn8834IdcQjJF0kvB6+rqAPSOU997771QKBTxn0UiEXz22WeYNGlS2gMkhBwZmUwGo9GIpqYmyOVysGzKn2VyxjWn6PFxswfb23zYbe7BhCIayhsMwzBQKBTo7OxEfn4+DXuSjJR0ctPQ0ACgt+fmq6++6lNDQyKRYOLEibjjjjvSHyEh5IgZDIb4zs5qtVrocDJWqUaCi8Zq8c53Tjz9eTcev7iMNtVMQC6Xw2q1wmw203QEkpGSTm5iq6Suu+46PP7447T7MCEjgEQiQXFxMQ4cOIBIJAKO44QOKWP9YlIe/nPAhX3dAWxu8uCcGkoGB8MwTHxTTYPBQIkzyTgp91OvWLEC4XD/glc2m+2IivitXLkS1dXVkMlkqK2txZYtW5K63yeffAKRSERDYYQMIT8/nzbVTEKeXISfntRbq2vNTiuCtKlmQjKZDKFQCGazmSatk4yTcnLzs5/9DH//+9/7HX/11Vfxs5/9LKVrrVu3DgsXLsTixYvR0NCAadOmYdasWWhpaUl4P6fTiblz5+K8885L6fEIyUWxTTWj0eiAH0zIf/3PCTro5Rw6PWG8s8chdDgZT6PRwGq1wuFwCB0KIX2knNx89tlnmD59er/j55xzDj777LOUrrV8+XLMmzcPN9xwA8aPH48VK1agvLwcq1atSni/G2+8EXPmzMGUKVNSejxCclVsU03aIiUxmZjFNd9vy/Dyl3a4A7TVQCKxQn4mk4m2ZSAZJeWNMwOBwICf/kKhEHp6epK+TjAYxM6dO3H33Xf3OT5z5kxs3bp10Ps999xzOHjwIF588UU8+OCDScV7eD2G2It7KBRCKBRKOt5kxK6X7utmI2qr5KWrrQoKCmC329HT05O1m2rG3mCP5o323GoF3vxGjEOOEF7eZcUNk7N3o8h0tJdSqYTNZoPFYoHBYEhXaBmHXrOSN1xtlcr1Uk5uTj31VKxevRp/+ctf+hx/8sknUVtbm/R1uru7EYlEYDT2rQhqNBphNpsHvM/+/ftx9913Y8uWLRCJkgt92bJlWLp0ab/jGzZs6LOcPZ3q6+uH5brZiNoqedRWyWtsbDyq+19YxOApB4e39zhxotSG/Cwv4XW07QUgZ/acoudh8tLdVqnU7Eo5uXnooYdw/vnn48svv4zPefnggw/w+eefY8OGDalebsBtHAYqEhiJRDBnzhwsXboUxx13XNLXX7RoUbxGD9Dbc1NeXo6ZM2emfcVXKBRCfX09ZsyYQfuuDIHaKnnpbCufz4e9e/dCJBJlZdHNSCSCxsZG1NTUHNXKsNE8j88cZuwy+7HJrsVvzypMY5SZI13tFY1GYbfbUVVVhaKiojRGmDnoNSt5w9VWqQyrp5zcnHHGGfj000/xpz/9Ca+++irkcjkmTJiAZ599FmPGjEn6OgaDARzH9eulsVgs/XpzAMDtdmPHjh1oaGjAzTffDKD3CcXzPEQiETZs2IBzzz233/2kUumAmweKxeJh+wMdzmtnG2qr5KWjrbRaLYqLi9HS0gKFQjHgB4lswHHcUS97/9WPCnDzP1uxqdmLK08K4ThD9iWDMUfbXhzHQaFQoLu7G4WFhVm9YSu9ZiUv3W2VyrVSTm4AYNKkSXjppZeO5K5xEokEtbW1qK+vxxVXXBE/Xl9fj8suu6zf+RqNBl999VWfYytXrsSHH36I119/nQpJEZKkwsJCdHd3w+v1QqVSCR1OxhqdL8W5o9T44KAbT3/ejT9eWJq1yWA6KJVKWK1WdHZ2oqKiQuhwSI47ouQmpqenp98En1SGeurq6nD11Vdj8uTJmDJlClavXo2WlhbMnz8fQO+QUnt7O55//nmwLIsTTzyxz/0LCwshk8n6HSeEDE4qlaKoqAiNjY1QKBS0LUMC156ix+ZmD3ab/fis1YfTK5RCh5SxGIaBUqmExWJBfn4+lEpqKyKclF/VfD4fbr75ZhQWFkKlUiEvL6/PVypmz56NFStW4P7778ekSZOwefNmrF+/HpWVlQB6lxcOVfOGEJK6goICaDQaKuw3hEKVGFccrwUAPLOjG5EoFatLRC6XIxgMUmE/IriUk5s777wTH374IVauXAmpVIpnnnkGS5cuRUlJCZ5//vmUA1iwYAGam5sRCASwc+dOnHXWWfGfrV27Fhs3bhz0vkuWLMGuXbtSfkxCcp1IJEJxcTHC4TAV9hvCzybkQSNl0eoM4b19VCdoKGq1Gt3d3VRTiQgq5eTmnXfewcqVK3HllVdCJBJh2rRpuOeee/Dwww8f9TwcQsixo9frodfrqfdmCEoJh19M0gMAnm+wwRuMChxRZpNIJOB5HmazGdEotRURRsrJjc1mi0/e1Wg0sNlsAIAzzzwTmzdvTm90hJBhw7JsfNluMBgUOJrM9uOxWpRqxHD4I/j7bpvQ4WQ8tVoNm80Gu90udCgkR6Wc3NTU1KC5uRkAcPzxx+PVV18F0Nujo9Pp0hkbIWSYabVaFBQUUO/NEMQcg/89tbf67lvfONDhoiq1iYhEInAch46ODqroSwSRcnJz3XXX4csvvwTQu5opNvfm9ttvx5133pn2AAkhw4dhGBQXF0MqlaZU/TMXnVauwCklcoSiwNOfdwsdTsZTq9VwuVzo6uoSOhSSg1JeCn777bfH/z19+nR899132LFjB0aNGoWJEyemNThCyPBTKBQwGo04dOgQZDIZLQ0fBMMwuPFHBtz0diu2tnixq8OHSSXDs4VLNmBZFnK5HCaTCXl5eZDL5UKHRHJISq9ioVAI06dPx759++LHKioq8D//8z+U2BAyghUWFkKtVsPj8QgdSkarypPi4nG9S8Of3E5Lw4eiUCjg9/tpaTg55lJKbsRiMb7++muq0klIlpFIJCguLkYoFDqqHaJzwdWT9FBJWDTZg3iXloYnxDAM1Go1urq6aF4XOaZS7n+eO3cunn322eGIhRAioPz8fOj1eqpPMgSNjMPck79fGv6FFZ4AJYOJSKVSRKNRdHR00NJwcsykPOcmGAzimWeeQX19PSZPntyvxPby5cvTFhwh5NhhWRbFxcVwOp0IBoOQSCRCh5SxfjxOi39950SLM4SXvrThxh8VCB1SRostDbfZbDAYDEKHQ3JAysnN119/jVNOOQUA+sy9AUDDVYSMcBqNBgUFBTCbzdDr9fScHoSIZXDjaQYs3mDC2986cdFYLcq1lAwORiQSQSwWo6OjA1qtlnbVJsMu6eSmsbER1dXV+Oijj4YzHkKIgGJLw+12O3p6eqBQ0GqgwUwuVeJHZQpsb/Nh9fZuPDCjROiQMppKpYLNZoPFYkFpaanQ4ZAsl/ScmzFjxvSpVzB79mx0dnYOS1CEEOHI5XIUFRXB5/PRHIkh/O+PDOAYYHubD5+3eYUOJ6OxLAuFQgGz2Uw1lciwSzq5+eEyvvXr18PrpSczIdmIloYnp1wrwWXf7xr+1PZuhGlpeEJyuRyBQICWhpNhR9W6CCH9iMVilJSUIBQK0a7hQ/jFRD203+8a/s89TqHDyWiHLw2nVXlkOCWd3DAM029yIU02JCR76fV65Ofn05vQEFRSDtfW5gMAXmiwwuqjZDCR2K7hHR0dVFOJDJukJxTzPI9rr70WUqkUAOD3+zF//vx+S8HffPPN9EZICBFEbGm4y+WC3++HTCYTOqSMdcEYDd7b58Le7gCe/rwbd59dJHRIGU2j0cBut6Orqyu+Mz0h6ZR0z80111yDwsJCaLVaaLVa/PKXv0RJSUn8duyLEJI9NBoNjEYjPB4PTS5OgGMZ3DKlAAyAjxo92NVBE2YT4TgOUqkUJpMJfr9f6HBIFkq65+a5554bzjgIIRmqqKgIDocDHo8HGo1G6HAy1hiDDBeP0+Kd75z467YurLqsAmKOhu4Ho1QqYbVaYTKZUFVVRdMcSFrRhGJCSEISiYQmFyfp2lP00Mk4tDpDePMbh9DhZDSGYaDRaGCxWOB00kRskl6U3BBChqTX62EwGOhNaAgqKYdfndo7ufilL23o9IQEjiizxSYXt7e3U+JM0oqSG0LIkFiWRUlJCSQSCXp6eoQOJ6OdN0qNk4wyBMI8nvysW+hwMp5Go4HD4ehTJJaQo0XJDSEkKSqVCkVFRfB6vTS5OAGGYXDzlAJwDLC1xYvPWqnYaSIcx0Eul8NkMlHiTNKGkhtCSNKMRiPUajXcbrfQoWS0qjwprjhBBwBYua0LgTAlg4koFAr4/X50dHRQ5WKSFoInNytXrkR1dTVkMhlqa2uxZcuWQc/9+OOPccYZZyA/Px9yuRzjxo3DY489dgyjJSS3icVilJaWIhKJIBSi+SSJ/HKSHgYFB7MnjHVf2YUOJ6MdXrnY4XAIHQ7JAoImN+vWrcPChQuxePFiNDQ0YNq0aZg1axZaWloGPF+pVOLmm2/G5s2bsWfPHtxzzz245557sHr16mMcOSG5Kza5mCoXJyYXs5h/WgEA4NXddrQ7gwJHlNkkEgkYhqHJxSQtkq5zMxyWL1+OefPm4YYbbgAArFixAu+//z5WrVqFZcuW9Tv/5JNPxsknnxy/XVVVhTfffBNbtmzB//7v/w74GIFAAIFAIH479oIcCoXS/skzdj36RDs0aqvkZWJbFRYWxmvfyOVyocOJi5Xzz5Sy/lPKZKgtkWNnRw/++mkXHjjfmFH1XDKtvZRKJex2Ozo6OlBcXCx0OH1k4vMwUw1XW6VyPYYXaIAzGAxCoVDgtddewxVXXBE/ftttt2HXrl3YtGnTkNdoaGjArFmz8OCDD8YTpB9asmQJli5d2u/4yy+/DIVCceS/ACGEJKGrB/jDlxzCPINfjo7g1AKaU0LIkfD5fJgzZw6cTueQBUUF67np7u5GJBKB0Wjsc9xoNMJsNie8b1lZGbq6uhAOh7FkyZJBExsAWLRoEerq6uK3XS4XysvLMXPmzLRXWw2FQqivr8eMGTMgFovTeu1sQ22VvExtq3A4jAMHDsDtdkOn0wkdDoDeHojGxkbU1NSA4zihwwEAjAHwC8aBvzXY8XaLGBedXAadPDNiy8T2AgCbzQa9Xo9Ro0aBZQWfGgogc5+HmWi42iqVoXBBh6WA/juL8zw/ZLftli1b4PF4sG3bNtx9990YPXo0fv7znw94rlQqjW/2eTixWDxsf6DDee1sQ22VvExrK7FYjIqKCuzduxehUCijNtbkOC6j3qx/OkGPjw95cdAWxFM7bPjdOZm1WWSmtVdeXh4cDgecTicKCwuFDqePTHseZrJ0t1Uq1xIsJTYYDOA4rl8vjcVi6deb80PV1dU46aST8Ktf/Qq33347lixZMoyREkIGo9VqUVRURBtrDkHEMrj9jEKwDLCpyYOthzxCh5TRRCIRpFIp2tvbqfYNOSKCJTcSiQS1tbWor6/vc7y+vh5Tp05N+jo8z/eZMEwIObaKioqg0Who9dQQxhhkuOpEHQDgL592wRPIjEm8mUqpVMLn86G9vZ1q35CUCTqYWVdXh2eeeQZr1qzBnj17cPvtt6OlpQXz588H0DtfZu7cufHzn3jiCbzzzjvYv38/9u/fj+eeew6PPPIIfvnLXwr1KxCS8yQSCcrKyuiDRhJ+MUmPMo0Ytp4Int5hFTqcjMYwDLRaLbq6umC1UluR1Ag652b27NmwWq24//77YTKZcOKJJ2L9+vWorKwEAJhMpj41b6LRKBYtWoSmpiaIRCKMGjUKf/jDH3DjjTcK9SsQQgDodDoYjUZ0dHTE65WQ/qQiFgvPKMQd77bjvX0unFOtwskltGpzMGKxGCKRCO3t7VCr1QPOnyRkIIJPKF6wYAEWLFgw4M/Wrl3b5/Ytt9yCW2655RhERQhJBcMwKCkpgdvthtvtTvtKxGxyUpEcl4zT4p3vnFjxiQVPXV4BmTgzVgRlIrVaDavVivb2dlRXV1PiTJJCzyhCSFpIpVKUlpYiHA5TobMhXD85HwVKEcyeMP7WQEMuiTAMA41GA4vFArudtrEgyaHkhhCSNnq9HoWFhXA6nTQJNAGFmMVtU3u3ZnjrGyf2WPwCR5TZJBIJWJZFW1sbgkHaxoIMjZIbQkjaMAyD0tJSKJVKeDy03DmRU8uUOG+UGjyA5Z9YEIxQMphIbDd6k8lEiTMZEiU3hJC0kslkKC0tRTAYpOGpIcz/kQE6GYcWRxAv7bIJHU5GY1kWarUaZrMZTqdT6HBIhqPkhhCSdvn5+SgoKIDL5aJP2QloZBxunvL9zuFf2fF1JxWsSyS2WoqGp8hQKLkhhKQdy7IoKyuDQqGA2+0WOpyMNq1KhRmj1YjywB83d8IbpErPiWg0GjidTiruRxKi5IYQMizkcjnKysoQiUToU/YQbjqtAEaVCJ2eMFZ91iV0OBmNZVloNBqYzWYq7kcGRckNIWTY5Ofnw2g0wuVy0d5TCSglLO46ywiWAeoPuLG5iSZjJyKRSCASidDW1kZ7T5EBUXJDCBk2sdVTWq2WJoEO4USjHLNPygMAPL7Vgi5vWOCIMptarYbX60VbWxslzqQfSm4IIcNKIpGgvLwcHMfRp+wh/PJkPY4zSOEJRvHIlk5EaU7JoBiGgU6ng8ViQWdnp9DhkAxDyQ0hZNhptVoUFxfD6/UiEqHdsAcjYhncdZYRUhGDXaYe/ONb6u1KRCQSQS6Xo6Ojg+oqkT4ouSGEHBPFxcXIz8+Hw+EQOpSMVq6V4MZTDQCANTu60WSjndYTUSqVCAaDaGlpQThMQ3mkFyU3hJBjguM4lJeXQyqV0qfsIVw0VoPTyhUIRYH/29yJYJjmlCSi1Wpht9vR0dFBy8MJAEpuCCHHkFKpRGlpKQKBAFUvToBhGNSdUQidjEOTPYjnvqAlz4lwHAeVSgWTyUQ9gwQAJTeEkGOssLCQNtdMgk4uwu1nFgIA3vzGiU9bvAJHlNlkMhkAoLW1FYEADeXlOkpuCCHHVKx6sVqtpuXhQzi9XInLx2sBAH/a0gmTm3q7EtFoNHC73WhtbaXl4TmOkhtCyDEnk8lQUVEBlmXh8/mEDiej3XCqAeMLZPAGo3jwQzPNv0mAZdn48nCTySR0OERAlNwQQgSh0+lQWlqKnp4emn+TgJhjsHh6EbRSFgdsATzxWbfQIWU0kUgEhUKB9vZ22O12ocMhAqHkhhAimKKiIhiNRjidThpGSKBAKcJvzy4CA+C9fS5s2O8SOqSMJpfLAQAtLS1UODJHUXJDCBEMy7IoLy+HTqejVS5DqC1V4OqT9QCAv3zahYNU/yYhjUYDr9dL9W9yFCU3hBBBSSQSVFRUQCKRUP2bIfx8Yh5OLVUgGOHx4IdmeINU7Xkwse0Zuru7qf5NDqLkhhAiOLVajfLycgSDQVrGmwDL9G7PUKgUocMdwiNbLPSmnQDHcVCr1ejo6IDVSrWCcongyc3KlStRXV0NmUyG2tpabNmyZdBz33zzTcyYMQMFBQXQaDSYMmUK3n///WMYLSFkuBQUFKCkpARut5v2n0pAI+Nwz/QiiFlga4sXr3/tEDqkjCaVSiESidDS0gKvl2oF5QpBk5t169Zh4cKFWLx4MRoaGjBt2jTMmjULLS0tA56/efNmzJgxA+vXr8fOnTsxffp0XHLJJWhoaDjGkRNC0o1hGJSWlsb3n6IeicGNLZDhxtMKAABrdlqx20yTZhNRqVQIBAI4dOgQrczLEYImN8uXL8e8efNwww03YPz48VixYgXKy8uxatWqAc9fsWIF7rrrLpx66qkYM2YMHn74YYwZMwbvvPPOMY6cEDIcRCIRKioqoFAo4HLRiqBELh6rwXmj1IjywAMfmtDuCgodUsaKzb+x2+1U4C9HiIR64GAwiJ07d+Luu+/uc3zmzJnYunVrUteIRqNwu93Q6/WDnhMIBPqM4cdeMEOhUNoz+Nj16JPB0KitkpdrbSUWi1FaWorGxka43W4oFIqk7xsbzsqVYa1fn6ZHqyOAfdYg7q03YfmsYqilXNL3z7X2is2/EYlEKCoqAsMwSd83156HR2O42iqV6wmW3HR3dyMSicBoNPY5bjQaYTabk7rGo48+Cq/Xi5/+9KeDnrNs2TIsXbq03/ENGzak9KKZivr6+mG5bjaitkoetVXyGhsbhQ7hmLm6CnjUzaHdFcI97zZj/vgoRCn2yedSewE4qsnF9DxMXrrbKpVq5oIlNzE/zJx5nk8qm37llVewZMkSvP322ygsLBz0vEWLFqGuri5+2+Vyoby8HDNnzoRGoznywAcQCoVQX1+PGTNmQCwWp/Xa2YbaKnm52lY8z6O9vR1tbW3QarUQiYZ+uYpEImhsbERNTQ04LvkejJHuoZIg7ni3A/tdLN63anDbFENSr6O52l6xicU1NTVJvw/k6vPwSAxXW6UyVC1YcmMwGMBxXL9eGovF0q8354fWrVuHefPm4bXXXsP555+f8FypVAqpVNrvuFgsHrY/0OG8drahtkpeLrZVRUUFIpEIOjs7kZeXl/QbMMdxOfVmPdogx+/OKcLvPzBhwwEPynVS/PSkvKTvn2vtpdFo4HA40NbWhjFjxqTUi5+Lz8Mjle62SuVagk0olkgkqK2t7ddtVV9fj6lTpw56v1deeQXXXnstXn75Zfz4xz8e7jAJIQLiOA6VlZXQ6/W0gmoIPypXYv6PDACANTus+LiZCiImotVq4fV60dzcjGCQJmNnG0FXS9XV1eGZZ57BmjVrsGfPHtx+++1oaWnB/PnzAfQOKc2dOzd+/iuvvIK5c+fi0Ucfxemnnw6z2Qyz2Qyn0ynUr0AIGWZisRhVVVVQqVS0RcMQLjteh8vGa8ED+OPmTuzr9gsdUsZiGAZ5eXmw2+04dOhQzkyqzhWCJjezZ8/GihUrcP/992PSpEnYvHkz1q9fj8rKSgCAyWTqU/PmqaeeQjgcxq9//WsUFxfHv2677TahfgVCyDEgl8tRVVUFsVgMt9stdDgZ7cYfGXBqqQKBCI/f/8cEi4dW9wyGZVnodDpYLBa0tbVRz2AWEXxC8YIFC7BgwYIBf7Z27do+tzdu3Dj8ARFCMpJGo0FVVRUOHDiAnp6e+M7PpC+OZbDonCL8Zn0bmuxB3PcfEx6ZVQpVCkvEc4lIJIJarUZ7ezvEYjFKSkqEDomkgeDbLxBCSLLy8/NRXl4On89He1AloJSwuP/8YujlHJrsQdxTb0JPiArXDUYqlUIul6O1tRXd3d1Ch0PSgJIbQsiIUlxcjLKyMng8HpoImkChSoyHZ5ZAJWGxp8uPJR+YEAxTgjMYhUIBjuPQ1NQEm80mdDjkKFFyQwgZURiGQVlZGUpKSuByuahibALVeikemlkCuYjBLlMPHtpoRjhK80oGo1arAQBNTU00eX2Eo+SGEDLisCyL8vJyFBcXw+l0UoKTwLgCGe6fUQIJx2Bbqw9/2tyJCCU4g9JoNPHihrQSd+Si5IYQMiLFauAUFRXB6XQiHA4LHVLGmlAkx33nFkHEAhubPPjzp120MigBrVaLYDCIpqYmWp03QlFyQwgZsWIJTmFhIRwOB9UqSeDUMiXuPrsILAO8t8+Fp7Z3U4KTgE6ng9/vR2NjY3y7BjJyUHJDCBnRRCIRqqqqYDAYaJ7EEKZVqVB3Ru9efG9968SLXzqEDSiDMQwDnU4Hn8+HgwcPoqenR+iQSAoouSGEjHhisRjV1dXIy+vdTykapVVBg5kxRoNfn967TcMrux14v42hHpxBxKoYezweNDc3Cx0OSQElN4SQrCCRSFBdXQ0AcDgclOAkcOl4Ha6vzQcArG/lsOYLOyU4g4glOLEdqf1+2tJiJKDkhhCSNSQSCYDe+RI2m43m4CQwe0IefjVZDwB44xsnHt/aRauoBhHbpgHoXSZOc3AyHyU3hJCsU11dDYPBALvdTquoErjieC1+PioClgHe3efC/23upDo4g2DZ3rdLl8uFgwcPwuOhXdczGSU3hJCsI5FIUFNTE19FRXVwBnd6IY/fTiuEiAU2NXmw9AMTAlTJeFB5eXnw+Xw4cOBAfKiKZB5KbgghWSk2yThWB4e2ahjctCollpxXDCnHYHubD4s3dMAbpARnILFVVH6/HwcOHIDdbhc6JDIASm4IIVkrtky8pKQEbrebNttM4NQyJR6+oAQKMYuvOv347XvtcPppztJAYpOMw+EwDh48CKvVKnRI5AcouSGEZLVYob/YZpu02mVwJxrl+OOFJdBKWey3BvCb9W0wuWlIbzBarRYA0NjYiK6uLoGjIYej5IYQkvVie1FVVFTA5/PRapcExhhkeOSiMhgUIrQ6Q7j1nVbsNlMBu8Go1WqwLIvGxkZ0dHTQkvoMQckNISQnMAyD0tJSVFdXIxKJwOl00hvRICp0Ejx+cRmOM0jhCkRx93vtWL+XNpEcjEqlglQqRXNzM5qbm2mFXgag5IYQkjMYhkFRURHGjBkDiUQCm81Gxf4GYVCK8MisUpxTrUKEBx7f2oVV26gWzmDkcjnUajU6Ojpw8OBBGv4UGCU3hJCco9PpMGbMmHixP/qkPTCpiMXdZxtxzcm9xf7+sceJe+s74AnQROOBSCQS5OXlwWq1Yv/+/bSjuIAouSGE5CSlUonRo0fDaDTC4XDQJ+1BMAyDOZP0uHd6EaQiBjs7enDbv9vQ5qSl9QPhOA55eXnwer3Yt28frFYrDX8KgJIbQkjOiu1HVV5eDp/PR1VnEzizSoXHLipDgVKENmcIt/2rDZ+30cTsgbAsi7y8PPA8jwMHDqCjo4OGP48xSm4IITmN4ziUl5dj1KhR4Hkedrud3ogGMSpfir9cUobxBTJ4glHcU2/C6u3dCEWoZ2IgarUaUqkUhw4dQnNzMxWSPIYouSGE5DyGYVBYWIgxY8ZAqVTCZrPRlg2DyJOL8MdZpbhsfG+Nlze+ceD2f7ehnYapBiSXy6HRaGAymbBv3z7asuEYETy5WblyJaqrqyGTyVBbW4stW7YMeq7JZMKcOXMwduxYsCyLhQsXHrtACSFZT6vVYuzYsSgqKoLL5aJhqkFIOAYLTi/AkvOKof6+4N+Cf7ai/gC9cQ9ELBZDr9fD4/Fg7969NEx1DAia3Kxbtw4LFy7E4sWL0dDQgGnTpmHWrFloaWkZ8PxAIICCggIsXrwYEydOPMbREkJyQWzTzVGjRgEAbDYbIhFaHTSQKRVKPHlZBSYUyeAP83hkiwV/3NxJ+1INIDYPRyQSobm5mZaLDzNBk5vly5dj3rx5uOGGGzB+/HisWLEC5eXlWLVq1YDnV1VV4fHHH8fcuXPjZa8JISTdGIaB0WjE2LFjodVqYbfbaV+qQRiUIvzhglJcc7IeLAN8cNCNX/+zFXu76I17IAqFAlqtFl1dXdi7dy9sNhutphoGIqEeOBgMYufOnbj77rv7HJ85cya2bt2atscJBAJ9XpRi452hUCjtY+qx69FY/dCorZJHbZW8dLeVVCqND5ubzWb09PRApVKBYZi0XF9osR6pdPRMzT5JixONUvxxiwUmdwgL/92GnxyvxZyJOkhFgs+AOGrpbCuGYaDVauF2u7F3714UFRWhqKgIIpFgb8lpNVyvWalcT7CW7O7uRiQSgdFo7HPcaDTCbDan7XGWLVuGpUuX9ju+YcMGKBSKtD3O4err64flutmI2ip51FbJo7ZKTWNjY1quIwFQdzzwaiOLBiuL175x4qMDDvy0JoqxuuzonUhXWx3OYrFg9+7dab+u0NL9PPT5fEmfK3ia+MNPQDzPp/VT0aJFi1BXVxe/7XK5UF5ejpkzZ0Kj0aTtcYDerLK+vh4zZsyAWCxO67WzDbVV8qitkjfcbRUIBGAymdDV1QWGYeKbJo5UkUgEjY2NqKmpAcdxabvuxPHAtlYvnvjMim5fBCv3cDh/lAo31OqhkaXvcY6l4WorAIhGo3C5XGBZFoWFhTAajZBIJGl9jGNpuJ6Hqaw0Eyy5MRgM4DiuXy+NxWLp15tzNKRSKaRSab/jYrF42N4ohvPa2YbaKnnUVskbrrYSi8UYPXo0DAYD2tvb4XQ6oVKpIJPJ0v5YxxLHcWl/wz6jSoNJJSo8t9OKf33nxH8OevB5mw/zTyvA9JqRO7Q3HG3FcRzy8/PjybPH40FpaSn0ev2IbScg/c/DVK4l2EcOiUSC2traft1W9fX1mDp1qkBREUJIYgzDIC8vD2PHjkVFRQWCwSDsdjutqBqAUsLi5ikFWP7jUlTqJHAGovi/zZ1YXG9Ch4vmkP2QVCqFXq9HIBDAvn370NjYiJ6eHqHDGpEE7U+tq6vDM888gzVr1mDPnj24/fbb0dLSgvnz5wPoHVKaO3dun/vs2rULu3btgsfjQVdXF3bt2oVvv/1WiPAJITlMLBajvLy8z4oqj8dDK18GcHyhHE9cWo5rTtZDzAI723341VuH8ORnXXD5KSk8HMuy0Gg0UKvV6OzsxJ49e2A2myl5TpGgc25mz54Nq9WK+++/HyaTCSeeeCLWr1+PyspKAL1F+35Y8+bkk0+O/3vnzp14+eWXUVlZiebm5mMZOiGEAAA0Gg2USiUsFgvMZjOsViuUSiXkcrnQoWUUMde7Aee0ahVWfdaNne0+vPWtExsOuPHzCXm4bLwWkixYVZUuscJ/Xq8XjY2NsFqtKCoqQl5e3oie53WsCD6heMGCBViwYMGAP1u7dm2/Y/SpiBCSaTiOQ3FxMfLy8mCxWGCxWGC1WqFWq0f0xNDhUK6V4OGZJdjZ7sPTn3ejyR7EMzuseOc7J649JR/n1KjAjuB5JunEMAxUKhXkcjncbjf27dsHvV6P4uJiqNXqET0fZ7gJntwQQki2kMlkqKioQH5+PsxmM7q7u+H1eqFWq7Omhkm61JYqMKm4HB8cdGPtF1Z0esL4v82dePMbB244NR+TioenVMdIxHEcdDodwuEwbDYbnE4nDAYDjEYjlEql0OFlJHq2EUJImimVStTU1MBgMMBkMsFut4PjOKhUqrSvtBnJOJbBzDEanFWtwpvfOPDqbjv2WwP47XsdOMEow88n5GFyqYJ6KL4nEoniE45jf1eFhYUoKCgY8Sv20o2SG0IIGQaxKrRqtRp2ux0mkwlOpxMsy0KlUlFPzmFkIhZzJupx0XEavLjLjvf2OfFNpx/31JswWi/FzybmYWqFEhxLSQ7w3xInPT09aGlpQVdXFwoKCmAwGGiu1/fo2UUIIcOIZVnk5+dDp9PBbrfDYrHA6XTGiwBSkvNfOrkIN08pwM8m5OGNbxz4914nDtgCePAjM8q0YvzspDxMH6WGiJIcAIBcLodMJosnORaLBQaDAQaDIeeHq+hZRQghxwDHcTAYDNDr9XA4HLBYLHA4HAAAlUpFBRoPY1CKcOOPDPjZhDz841sH/rnHiTZnCI98bMHzDTZcfrwWM8ZooJHSEB/DMFAoFJDL5fD7/Whvb4fFYkF+fj4KCgqyai+0VFByQwghxxDLstDr9cjLy4PT6YTFYokXAVQqlZBKpTn5ZjQQrYzDNafk48oT8/DvvU68+Y0DFm8Yqz+3Yu0XNpxdrcLF47QYa6A2YxgGcrkccrkcgUAgPqFdr9cjPz8fGo0mp+Z7UXJDCCECYBgGOp0uvju01WqFzWaDx+OBVCqFQqHIqTejRJQSFj89qbcWzgcH3XjnOycabUHUH3Cj/oAbY/KluHicFufUqCCjWjnxOTnBYBDd3d3o6uqCWq2GwWBAXl7egFsSZRtKbgghREAMw0Cj0UCj0aC4uBh2ux1dXV3xeTkKhSIn3oySIRWxuGisFrOO02BPlx//+s6FzU1u7LcG8NgnFjz9eTfOG6XGeaPUOI56cyCRSCCRSBCJROD1enHw4EHIZLJ4b04218qh5IYQQjKETCZDcXExCgsL4XK50N3dDYfDAbfbDalUCrlcThOQ0ZsQHl8ox/GFcvzvjwzYsN+Ff3/nhNkTxtt7nHh7jxOlGjGm16hwbo0apdrcLqTIcRw0Gg14noff74fZbEZnZyc0Gg0MBgM0Gk3WLSWnZwkhhGQYjuOQl5eHvLw8+Hw+2O12WK1WuN1uRKNRyGQyyOVyKsMPQCfj8NOT8nDliTrsbPfhPwfc+LTFi3ZXCC/usuPFXXYcZ5Bieo0a51SroFfk7tve4fNywuEwPB4PHA4HJBIJtFot8vLyoNFosqKqdu7+LxNCyAigUCigUChQXFwMt9sNp9MJm80WX2kll8shlUpzPtFhGQanlilxapkSvlAUWw958FGjB190+LCvO4B93QE8/Xk3TiiU4fQKJaaUK3O6R0ckEkGr1YLneQSDQVitVnR1dUEmk0Gn0yEvL29ElyoYmVETQkiOYVkWWq0WWq0WJSUlcLvdcDgcsNvtsNvtYBgGUqkUMpks5yciK8Qszh+twfmjNbD3hLG5yYOPGt3Y0xXAV51+fNXpx9OfW1GuFeP0ciVOr1BifIEsJ4sExv5upFIpotEoAoEAOjs7YTabIZfLodPpoNFooFKpRtTcL0puCCFkhBGJRPFhq9LSUng8HrhcLjidTrhcLkSjUUgkEshkspyvn5MnF+Gy43W47HgdzO4QPmv14tNWL3abetDqDKHV6cBrXzuglbKYXKbEySVyTCpWoECZe2+PLMvGh62i0Sh6enpgNpthMpkglUqhVquh0+nim3lm8mTk3PvfI4SQLCKRSKDX66HX6xGJRODxeODxeGC32+Hz+RAOh8FxHKRSKSQSSU736hSpxfFExxuMYEe7D9tavNje5oMzEMUHB9344KAbAFCmEWNisRwnlygwoUgOVY7liCzLQqlUQqlUgud5BAKB+Eo+iUQChUIBrVYLlUoFhUKRcfN0KLkhhJAswXFcn6Ern88X79WJfY9GoxCJRJBKpTmd6CglHM6uVuPsajXCUR7fdPZgZ7sPu0w92G8NoM0VQpsrhH/vdQEAavIkKJWyOJ3z4MQiBYwqUUb3XKQTwzCQyWTxFVWhUAg+ny9erkAikUCpVEKj0UCpVGZEbyElN4QQkoUYhol/8jYajQiHw/D5fPE3Ja/XC4/HAwBwOp2QyWSQSCQjdgLp0RCxDCYWKzCxWAEA8AQi2G3uwZemHjSYenDIEUSjPYhGsNhi7gIA5Mk5jCuQYXyBDOMLZTguXwqZODcmdYvF4ngCE41GEQwG4XK5YLVawbJs/Gc8zwsWY+79FRNCSA4SiUTxYoFFRUUIhUJwuVz4+OOPodPp4Pf74fF4EIlEwDAMxGIxJBIJxGJxzq3EUkk5TK1UYWqlCgBg7wmjod2LT/d3whyS4aAtCHtPBJ+2ePFpixcAwDK9Q1mj8qUYpZeiRi/F6HwptLLs7h1jWbZPr040GoXL1dvbRckNIYSQY0osFkOj0QAARo8eDZZl0dPTA7/fD5/PB5fLhUAgAK/XC57n4wlP7CuXhrTy5CKcXa1CSdiEMWNKEeYZHLAG8K3Fj++6/NjT5YfVF0GLM4QWZwgfNXri9zUoONTopajKk6BCK0G5ToIKnQSKLO3lObznRkiU3BBCCAHHcVCpVFCpensrYvVPAoEA/H4/enp64PF4EAgE0NPTg2g0Gr+fSCSCWCyGSCTKiaRHKmJxglGOE4xyAL1tZeuJoNEWwAFrAAdtQTTaAmh3hdDti6Db58P2Nl+faxgUIlToxKjQSlCmk6BELUaxWgyjSpSTS9LTjZIbQggh/Rxe/yTWwwP0TiaNJTyBQCA+jyc2yTSW9LAsC5FIFE94RCJR1g5vMQyDfIUI+QoRTi1Txo/7QlE02QI4aAugxRFCiyOIVmcQtp4Iun1hdPvC+KKjp8+1WAYwqkQo/j7ZKVaLUaQWo0ApQoFShDw5BzZHJjIfDUpuCCGEJC02LBXr4QF6ey5iSU8wGEQwGERPTw96enoQDAbh9/sRiUTA83x8iCuW9Bz+lW3Jj0Lct4cnxh2IoNUZRKsjhBZnEG3OIDrcYZjdIQQjPEzuMEzuMICeftcUsYgnOgXK3qRHL+egV/R+z//+uyTHd0en5IYQQshRiS0HHqjWSSQSQSgUQjAYRCgU6tPzE0t6gsEgIpFIvNcH6O35iSU8hyc/HMeN+CXYaikX3/jzcFGeh80XgckdQoc7BJMrBJM7BIs3DIsnDFtPGOEoDkt+/IM+hkrCQq8QIU/GQSfnoJX1ful+8F0tZaGScBBzI7tNf4iSG0IIIcMmlpgMtOs0z/MIh8MIh8MIhULxf8d6fwKBAEKhECKRCAKBACKRCCKRCIDehCrWC8Sy7JBfIwHLMDAoRTAoRTipSN7v55EoD6svDIs3jK7vE55uXxg2XxhWXwS2njBsPRGEIjw8wSg8wSBaknxsuYiBWtqb7MS+KyUcVBIWSjELpST2xfV+F7OQi1koxAzkYhYSjsmopFPw5GblypX405/+BJPJhBNOOAErVqzAtGnTBj1/06ZNqKurwzfffIOSkhLcddddmD9//jGMmBBCSDocvgJLLu//Zg70JkCxpCaW/ITD4T7HYr1BsUQoHA4jGo2C53lEo1FEo9H4G29seXIsKUr0/fB/C7msOYZjGRSqxChMUC6Z53sTG6uvN9Fx9ETg9Efg8Pd+/+G/vcEoeAA9YR494TAs3iOLjWV6h+HkYhYyDlAwLC688MiulQ6CJjfr1q3DwoULsXLlSpxxxhl46qmnMGvWLHz77beoqKjod35TUxMuuugi/OpXv8KLL76ITz75BAsWLEBBQQF+8pOfCPAbEEIIGU6x+TmxqspDiQ1vxZKfw//9w9uHJ0uHJ0ShUCg+Pyh2LBwOAwBsNhs4juuX7MSSoWS+fnj+4bcH+lmq7dXb88KhKm/o8yNRHt5QFO5ABO5A3+/eYBSeYBTeUBTeYBTeYOT7773HfKEoAuHedojy+L63qHdoUScRthdH0ORm+fLlmDdvHm644QYAwIoVK/D+++9j1apVWLZsWb/zn3zySVRUVGDFihUAgPHjx2PHjh145JFHBk1uAoEAAoFA/HasuFBs7DedYtdL93WzEbVV8qitkkdtlZpsbq/YcFgqYr08h/f6xL4Hg0FYrVZUV1eD47j4ecDACdXh14lNpo6dH/seS6AO//cPb8f8sOfp8OOD9SodniwN9P3wfysBqCQMiiUA1ADDsADYfuf/8DqRKI9AhEdPiIf/++9OXwB2hwPhcDitQ4Kp/J0yvEB9bcFgEAqFAq+99hquuOKK+PHbbrsNu3btwqZNm/rd56yzzsLJJ5+Mxx9/PH7srbfewk9/+lP4fL4BCwctWbIES5cu7Xf85ZdfhkKhSNNvQwghhJDh5PP5MGfOHDidzj7lCQYiWM9Nd3c3IpEIjEZjn+NGoxFms3nA+5jN5gHPD4fD6O7uRnFxcb/7LFq0CHV1dfHbLpcL5eXlmDlz5pCNk6pQKIT6+nrMmDEjIyo0ZjJqq+RRWyWP2io11F7Jo7bq22t0eA/TD4+Fw2F8+OGHaW+r2MhLMgSfUPzD7q7Y7PdUzh/oeEysCNUPHb7xV7oN57WzDbVV8qitkkdtlRpqr+RRWw0tNnyU7rZK5VqCrY8zGAzgOK5fL43FYunXOxNTVFQ04PkikQj5+fnDFishhBBCRg7BkhuJRILa2lrU19f3OV5fX4+pU6cOeJ8pU6b0O3/Dhg2YPHkyZdKEEEIIASBgcgMAdXV1eOaZZ7BmzRrs2bMHt99+O1paWuJ1axYtWoS5c+fGz58/fz4OHTqEuro67NmzB2vWrMGzzz6LO+64Q6hfgRBCCCEZRtA5N7Nnz4bVasX9998Pk8mEE088EevXr0dlZSUAwGQyoaXlv/UVq6ursX79etx+++144oknUFJSgj//+c9U44YQQgghcYJPKF6wYAEWLFgw4M/Wrl3b79jZZ5+NL774YpijIoQQQshINTI23CCEEEIISRIlN4QQQgjJKpTcEEIIISSrUHJDCCGEkKxCyQ0hhBBCsgolN4QQQgjJKpTcEEIIISSrUHJDCCGEkKwieBG/Yy22i3gqW6cnKxQKwefzweVy0V5XQ6C2Sh61VfKorVJD7ZU8aqvkDVdbxd63Y+/jieRccuN2uwEA5eXlAkdCCCGEkFS53W5otdqE5zB8MilQFolGo+jo6IBarQbDMGm9tsvlQnl5OVpbW6HRaNJ67WxDbZU8aqvkUVulhtoredRWyRuutuJ5Hm63GyUlJWDZxLNqcq7nhmVZlJWVDetjaDQa+uNPErVV8qitkkdtlRpqr+RRWyVvONpqqB6bGJpQTAghhJCsQskNIYQQQrIKJTdpJJVK8fvf/x5SqVToUDIetVXyqK2SR22VGmqv5FFbJS8T2irnJhQTQgghJLtRzw0hhBBCsgolN4QQQgjJKpTcEEIIISSrUHJDCCGEkKxCyc0wufTSS1FRUQGZTIbi4mJcffXV6OjoEDqsjNPc3Ix58+ahuroacrkco0aNwu9//3sEg0GhQ8tYDz30EKZOnQqFQgGdTid0OBll5cqVqK6uhkwmQ21tLbZs2SJ0SBlp8+bNuOSSS1BSUgKGYfCPf/xD6JAy0rJly3DqqadCrVajsLAQl19+Ofbu3St0WBlp1apVmDBhQrxw35QpU/Duu+8KFg8lN8Nk+vTpePXVV7F371688cYbOHjwIK688kqhw8o43333HaLRKJ566il88803eOyxx/Dkk0/id7/7ndChZaxgMIirrroKN910k9ChZJR169Zh4cKFWLx4MRoaGjBt2jTMmjULLS0tQoeWcbxeLyZOnIi//vWvQoeS0TZt2oRf//rX2LZtG+rr6xEOhzFz5kx4vV6hQ8s4ZWVl+MMf/oAdO3Zgx44dOPfcc3HZZZfhm2++ESQeWgp+jPzzn//E5ZdfjkAgQDvKDuFPf/oTVq1ahcbGRqFDyWhr167FwoUL4XA4hA4lI5x22mk45ZRTsGrVqvix8ePH4/LLL8eyZcsEjCyzMQyDt956C5dffrnQoWS8rq4uFBYWYtOmTTjrrLOEDifj6fV6/OlPf8K8efOO+WNTz80xYLPZ8NJLL2Hq1KmU2CTB6XRCr9cLHQYZQYLBIHbu3ImZM2f2OT5z5kxs3bpVoKhItnE6nQBAr09DiEQi+Pvf/w6v14spU6YIEgMlN8Pot7/9LZRKJfLz89HS0oK3335b6JAy3sGDB/GXv/wF8+fPFzoUMoJ0d3cjEonAaDT2OW40GmE2mwWKimQTnudRV1eHM888EyeeeKLQ4WSkr776CiqVClKpFPPnz8dbb72F448/XpBYKLlJwZIlS8AwTMKvHTt2xM+/88470dDQgA0bNoDjOMydOxe5MgqYalsBQEdHBy688EJcddVVuOGGGwSKXBhH0l6kP4Zh+tzmeb7fMUKOxM0334zdu3fjlVdeETqUjDV27Fjs2rUL27Ztw0033YRrrrkG3377rSCxiAR51BHq5ptvxs9+9rOE51RVVcX/bTAYYDAYcNxxx2H8+PEoLy/Htm3bBOumO5ZSbauOjg5Mnz4dU6ZMwerVq4c5usyTanuRvgwGAziO69dLY7FY+vXmEJKqW265Bf/85z+xefNmlJWVCR1OxpJIJBg9ejQAYPLkyfj888/x+OOP46mnnjrmsVByk4JYsnIkYj02gUAgnSFlrFTaqr29HdOnT0dtbS2ee+45sGzudSgezd8W6X1Rra2tRX19Pa644or48fr6elx22WUCRkZGMp7nccstt+Ctt97Cxo0bUV1dLXRIIwrP84K951FyMwy2b9+O7du348wzz0ReXh4aGxtx3333YdSoUTnRa5OKjo4OnHPOOaioqMAjjzyCrq6u+M+KiooEjCxztbS0wGazoaWlBZFIBLt27QIAjB49GiqVStjgBFRXV4err74akydPjvcAtrS00PytAXg8Hhw4cCB+u6mpCbt27YJer0dFRYWAkWWWX//613j55Zfx9ttvQ61Wx3sGtVot5HK5wNFllt/97neYNWsWysvL4Xa78fe//x0bN27Ee++9J0xAPEm73bt389OnT+f1ej0vlUr5qqoqfv78+XxbW5vQoWWc5557jgcw4BcZ2DXXXDNge3300UdChya4J554gq+srOQlEgl/yimn8Js2bRI6pIz00UcfDfg3dM011wgdWkYZ7LXpueeeEzq0jHP99dfHn3sFBQX8eeedx2/YsEGweKjODSGEEEKySu5NbiCEEEJIVqPkhhBCCCFZhZIbQgghhGQVSm4IIYQQklUouSGEEEJIVqHkhhBCCCFZhZIbQgghhGQVSm4IIYQQklUouSGEkEE0NzeDYZj4FhcbN24EwzBwOBzxc/7xj39g9OjR4DgOCxcuHPTYQD788EOMGzcO0WgUQO/u8JMmTRqeXwb9f59j6dprr8Xll18ev33llVdi+fLlxzwOkhsouSEkAYZhEn5de+21aXmMf/zjH0d9HTL8pk6dCpPJBK1WGz9244034sorr0RrayseeOCBQY8N5K677sLixYsF3yz20KFDkEqlcLlcx+wx77vvPjz00EPH9DFJ7qCNMwlJwGQyxf+9bt063Hfffdi7d2/8WC5snhcKhSAWi4UOI07IeCQSSZ8NXT0eDywWCy644AKUlJQMemwgW7duxf79+3HVVVcNe9xDefvtt3HOOedAo9Ecs8ecMGECqqqq8NJLL+Gmm246Zo9LcgP13BCSQFFRUfxLq9WCYZj4bbFYjPnz56OsrAwKhQInnXQSXnnllT73r6qqwooVK/ocmzRpEpYsWRL/OQBcccUVYBgmfhsAVq1ahVGjRkEikWDs2LF44YUX+lyHYRg89dRTuPjii6FQKDB+/Hh8+umnOHDgAM455xwolUpMmTIFBw8e7HO/d955B7W1tZDJZKipqcHSpUsRDof7XPfJJ5/EZZddBqVSiQcffHDAtgkGg7jrrrtQWloKpVKJ0047DRs3boz/fO3atdDpdHj//fcxfvx4qFQqXHjhhX0SRgB47rnnMH78eMhkMowbNw4rV66M/yw2jPLqq6/inHPOgUwmw4svvohwOIxbb70VOp0O+fn5+O1vf4trrrkmPuzx/PPPIz8/H4FAoM9j/eQnP8HcuXMH/H0AYPv27Tj55JMhk8kwefJkNDQ09Pn54cNSGzduhFqtBgCce+65YBhm0GMD+fvf/46ZM2dCJpMNGk9TUxNGjx6Nm266KT509cknn+Dss8+GQqFAXl4eLrjgAtjtdgDAe++9hzPPPDPeLhdffHG///+BvP3227j00ksB/Hf46OGHH4bRaIROp4v/jdx5553Q6/UoKyvDmjVr+lzjq6++wrnnngu5XI78/Hz87//+LzweT8LHvfTSS/s9ZwhJC8G27CRkhHnuued4rVYbv93W1sb/6U9/4hsaGviDBw/yf/7zn3mO4/ht27bFz6msrOQfe+yxPteZOHEi//vf/57neZ63WCzxXYZNJhNvsVh4nuf5N998kxeLxfwTTzzB7927l3/00Ud5juP4Dz/8MH4dAHxpaSm/bt06fu/evfzll1/OV1VV8eeeey7/3nvv8d9++y1/+umn8xdeeGH8Pu+99x6v0Wj4tWvX8gcPHuQ3bNjAV1VV8UuWLOlz3cLCQv7ZZ5/lDx48yDc3Nw/YHnPmzOGnTp3Kb968mT9w4AD/pz/9iZdKpfy+ffvi7SUWi/nzzz+f//zzz/mdO3fy48eP5+fMmRO/xurVq/ni4mL+jTfe4BsbG/k33niD1+v1/Nq1a3me5/mmpiYeAF9VVRU/p729nX/wwQd5vV7Pv/nmm/yePXv4+fPn8xqNhr/ssst4nud5n8/Ha7Va/tVXX40/VldXFy+RSPq04eE8Hg9fUFDAz549m//666/5d955h6+pqeEB8A0NDTzP/3c3bbvdzgcCAX7v3r08AP6NN97gTSbToMcGMnHiRP4Pf/hDn2O///3v+YkTJ/I8z/NfffUVX1xczN99993xnzc0NPBSqZS/6aab+F27dvFff/01/5e//IXv6urieZ7nX3/9df6NN/6/vbuPqSn+4wD+7laTymPsKpKw7KaH3QhpHibc1TzN0yTKJEZr5HGtMkZkeShpSbgtQ/74sVCopT+Qh1akdPVAk7UlqRU90b2f3x/tnl+ne0vKj2mf19bW93O+3+/5fs+52/l0zvee/kOlpaX08uVLWrp0KTk6OpJarRYdT+18iIjq6+vJ2NiYKisriajjv84PGTKEAgMD6e3bt3Tp0iUCQAqFgiIiIqi0tJSOHDkiatPU1ERWVla0cuVKKiwspKysLLK1tRX9l3E/Pz/h/Gilp6fToEGDqLW1Ve8xYqyvOLlhrJe6Jjf6eHl50Z49e4Tyz5Iboo5k4tatW6I6s2fPpoCAAFFszZo15OXlJWoXFhYmlJ8+fUoA6NKlS0Ls+vXrZGJiIpTnzJlDx44dE/V75coVsrS0FPW7a9euHudZXl5OBgYGVFVVJYp7eHhQSEgIEXUcLwBUXl4ubI+LiyOpVCqUra2t6dq1a6I+jhw5Qm5ubkT0v4txdHS0qI5UKqWoqCih3N7eTuPHjxddPLdv306enp5COTo6miZOnEgajUbvnBISEmjkyJHU1NQkxOLj47tNbog6EgMAlJ2dLbTRF9Nn2LBhlJycLIppk5ucnBwaOXKkaI5ERN7e3uTu7t5jv51pk+fCwkIi0p/cXL16lVxcXISyn58f2djYCAkREdGUKVNozpw5Qrm9vZ3MzMzo+vXrRNSRpI4YMYK+ffsm1ElLSyOJRELV1dVCv12Tm4KCAgLQbQLNWF/xmhvG+kitViMyMhI3btxAVVUV2tra0NbWBjMzs373rVKpsHXrVlHM3d0dMTExopiTk5Pwu1QqBQA4OjqKYq2trWhsbMTQoUORl5eH3NxcREREiObR2tqK5uZmmJqaAgCmT5/e4/jy8/NBRLCzsxPF29raYGFhIZRNTU0xadIkoWxpaYmamhoAwOfPn/Hx40f4+/sjICBAqNPe3i5asNt1PA0NDfj06RNmzJghxAwNDTFt2jTh0Q0ABAQEwNXVFVVVVRg7diyUSiU2bdoEAwMDvXNSqVRwdnYWjgEAuLm59Xgc+qOlpUXvI6nKykosXLgQR48eRXBwsGjbq1evelyj8+7dO4SHh+PZs2eora0VjkdlZSUcHBz0tun8SEpr6tSpokXOUqlU1N7Q0BAWFhbCudQeu86ffXd3d2g0GpSUlAifza60a9aam5u7nRNjfcHJDWN9dOrUKZw5cwbR0dFwdHSEmZkZdu3ahe/fvwt1JBIJiEjU7sePH73qv+tFmIh0Yp0X1mq36YtpL3IajQaHDx/GypUrdfbX+UL7swRNo9HA0NAQeXl5MDQ0FG0zNzfXOz7teLTHQzumxMREzJw5U1Sva5/6xqPv+HQml8vh7OyM5ORkKBQKFBYW4s6dO93OqWv7/7dRo0YJa2U6Gz16NKysrJCSkgJ/f3/RIt+fLWBfunQprK2tkZiYCCsrK2g0Gjg4OIg+k539+PED9+/fR0hIiCiu77zpi2nPob7PZud63amrqwPQMWfGfideUMxYHz169AjLly/Hhg0b4OzsjIkTJ6KsrExUZ/To0aIFtI2NjaioqBDVMTY2hlqtFsVkMhkeP34siuXk5EAmk/VrzC4uLigpKcHkyZN1fn7l68hyuRxqtRo1NTU6/XT+NlFPpFIpxo4di/fv3+v0YWtr2227YcOGQSqV4sWLF0JMrVbrLP4FgC1btkCpVOLy5ctYuHAhrK2tu+3X3t4eBQUFaGlpEWLPnj3r1Vz6Qi6Xo7i4WCc+ePBg3L17FyYmJlAoFPj69auwzcnJCVlZWXr7+/LlC1QqFcLCwuDh4QGZTKY3eeosOzsbw4cP7/e7dezt7fHq1Ss0NTUJsSdPnkAikejc3eusqKgI48aNw6hRo/q1f8a64uSGsT6aPHkyMjMzkZOTA5VKhW3btqG6ulpUZ8GCBbhy5QoePXqEoqIi+Pn56dyVmDBhArKyslBdXS1cjPbt24ekpCScP38eZWVlOH36NG7evIm9e/f2a8wHDx5EcnIyDh06hDdv3kClUuHGjRsICwv7pX7s7Ozg4+MDX19f3Lx5ExUVFcjNzcWJEyeQnp7e634OHTqE48ePIyYmBqWlpSgsLIRSqfzpy92CgoJw/PhxpKamoqSkBDt37kR9fb3OXQIfHx9UVVUhMTERmzdv7rHP9evXQyKRwN/fH8XFxUhPT8fJkyd7PZdfpVAodBJYLTMzM6SlpcHIyAienp7Ct45CQkKQm5uLHTt24PXr13j79i3i4+NRW1uLESNGwMLCAhcuXEB5eTkePnyI3bt39ziG27dv6zyS6gsfHx+YmJjAz88PRUVFyM7ORlBQEDZu3NjtIymg4w+ExYsX93v/jHXFyQ1jfRQeHg4XFxcoFArMnz8fY8aMEb2BFei4GM2dOxdLliyBl5cXVqxYIVqDAnQ83srMzIS1tTXkcjkAYMWKFYiJiUFUVBSmTp2KhIQEKJVKzJ8/v19jVigUuHv3LjIzM+Hq6opZs2bh9OnTsLGx+eW+lEolfH19sWfPHkyZMgXLli3D8+fPe7w70tWWLVtw8eJFJCUlwdHREfPmzUNSUlKPd24A4MCBA/D29oavry/c3Nxgbm4OhUKhs4Zl6NChWLVqFczNzXXOTVfm5ua4c+cOiouLIZfLERoaihMnTvR6Lr9qw4YNKC4uFr03qet47t27ByKCl5cXmpqaYGdnh4yMDBQUFGDGjBlwc3NDamoqjIyMIJFIkJKSgry8PDg4OCA4OBhRUVE9juH27dtYvnx5v+diamqKBw8eoK6uDq6urli9ejU8PDxw7ty5btu0trbi1q1bovVWjP0uBvSnHzQzxthvptFoIJPJsHbtWp03Ai9atAgymQxnz579S6Pr3v79+9HQ0ICEhIQ/vu/8/HwsWLAAnz9//isvRYyLi0NqaioyMjL++L7ZwMd3bhhj/5wPHz4gMTFReJS1fft2VFRUYP369UKduro6pKSk4OHDhwgMDPyLo+1eaGgobGxsdNZc/Qnt7e2IjY39a297NjY2Rmxs7F/ZNxv4+M4NY+yf8/HjR6xbtw5FRUUgIjg4OCAyMhJz584V6kyYMAH19fUIDw/v91olxti/hZMbxhhjjA0o/FiKMcYYYwMKJzeMMcYYG1A4uWGMMcbYgMLJDWOMMcYGFE5uGGOMMTagcHLDGGOMsQGFkxvGGGOMDSic3DDGGGNsQPkvb0vquhMrLEwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x, weights[:,1])\n",
    "plt.xlabel('Tautomer energy diff (kcal/mol)')\n",
    "plt.ylabel('Fraction existing')\n",
    "plt.fill_between(x, \n",
    "                 weights_down[:,1], \n",
    "                 weights_up[:,1], \n",
    "                 color='grey', \n",
    "                 alpha=0.3)\n",
    "plt.grid()\n",
    "plt.yticks(np.linspace(0,1, 11));"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d90d970-32f6-422d-a550-ef05ae759be5",
   "metadata": {},
   "source": [
    "# weighting binding constants. \n",
    "\n",
    "Imagine you have a compound with just two species - a major and minor tautomer. The major is 90% populated. The minor is 10% populated. Each of these can bind a pocket, but with a different dissociation constant due to the change in HBD/HBA properties. \n",
    "\n",
    "To calculate the combined dissociation constant, we convert to equilibrium binding constants, sum the weighted contributions of each tautomer, and convert back to dissociation constant. Note that the measured Ki (pKi=7) is less than the best Ki (pKi=8)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ef8dfdb8-a689-4212-8590-594c10080588",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Combined Ki: 9.910802775024777e-08\n",
      "Combined pKi: 7.003891166236911\n"
     ]
    }
   ],
   "source": [
    "major_ki = 10 #µM\n",
    "minor_ki = 0.01\n",
    "\n",
    "major_weight = 0.9\n",
    "minor_weight = 0.1\n",
    "\n",
    "#convert to keq:\n",
    "major_keq = 1/(major_ki/1e6) #molar\n",
    "minor_keq = 1/(minor_ki/1e6) #molar\n",
    "#combine:\n",
    "combined_keq = major_weight*major_keq + minor_weight*minor_keq\n",
    "#combined ki:\n",
    "combined_ki = 1/combined_keq\n",
    "print(f'Combined Ki: {combined_ki}')\n",
    "print(f'Combined pKi: {-np.log10(combined_ki)}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49d5e890-f98c-46a2-a143-aa665f0aaf54",
   "metadata": {},
   "source": [
    "# now, use _estimated_ tautomer populations\n",
    "\n",
    "imagine your tautomer populations are _estimated_ tautomer populations. Our estimator gives a 95% confidence interval of exactly 1kcal/mol. \n",
    "\n",
    "We can model the range of expected affinities given the uncertainty in the tautomer populations. This model will assume that we have a perfect estimate of the dissociation constant for each species (thus removing that source of error from our predictions). In practice, of course, our model of dissociation constants is imperfect as well. \n",
    "\n",
    "The 90%/10% populations estimate corresponds to an estimated energy difference of +1.2963625 kcal/mol for the minor tautomer with respect to the major tautomer. Assuming an evenly distributed confidence interval, the true energy difference is somewhere between 1.2963625-0.5 and 1.2963625+0.5, i.e. +/- 0.5 kcal/mol, which would be:\n",
    "\n",
    "- 1.2963625+0.5 (aka upper estimate): 95% / 5%\n",
    "- 1.2963625-0.5 (aka lower estimate): 79% / 21%\n",
    "\n",
    "Assuming we know the exact dissociation constant for each species, upon measuring experimentally, we would expect the _combined_ dissociation constant to be in the range:\n",
    "\n",
    "**6.7 < True pKd < 7.3**\n",
    "  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "db4ebb7c-6a1e-493d-816b-e81cd99a7145",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "populations upper estimate: [0.95455239 0.04544761]\n",
      "populations lower estimate: [0.79409152 0.20590848]\n"
     ]
    }
   ],
   "source": [
    "p = -np.log([0.9, 0.1])*RT\n",
    "p -=p.min()\n",
    "p[1]+=0.5\n",
    "weights = np.exp(-1*p/(RT))\n",
    "weights = weights/weights.sum()\n",
    "print('populations upper estimate:', weights)\n",
    "p = -np.log([0.9, 0.1])*RT\n",
    "p -=p.min()\n",
    "p[1]-=0.5\n",
    "weights = np.exp(-1*p/(RT))\n",
    "weights = weights/weights.sum()\n",
    "print('populations lower estimate:', weights)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d2d79d60-0336-43fb-8265-bee2868e9d6b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Combined Ki (upper): 2.1550719804719858e-07\n",
      "Combined pKi (upper): 6.666538219610108\n",
      "\n",
      "Combined Ki (lower): 4.837869179112959e-08\n",
      "Combined pKi (lower): 7.315345879578421\n"
     ]
    }
   ],
   "source": [
    "major_ki = 10 #µM\n",
    "minor_ki = 0.01\n",
    "\n",
    "\n",
    "major_weight = 0.95455239\n",
    "minor_weight = 0.04544761\n",
    "\n",
    "#convert to keq:\n",
    "major_keq = 1/(major_ki/1e6) #molar\n",
    "minor_keq = 1/(minor_ki/1e6) #molar\n",
    "#combine:\n",
    "combined_keq = major_weight*major_keq + minor_weight*minor_keq\n",
    "#combined ki:\n",
    "combined_ki = 1/combined_keq\n",
    "print(f'Combined Ki (upper): {combined_ki}')\n",
    "print(f'Combined pKi (upper): {-np.log10(combined_ki)}')\n",
    "\n",
    "print()\n",
    "\n",
    "major_weight = 0.79409152\n",
    "minor_weight = 0.20590848\n",
    "\n",
    "#convert to keq:\n",
    "major_keq = 1/(major_ki/1e6) #molar\n",
    "minor_keq = 1/(minor_ki/1e6) #molar\n",
    "#combine:\n",
    "combined_keq = major_weight*major_keq + minor_weight*minor_keq\n",
    "#combined ki:\n",
    "combined_ki = 1/combined_keq\n",
    "print(f'Combined Ki (lower): {combined_ki}')\n",
    "print(f'Combined pKi (lower): {-np.log10(combined_ki)}')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}

relying on this as the method of combining binding constants. It's from https://pubs.acs.org/doi/10.1021/acs.jctc.3c00870 which, in turn, references Gallicchio, equation 22, https://pubs.acs.org/doi/10.1021/ct1002913 . Note that K_i is not inhibition constant, but the i-th Keq