dwhswenson/contact_map

View on GitHub
examples/concurrences.ipynb

Summary

Maintainability
Test Coverage
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Concurrences\n",
    "\n",
    "One of the tools in Contact Map Explorer is the ability to look at simultaneous contacts. The idea is that you might have a set of contacts that is likely to happen concurrently, and that this set of contacts might help you define a stable state. This is managed in Contact Map Explorer by what we call contact concurrences.\n",
    "\n",
    "To start, we'll look at a trajectory of a specific inhibitor during its binding process to GSK3B."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<mdtraj.Trajectory with 100 frames, 5704 atoms, 360 residues, and unitcells>\n"
     ]
    }
   ],
   "source": [
    "from __future__ import print_function\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from contact_map import ContactFrequency, ResidueContactConcurrence, plot_concurrence\n",
    "import mdtraj as md\n",
    "traj = md.load(\"data/gsk3b_example.h5\")\n",
    "print(traj)  # to see number of frames; size of system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we'll use MDTraj's [atom selection language](http://mdtraj.org/latest/atom_selection.html) to split out the protein and the ligand, which has residue name YYG in the input files. We're only interested in contacts between the protein and the ligand (not contacts within the protein). We'll also only look at heavy atom contacts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "topology = traj.topology\n",
    "yyg = topology.select('resname YYG and element != \"H\"')\n",
    "protein = topology.select('protein and element != \"H\"')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll make a list of all the residue contacts that are made, keeping only those that occur more than 20% of the time. We'll put that into a `ResidueContactConcurrence` object, and plot it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 25.8 s, sys: 37.2 ms, total: 25.8 s\n",
      "Wall time: 2.51 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "contacts = ContactFrequency(traj, query=yyg, haystack=protein)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "contact_list = [(contact_pair, freq)\n",
    "                for contact_pair, freq in contacts.residue_contacts.most_common()\n",
    "                if freq >= 0.2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 361 ms, sys: 4.04 ms, total: 365 ms\n",
      "Wall time: 365 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "concurrence = ResidueContactConcurrence(traj, contact_list, select=\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# optionally, create x-values... since we know that we have 1 ns/frame\n",
    "times = [1.0*(i+1) for i in range(len(traj))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfMAAAEGCAYAAABxUHzhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABwhElEQVR4nO3dd1xTV/sA8OcmzMjeEAggIWQwBJx11NFqHVhX1YLY2lerfbW1jrpatVq17rdqa92taLUoWhVfZ621P6tFcTFCiMgeEWRDAmTc3x94eSOSsBIh+Hw/Hz7mnnvuOedeIk/uzb3nIUiSBIQQQggZLlpHDwAhhBBC7YPBHCGEEDJwGMwRQgghA4fBHCGEEDJwGMwRQgghA2fUmsoODg6kl5eXnoaCEEJd0717956RJOnYju2djIyMDgCAP+BJ2OtIBQBJCoViZmhoaGFTFVoVzL28vCA+Pl4nI0MIodcFQRBZ7dneyMjogIuLC8/R0bGURqPh88SvGZVKRRQVFfElEskBABjbVB38hIcQQp2fv6OjYwUG8tcTjUYjHR0dy6H+ykzTdV7heBBCCLUNDQP56+35719jzMZgjhBCCBk4DOYIIYSQgcNgjhBCSKvU1FQTMzOzEC6Xy6fK6HR6KJfL5fv6+gpGjhzZvbKykgYAwGAwgtW33blzp/306dNZAAALFy50c3JyCuRyuXzq59mzZ/Tr168zqGU/Pz9+VFSUDbX9wIEDff38/PhsNlsQHh7OUigUGseZl5dnxGQyA7Kzsxtu7p42bRpr8eLFrk2Vr1ixwgUA4Pr164zevXv7eXp6+vP5fN7gwYPZd+7cMQcA2Lx5syOHw+FzuVx+aGio371798waHwMul8sfOnQomyrfsGGDI4vF8icIIrSgoKChz/3799uyWCz/IUOGNNTVFQzmCCGEmuXh4VErEomE1LKpqalKJBIJHz9+nGxsbExu27atRY/ezZkz56lIJBJSPw4ODsqePXvWJCYmCkUikfDKlSuPP//8c0+5XA4AAGfPnn2SmpoqFIvFycXFxcaHDh2y1dQ2k8lUzJ8/X/Lpp596AADcvHmTcefOHYv169dLmipfvXr105ycHKNp06b5rF+/PjcrKytJKBSmLF++XJKammoKADBz5sxisVgsFIlEwoULF0o+//xzj8bHQCQSCf/44480qvzNN9+sunr1qtjNza1OfXyzZs0q3b17d7uebNAEgzlCCHVBKpKEe1ml5qpXkBlzwIABVWlpaaZt3d7S0lJlbGwMAAAymYwgCKJhnZ2dnQoAQC6XE3K5/IV1TVm0aFFRVlaWaWxsrOWnn37K2rFjR7apqSmpqXzr1q1OkydPLn777berqTZGjBhRFRkZWabePwBAVVUVvbn+AQD69+8v8/Pzq2u2og5hMEcIoS7oQXaZ+dxj930eZJeZ67MfuVwOly9ftgoICJABANTW1tLUL6N/++23bur19+zZ40yt69OnD4cq/+OPP7qx2WxBSEiI4D//+U8WFdwBAAYMGODr6OgY1K1bN+WMGTNKtY2HTqfD7t27s6ZNm+bTvXv3mpEjR1ZpK09JSTEPDQ2Vamvz22+/dfTw8PBfvXq1+w8//JBNldfV1dH8/f15QUFB3CNHjti0/KjpHgZzhBDqgoJZNrIfwkOeBLNsZPponwraAQEBfHd397r58+c/A3jx0rNIJBIuX748X3079cvscXFxYqp86NCh1Wlpack3b95M2bJli6tUKm04Bb558+ZjiUTyqK6ujhYbG2vV3NjeeOMNma+vr2zevHmFLSlXFxgYyO3evbtgxowZDZfTly9fXpSTk5P09ddf565evdqVKk9LS0tISkpKOX78ePqyZcs8kpOT23x1or0wmCOEUBdEIwgI9bSV0VpwWbgt1IP24cOHc8zMzHRyPT8kJKSGwWAo4+PjX7iiwGAwyDFjxpT99ttvNi1ph0ajAZ1Ob7acx+PJ7t27x6CWExISRCtXrsyvqKh4aeNZs2aVXL16taF/Ly8vOQAAn8+v69u3b+WdO3cYjbd5VXQezEmShOT8ciBb8T1NW7bRh8bjUF/W9FoX9bSNozVtvMrxdpY2OoPOeGy6al9tbaMtvy9djpcwNtPrpW5DJxKJTKgb3sRisUlGRoaZr69vXXl5OS0rK8sYoP5y/qVLl6y5XK4MACAqKspm7ty5zPb2vWjRosLo6Gj7q1evdqPKqqurG2JjYmJiw9l2dHS0taenZy0AQFFREV0mkxEAAAUFBUbx8fEWgYGBerkK0hI6D+bCggr45Oh9EBZU6HUbfWg8DvVlTa91UU/bOFrTxqscb2dpozPojMemq/bV1jba8vvS5XiNbFx8dPR2M3jq35lzuVx+amqqybVr1yx4PJ6Ay+Xyx40b57Nt27ZsV1dXRUVFBW306NFsDofD5/F4AgcHB/kXX3xRBACQlpZmamVlpWzveFgsluLIkSPpK1ascGexWP7BwcHc06dP286fP78QAGD79u1ObDZbwOVy+Tt27HD++eefMwAAHj58aBYUFMTz8/Pjv/nmm5zPP/9cEhoaWgMAsG7dOidnZ+fAp0+fmgQFBfGnTJni2d5xNodozRlOz549yeYSrZAkCcKCCuC7WkFL7vpr6zb60Hgc6ssA0ORrXdRrvM/a6rVlnT7G21na6Ax0/b4xtN/lq+yrrW2ov1da+vvS5Xjppgyhqk4maOt77NGjR5lBQUHP2rp9e6WmppqMGTPG9/Hjx8kdNYbG3n33Xe8ff/wxx83NTfOD553Q+fPnLbdt2+Z8/fr1tOZrv+jRo0cOQUFBXk2t03kwRwgh9CKCIO6RJNmzrdt3dDBPS0szHjBgAM/Gxkah/qw5ap39+/fbbty40S0gIEB65syZjNZury2YtyoFKkIIodcPm82WSySShI4eh6GbNWtW6axZs7Q+WtdWeDc7QgghZOAwmCOEEEIGDoM5QgghZOAwmCOEEEIGDoM5QgghrZpKgUqlOk1NTTXx9fV96bG7iRMnejGZzADqefLg4GAuAMCPP/5ox+Fw+BwOhx8cHMy9fft2w4Q6a9ascWKz2QJfX19BWFiYt/qUro3Fx8ebeXl5+VdVVTXUGTx4MHvp0qUuTZXv27fPFgAgJibGKiAggOft7S3gcrn80aNHd3/8+LEJAMD8+fPdqHSn/fv3983MzDRuvP9cLpcfHh7Ootr+9NNPmS4uLoGNU7+uWbPGydXVNYBK/6pvGMwRQgg1q3EK1JZYt25dLjXl64MHD0QAAGw2u/bvv/9OFYvFwuXLl+fPnj3bEwAgIyPDeN++fc4PHz4UPn78OFmpVBIHDhyw09R2z549a0aNGlW6YsUKVwCAI0eO2MjlcmLTpk2Spso//vjj0rt375otWrSIdfjw4YyMjIxkkUgkDA8PL05LSzMBAFi9erWESnc6cuTIcqoN9f0XiUTCY8eONSRbGTduXFlcXFxK4/GtXr26sPG89PqEj6YhhFBXRKoAcu+ag3svGRCd57xNPdXokCFDqufNm2dCLSuVSqK6uppmamqqlMlkNHd3d7m2tjZt2lQQEBDAnzp1aumqVauY586dS9NWvn79eteFCxcWhISE1FBtRERElFOv1dOdVldX01oyMdWwYcOqm630CnSe3zBCCCHdyb1rDic/9IHcux02L/xXX33lTl2aHjt2rHfj9bt27XIYMmRIOQCAt7e3fO7cuRJvb+9AJyenIEtLS+WECRO0zttsaWmp2rBhQ87w4cO548ePLw0ICKjVVi4Wi8169+6tNd0pddk8JibGfsuWLQ1n1rm5uSY8Ho/fq1cvv0uXLlm05Xjok36DOUkCFCTU/9t4WdNrbfV00UZX7QvHq/82EDIk7r1k8N7PT8C9V4cl/1C/zH7u3LkXZjyLjY21PHr0qMOOHTtyAeoTl/z3v/+1SUtLS5RIJAlSqZS2e/dujZfZKeHh4eWWlpaKRYsWFbaknCKRSOhcLpfv5eXlv2rVKmeqfNeuXXkSiSRh0qRJxVu2bHECAGCxWPKMjIyElJQU4fbt23M+/PDD7iUlJZ3qZFi/g5EkApyIrP+38bKm19rq6aKNrtoXjlf/bSBkSAgagEefTnWJnRIXF2f+73//2/PMmTNpLi4uSgCA2NhYKxaLVevm5qYwNTUlx40bV3br1q0WnQHTaDSg0V7ez8blHA6nhkpT6uLiohSJRMLp06cXVVVVvZTudMaMGSXnz5+3BQAwNzcnqXEOHDhQymKxapOSkszasu96Q6Xsa8lPaGgo2SoqFUnmP6r/t/Gyptfa6umija7aF45X/20g1EYAEE+24m9t45+HDx9mkiQZ31E/IpEogc1my9TLzM3NlZrWkSQZP2HChGeHDh160rhcLBYneHh41Fy5ciVFvfzatWspPj4+soqKivtKpTJ+/Pjxz9atW5dNkmT8+vXrs9avX5+laXxubm61+fn5D5srj4uLS/bw8Ki5d+9eElW2aNGivAULFuSTJBmfkJCQSJWvW7cue8SIESUkScbn5eU9lMvl8SRJxicnJyc4OjrWSSSSB00dD/WfHTt2ZERGRhbq6vfw/H3Q5HsEE60ghJCeGXqilaaypjEYjGCpVPogNTXVRCAQ+Nvb2zdkL/v2229z/vvf/1r/888/lpaWlg1pSh8+fJjywQcfsC5evGjr5uZWBwBgZGREJiUlpQAALFiwwO3MmTO2RkZGIBAIpMePH880Nzcnp0+fzurfv3/V7NmzS5oaH5PJDIiPj09xdXVVNFf+66+/Wq9bt86turqaZmtrq2QymbXr16/PDwwMrB0xYoRPenq6GUEQpLu7e93BgwezvL295T///LPNunXrmHQ6naTT6eRXX32VHx4eXg4AMGfOHPfffvvNrqioyNjR0VEeERHxbPv27fkAADt37rSPj4/vFhUVlQ06gFnTEEKoA3XFYP4qDRkyhH3x4sUnZmZmBnUDy6sM5p3vyxSEEEKdCp1OJysrK+nqk8a8StevX08ztEC+Zs0ap+3bt7taWVkpm6/dfvicOUIIIa0wBWrrrV69unD16tVN3kmvD3hmjhBCCBk4DOYIIYSQgcNgjhBCCBk4DOYIIYSQgcNgjhBCSKumUqDm5OQYhYWFebu7uwcIBAJejx49uFFRUTYAAOfPn7ccMmQIGwBAKpUS3t7egjt37jTMEf/VV185h4eHs27dumXeo0cPLpvNFnA4HP7+/fttqTrnzp2z5PP5PF9fX8GECRO85HLNOVe09aGpHAAgMTHRdMiQIWwPDw9/gUDA69OnD+fixYsWAABHjx61odKh+vv78y5fvtwwGx2TyQxQX0eVHzp0yJbNZgtoNFroX3/9xaDKL126ZOHj4yNoKlWsrmAwRwgh1Cz1FKgqlQrCwsLYAwcOrMrNzU1MTk5OOXHiRHpOTo5J4+0YDAa5ZcuWnE8++YSlUqkgIyPD+PDhw447duzIs7CwUB05ciQjLS0t+cqVK49XrFjh8ezZM7pSqYSPP/7Y+9dff01//PhxMovFqvv+++8dNI1NWx+ayqVSKREWFuY7c+bMopycnKTk5OSU77//Pvvx48emAABhYWEV1LzyBw8ezJwzZ46nep83btwQi0QiITXhDQBAjx49ZKdOnUrr2bNnlXrdd955p+rChQuP2/s70EavwZwkSRCViICamEZ9WdNrbfV00UZX7QvHq9/3nra+2vr/ASF9UpEqeFj40FxFqpqv3EqxsbGWxsbG5JIlS4qoMg6HU/fll182+SjWpEmTKpydneU//PCD/dy5cz2WLVuW7+joqAwMDKylMpp5eXnJ7ezsFAUFBUZPnz41MjExUQUGBtYCALzzzjsVZ86csdE2Jk19aCrfu3evfUhISJV6CtRevXrVfPbZZ8UAANbW1ipqXvfKysoWpUMNCQmpCQoKqm22oh7oNZinlqbCgusLILU09aVlTa+11dNFG121Lxyvft972vpq6/8HhPQpoSjBfPGNxT4JRQk6T4GamJhoHhgYqDWVaGO7d+/OWbduHbO4uNho7ty5L03Lev36dYZcLif4fH6ti4uLQqFQENSl6ujoaNuCgoKXzvpb2kdT5cnJyWbBwcFa9yEqKsrG29tbMHHiRN99+/Zlqq8bNmyYr0Ag4G3dulXjFYNXStOk7U39tDbRikqlIlOKU0jV80QV6suaXmurp4s2umpfOF79vve09dXW/w/o9QEdkGhFqVLGP3j6IFmpUuo80co333yT/dFHHz2llqdNm1bI4XCkAoGgmiTJ+NjY2NTBgweXNW5n3Lhxz/bv3/9S8pXMzMxHnp6eNb///ntD8pWrV6+mhISEVPr7+1d/+umn+VwuV9qSsWrqo3H5v/71L8natWuzqeW33nqrlM1my95+++3SxtteuHBB1K9fvwpqOSMj4xFJkvG5ubkPORyO9MKFCyL1+r169aq8ceOGUNsxbMuPtkQrej0zJwgCuHZcoC5PqC9req2tni7a6Kp94Xj1+97T1ldb/z8gpE80ggY9nHrIaHpIgRoQECBLSEhouMHryJEj2X/++ae4tLRU66yiNBoN6PQXs42WlJTQRo4cyV61alXesGHDqqnyt956q/revXupiYmJKYMHD67y9vauacnYmuqjqXKBQFDz4MGDhn24evXqk4MHD2aUlZW9tA8jR46sysrKMi0oKDACqP9KAACAyWQqRo8eXXb79u1uLRmbPuENcAghhFolLCyssra2lti0aZMjVVZVVdXqeFJTU0OMHj2aPXXq1OKPPvqoVH1dXl6eEQCATCYjtmzZ4jJnzpwigPrL8ePHj/dq5y7ArFmziuPj4y1++eUXa6qsurq6YR+SkpJMVar6+w1u3rzJkMvlhLOzs6KiooJWWlpKAwCoqKigXb9+3SowMFDW3vG0F87NjhBCqFVoNBrExsY+mTt3rsfOnTtd7OzsFAwGQ/n111/nUnVu375t5ezsHEgt//LLL08at3Po0CHbu3fvWpSWlhodO3bM4XlZxhtvvCFbu3aty9WrV61VKhXx0UcfFY4dO7YSACAzM9PU3Ny83XeRWlhYkGfPnk37/PPP3ZcuXcpycHCQd+vWTblixYp8AIDjx4/bRkdH2xsZGZFmZmaqI0eOpNNoNMjNzTUaP348GwBAqVQSEydOLJ40aVIFQP137F988QWrtLTUaPz48b48Hk968+ZNvd7FTsEUqAghpGeYAlV3Zs+e7f7RRx8V9+nTp8PPhltDF8cQU6AihBBqs45Ogapu7969uYYWyC9dumQRFhbGtrW1VeirD7zMjhBCSCtMgdo+77zzTpVYLBbqsw88M0cIIYQMHAZzhBBCyMBhMEcIIYQMHAZzhBBCyMDpPdFKTUrKC0ksqGVNr7XV00Ub2uppG7+2Nlqyja72S5/HVxfj1cf7QRfHt71ttOZ9g1BX0zgFau/evf1OnTplpV5n7dq1TtOmTWMBAOTn5xsZGRmFbNmy5YV5y5lMZgA1ixrlxx9/tONwOHwOh8MPDg7m3r592xwA4NGjR6ZcLpdP/VhYWASvXbvWSdMY4+Pjzby8vPyrqqoaplgcPHgwe+nSpS5Nle/bt88WACAmJsYqICCA5+3tLeByufzRo0d3f/z4sQkAwPz5892oVKf9+/f3zczMNG58PLhcLp9KqQoA8OmnnzJdXFwCGQxGsPr41qxZ4+Tq6howffp0FuiBXoN5rUgEuZ9+BrUi0UvLml5rq6eLNrTV0zZ+bW20ZBtd7Zc+j68uxquP94Mujm9722jN+wahrkg9Bep7771XfPz4cTv19adOnbKbNm1aCQBAVFSUbVBQUPXJkyftm2uXzWbX/v3336lisVi4fPny/NmzZ3sCAAQFBdVSKUiTkpKEZmZmqqlTp5Zpaqdnz541o0aNKl2xYoUrAMCRI0ds5HI5sWnTJklT5R9//HHp3bt3zRYtWsQ6fPhwRkZGRrJIJBKGh4cXp6WlmQAArF69WiIWi4UikUg4cuTIcqoN9eMhEomEx44dy6bKx40bVxYXF5fSeHyrV68uXL58eX5zx6PNNE3a3tRPWxKtyITCF5JYUMuaXmurp4s2tNXTNn5tbbRkG13tlz6Pry7Gq4/3gy6Ob3vbaM37BqHGoAMSraiUyvjq+/eTVUrdJ1opKCh4YGNjI5dKpfeo9S4uLrXK532FhIRU/vHHH0IPD4+a9PT0R9R2bm5utfn5+Q819VNYWPjA0dGxrnH5qVOnxMHBwVXNjbOiouK+p6dnzd9//53s5eUlS0hISNRW/u677xZ/9913GS05BsuWLcsNDw8vbOp4NPVjbm6ubFy2Y8eOjMjIyMK2/h46NNGKGY/3QhILalnTa231dNGGtnraxq+tjZZso6v90ufx1cV49fF+0MXxbW8brXnfINQZyB49Ms/7fIGP7NEjnadAdXFxUQYFBVWfOnXKGgDg8OHDdmPHji2l0WiQlpZm/OzZM+MhQ4ZIx44dW3r48GG75tqj7Nq1y2HIkCHljcuPHz9uN2nSpOLmtre0tFRt2LAhZ/jw4dzx48eXUrnSNZWLxWKz3r17a02DSl02j4mJsd+yZUvDmXVubq4Jj8fj9+rVy+/SpUsWLd1HfcEb4BBCqAsyDwqSMb/7zxPzoCC9zJY2efLkkujoaFsAgNOnT9tFRkaWAPwvsAMAREZGlsTExLQomMfGxloePXrUYceOHbnq5TU1NcTvv/9uHRkZWappW3Xh4eHllpaWikWLFhW2pJwikUjoXC6X7+Xl5b9q1SpnqnzXrl15EokkYdKkScVbtmxxAgBgsVjyjIyMhJSUFOH27dtzPvzww+4lJSUdGk8xmCOEUBdE0GjACA6WETT9/JmPiIgo+/vvv61u3rzJqKmpoQ0YMEAKUP/deXR0tD2TyQyYMGECOzU11TwxMdFUW1txcXHm//73vz3PnDmT5uLiolRfFxMTY83n86UeHh4tngqVRqMBrYn9blzO4XBq7ty5wwCov9ogEomE06dPL6qqqnoph+qMGTNKzp8/bwsAYG5uTlLjHDhwoJTFYtUmJSWZtXR8+oDBHCGEUKtZW1ur+vbtWzlz5kyvCRMmlADU34EulUrphYWFCXl5eYl5eXmJ8+bNk0RFRWk8O3/8+LHJe++953Po0KGMwMDA2sbrf/31V7vJkyeXqJdt2LDBccOGDY6N67bWihUrJNu2bXO9f/9+QyCWSqUNcVH9Q8jJkydtfHx8ZAD1d+srFPWfLYRCoUlmZqapn5/fS2N/lXBudoQQQm0yderUkg8++MDn+PHj6QAAhw8fth81alRpozql4eHh3bds2VIAABAUFMSn7jUJCwsrqayspJeVlRl9+umnngAARkZGZFJSUgoAQGVlJe3mzZtWhw8fzlJvUyQSmffv37+qvePv3bu3bPPmzTnTp0/3rq6uptna2iqZTGbt+vXr8wEAFi9e7J6enm5GEATp7u5ed/DgwSwAgCtXrlisW7eOSafTSTqdTn733XdZzs7OSgCAOXPmuP/22292NTU1NGdn58CIiIhn27dv199d7M9hClSEENIzTIGqW0OGDGFfvHjxiZmZmUFN9LBz5077+Pj4blFRUdnN134ZpkBFCCHUZp0pBSoAwPXr19MMLZCvWbPGafv27a5WVlbK5mu3Hl5mRwghpBWmQG2/1atXF65evbrJO+l1Ac/MEUIIIQOHwRwhhBAycB0WzEmShKKcyoap6KjX2uo1tdySetr6ass6XY9JWxutOW7t3Wd9H5tXOV5dvd8QQsgQdFgwf5ZbBZf2JsKz3KoXXmur19RyS+pp66st63Q9Jm1ttOa4tXef9X1sXuV4dfV+Qwghg6Bp0vamflqbaEUblUpFFmZXNCSxoF5rq9fUckvqaeurLet0PSZtbbTmuLV3n/V9bF7leHX1fkNIF6ADEq3o8kckEiWYmpqq/Pz8pOrlhw8fTgMA8v79+0lN1e3evbts3Lhxz2pqau5R6//44w9hr169KlksVg2Px6t+8803y+Li4pJJkoy/cOGCiMfjVdPpdPLQoUNP1PvauXNnBovFqmGxWDU7d+7MoMpDQkIq/fz8pH5+flJHR8e6YcOGlWraj+rq6nteXl4yqj+SJOO//PLLnPfff79QUzlJkvEJCQmJgwcPLnN3d6/h8/nVvXv3rrhw4YKIJMn4I0eOpPn6+kr9/PykAoGg+tKlSyKqDTc3t1r1dVT5wYMHn/j4+MgIgiBv3LghpMovXrwo6t69u0xTEhdtiVbwOXOEENKzrvqc+ahRo7o/ffrU+M0336ykJkZRr6tQKGDAgAGcDz744Nknn3xSkpOTY9S3b1/ezz//nP72229XAwBcvnzZorCw0CgyMrIsNTXVpKysjL5x40bnsWPHls+YMaMUAODp06f00NBQ/r1794Q0Gg2Cg4P5Dx48EDo6Or7wmNeIESN8wsLCyubNm6cxKUtMTIzVt99+63r37t3UrKws40GDBvndv38/5caNG92aKu/WrZuKz+cL1q9fnxMREVEOAHD37l2z27dvd/vss8+Ky8vLaZaWlioajQZxcXHmU6dO7Z6RkZEMUJ+/PT4+PsXV1fWFqWjv379vRqfTyVmzZnlt3bo1Z9CgQQ3JXrQ904/PmSOE0GuGJEkoeFJu3poTttYoLy+nxcfHW/z000+Zv/32m21TdYyMjCAkJKQ6Ly/PGABg69atTpMnTy6mAjkAwIgRI6oiIyPLAAD8/Pzq+vTpI2s8r/qZM2esBw0aVOHs7Kx0dHRUDho0qOL06dPW6nVKS0tpt2/ftgwPD9eakGXSpEkVzs7O8h9++MF+7ty5HsuWLct3dHRUairfu3evfUhISBUVyAEAevXqVfPZZ58VA9RPa0uNt7KyktaSTIohISE1QUFBOp3+FYM5Qgh1QZL0CvPL+5N8JOkVOk+BCgDwyy+/2AwePLg8MDCw1sbGRnnz5k1G4zpSqZS4d+9et7CwsAoAgJSUFPPQ0FCtKUebkpeXZ+zu7l5HLTOZzDrqA4LaeGzfeOONCjs7O1Vz7e3evTtn3bp1zOLiYqO5c+eWaCtPTk42Cw4O1jrmqKgoG29vb8HEiRN99+3bl6m+btiwYb4CgYC3detWhxbubptgMEcIoS7IpbuVbMQs/ycu3a30kgL1xIkTdu+//34pAMDEiRNLjhw50pBMJScnx5TL5fLt7e17MJnMuj59+jQ5hsDAQG737t0FM2bM8NDWV1NXFxqfAZ84ccJu6tSpJS9VbIKXl5f8jTfeqJg1a1ZRS8rVvf322z6+vr6C4cOH+1Bl06dPL8vIyEj+9ddf01atWsWkyv/++2+RUChMuXLlyuP9+/c7Xbx4UW95zzGYI4RQF0QQBLj6WMtactm3tSQSCf2ff/6xmjt3rieTyQz4/vvvXc6dO2erUtWfFHt4eNSKRCKhSCRKun//vsUvv/xiDQDA4/Fk9+7daziDT0hIEK1cuTK/oqLipZSj6tzd3eW5ubkm1HJeXp6Jm5ubXH08CQkJ3SZPnlzedAsvo9FoQKe/3G3jcoFAUPPgwYOGMV+9evXJwYMHM8rKyl6aQXXkyJFVWVlZpgUFBUYA9R8OAACYTKZi9OjRZbdv3+7W0vG1FgZzhBBCrXLkyBHbCRMmFOfn5yfm5eUlSiSSBHd397orV668cObp6ekpX7t2be6WLVtcAQAWLVpUGB0dbX/16tWGoFZdXd1sHBo3blz5jRs3rIqKiuhFRUX0GzduWI0bN64hcEdFRdkNHTq0jMFgNJzCX79+nTF+/Hiv9u7rrFmziuPj4xs+kDQec1JSkin1IebmzZsMuVxOODs7KyoqKmilpaU0AICKigra9evXrQIDA/VylQQAgzlCCKFWOnnypP2ECRNeuNHs3XffLVW/1E6ZNm1amUwmo126dMmCxWIpjhw5kr5ixQp3FovlHxwczD19+rTt/PnzCwEAbty4wXB2dg68cOGC7YIFCzzZbLYAAMDZ2Vn5xRdf5IeGhvJCQ0N5S5YsyadSjgIAxMTE2IWHh79wiT0zM9PU3Ny83Xf/WVhYkGfPnk3bt2+fo7u7e0CPHj2469atc12xYkU+AMDx48dtORyOgMvl8ufNm8c6cuRIOo1Gg9zcXKO+ffty/fz8+CEhIbzhw4eXTZo0qQKg/jt2Z2fnwIcPH3YbP36874ABA3zbO058NA0hhPSsqz6a1pnNnj3b/aOPPirW9H19Z4WPpiGEENKLzpYCtSX27t2ba2iB/NKlSxZhYWFsW1tbRfO1X4QpUBFCCGmFKVBfjXfeeadKLBYL27Jtq87MFXW1DVPHFWamv/QaANq0rq1tqHuVfTXW0jHpYrwtHUdL22hLvbZup+/fpS7a0PX7pq1tIN3S9/Ft7r1oQqfr5VlvhCitCuZlEgkUZWVAUVYGnNu24aXXANCmdW1tQ92r7Kuxlo5JF+Nt6Tha2kZb6rV1O33/LnXRhq7fN21tA+mWvo9vc+9FewuGTxObIaQzrboBrkdgAPngUf2VlqKsDHD09H7hNUEQQJJkq9e1tQ315ycbl+uzr8ZaOiZt27R0vNq0pQ1t+9+W8bamfV3/LnVxfHX9vmlrG0i3Wvp/WR/tkyQJZsbGwlqFQtDW9jv6BjjUOWi7AQ7vZkcIIT0z9LvZUeeAd7MjhBBqs9TUVBMzM7OQxnezR0VF2RAEEfrgwQMzqkypVMKHH37o4evrK+BwOHx/f3+eSCQyAajPIjZixIiGrxx++ukn24kTJ3oBABw9etSGw+HwuVwu39/fn3f58uWGCWhiYmKsvLy8/Fkslv+KFStctI01Pj7ezMvLy7+qqqrhEsngwYPZS5cudWmqfN++fbZUHwEBATxvb28Bl8vljx49uvvjx49NAADmz5/vRo2tf//+vpmZmcaNjwuXy+WHh4ezqLY//fRTpouLSyCDwQhWH9+aNWucXF1dA6ZPn84CHcJgjhBCqFnUFK3qZb/++qtdSEhIlfpkMQcOHLCTSCTGIpEoWSwWC8+ePZtmb2/fMMFLYmIiIz4+3gwaCQsLq3g+Bazw4MGDmXPmzPEEAFAoFLBgwQLWhQsXxGKxOPnUqVN29+7de2l7Ss+ePWtGjRpVumLFClcAgCNHjtjI5XJi06ZNkqbKP/7449K7d++aLVq0iHX48OGMjIyMZJFIJAwPDy9OS0szAQBYvXq1RCwWC0UikXDkyJHlVBvqx0UkEgmPHTuWTZWPGzeuLC4uLqXx+FavXl24fPny/JYd9ZbDYI4QQl0QSZKQL0555SlQCwoKjJ2dneXU/OY+Pj5y9bzjc+fOfbp27VrXxu1pSiX6559/dvP09Kzl8/l1ZmZm5IQJE0piYmJstI1t06ZNBefOnbO7deuW+apVq5h79uzJ1la+fv1614ULFxaEhITUUG1ERESUjxw5sgoAQD0TW3V1dYvSnA4bNqza09NT3mxFHcFgjhBCXVDBY5F57H82+hQ8Fr3SFKiRkZElv//+uw2Xy+XPmjXL/e+//36h/+nTp5ckJSUxkpKSTBu32VQq0ZycHBMmk9mQ/tTd3b0uLy/PpPG26iwtLVUbNmzIGT58OHf8+PGlAQEBtdrKxWKxWe/evbWmOaUum8fExNhv2bKl4cw6NzfXhMfj8Xv16uV36dIlvWVFaw4Gc4QQ6oJcfbmysAXLnrj6cl9pClQfHx95Wlpa0tq1a3NpNBqMGjXK7+zZs5bUdkZGRvDZZ59J1q5d+9J3302lEm3qygJBEM1ebggPDy+3tLRULFq0qLAl5RSJRELncrl8Ly8v/1WrVjlT5bt27cqTSCQJkyZNKt6yZYsTAACLxZJnZGQkpKSkCLdv357z4Ycfdi8pKemQuIrBHCGEuiCCIMCNw+uQFKjm5ubk5MmTK/bu3Zs7f/78gtOnT9uob//JJ5+UxMXFWWZlZTV5hq2eSpTFYr1wJp6bm/tC+lNtaDQaUJfutZVzOJyaO3fuMAAAXFxclCKRSDh9+vSiqqqql3Kkzpgxo+T8+fO21H66uLgoAQAGDhwoZbFYtUlJSRq/z9cnDOYIIYRaRVsK1Js3bzKou72VSiUkJiaae3p61qlvb2pqSn7yySdP9+7d60SVaUol+uabb1ZnZmaaiUQik5qaGuL06dN2EydOLAMA2LBhg+OGDRsc27s/K1askGzbts31/v37DYFYKpU2xMfExMSGrwROnjxp4+PjIwMAyM/PN1Io6qdRFwqFJpmZmaZ+fn617R1PW+Dc7AghhFrl5MmT9kuWLClQL6NSoI4fP75s9uzZnnV1dTQAgB49elQvW7bspUva8+fPf7Z9+/aGG+GOHz9uGx0dbW9kZESamZmpqFSiNBoNtm3blv3OO+9wlEolhIeHP+vZs2cNAIBIJDLv379/VXv3p3fv3rLNmzfnTJ8+3bu6uppma2urZDKZtevXr88HAFi8eLF7enq6GUEQpLu7e93BgwezAACuXLlisW7dOiadTifpdDr53XffZVGpWefMmeP+22+/2dXU1NCcnZ0DIyIinm3fvl3nd7FTcNIYhBDSM0OfNKazpkAdMmQI++LFi0/MzMwMKqnBzp077ePj47tFRUVlN1/7f3Q2aQwpVzUkhajLr3rpNQC8tIxQZ6Pt/dve9rS139J6umijq/ZlCONt6v3BMDYz6EQrnTUF6vXr19MMLZCvWbPGafv27a5WVlbK5mu3XOuyphXLQF5QDfKCaig+mvLSawB4aRmhzkbb+7e97Wlrv6X1dNFGV+3LEMbb1PuDZeNm0IlWqBSojSeNQa23evXqwszMzKTvv/8+T5fttuoye2hQCBn/8B4A1L9BjV27vfCaSh6hvoxQZ6P+HgWAdr9fG7/nNbXfuC9t42hvG121L0MYb1OJVixMGcLqOhkmWkHtgolWEEKoAxn6d+aoc8BEKwghhFAXhsEcIYQQMnAYzBFCCGmlnupTpVJBaGio34kTJ6yo9QcOHLAdMGCAb1PlAwcO9AUAyMnJMQoLC/N2d3cPEAgEvB49enCjoqJsAACuX7/OoNKI+vn58alyAIDevXv7eXl5+VPr8/LyjAAALl68aMHn83lGRkahP/30U0OiF22+++47ew6Hw+dwOHxfX1/B0aNHbQAAJk6c6MVkMgOoPoKDg7kA9Y+Q2draBnG5XL63t7dgzZo1DZPcbN682ZFKixoaGupHZXK7deuWeY8ePbhsNlvA4XD4+/fvbxjb2LFjva2trXu0dLytgZPGIIQQapZ6CtQ9e/ZkTZkyxWfMmDFChUJBfPPNN8wLFy48lkqlRFPlKpUKwsLC2OHh4cWxsbEZAABisdjk5MmTNgD1aUsTExOFxsbGkJWVZRwcHMx///33y4yNjQEAICoqKn3QoEEvJELp3r173U8//ZS5ceNGZ2iBJ0+eGG/bts314cOHKfb29sry8nJaQUFBQwxct25d7owZM0obbxcWFlYaFRWVLZFI6Dwezz8iIqKUzWbLZ86cWbxkyZIiAIBffvnF+vPPP/f4v//7v8cWFhaqI0eOZAQEBNRmZmYa9+rVizd+/PgKBwcH5blz5zKo/O26hsEcIYS6IJIkoS670tyEZanz+dl79epVM3z48PKVK1e6VFdX0ydPnlwsEAhqAQCaKj979qylsbExSQU/AAAOh1P35ZdfFgLUZzOjymUyGdGS8fr5+dUBQJNzrzeloKDAuFu3bipra2slQH3KVWtr67rmtqO4uLgoWSxWbU5OjjGbzZarp0WtqqqiU2MODAxsmM7Vy8tLbmdnpygoKDBycHDQ6XPljWEwRwihLqguu9K85FiKj10474mpp5XOM6dt3rw5PzAwkG9iYqJ69OhRirbyxMRE88DAQK0pRv/4449uH3/8sVd+fr7Jnj17MqizcgCAmTNnetFoNAgLCyvdtGlTQUsDuLq+fftKHRwc5B4eHgH9+/evnDBhQml4eHg5tf6rr75y37RpkysAAIfDkZ07dy5DffvHjx+b1NbW0vr06dNwLL/99lvH3bt3O8vlctrVq1dTG/d5/fp1hlwuJ/h8vt7na8fvzBFCqAsyYVnK7MJ5T0xYlnpJgWplZaUaN25cyeTJk4vNzc3J5srVRUZGsvz8/Pj+/v48qmzo0KHVaWlpyTdv3kzZsmWLq1QqJQAAoqOj08VisfD27duiW7duWezevdu+LeM1MjKCv/766/GxY8ee+Pr61ixbtsxj4cKFbtT6devW5YpEIqFIJBKqB/LY2FhbNpst4PF4AZ988slTBoPRsE/Lly8vysnJSfr6669zV69e7areX1ZWlvGMGTO679+/P5NOfyn5ms5hMEcIoS6IIAgw9bTSSwpUSktTjAYEBMgSEhIY1PKRI0ey//zzT3FpaelLV4dDQkJqGAyGMj4+3hwAwNvbWw4AYGtrq5oyZUrJnTt3urVnvEOGDJF+++23kqNHj6afP3/eprltwsLCStPS0pIvXbokWrVqlUd2dvZLY541a1bJ1atXG9oqKSmhjRw5kr1q1aq8YcOGvZLpUDGYI4QQ0quwsLDK2tpaYtOmTQ3pSquqqhrij0gkMpHL61OUi8Vik4yMDDNfX986uVwO1E1qtbW1xIULF6z9/f2bvdLg7e390mx7mZmZxjdv3mz4QBEfH89gMpkt/s78rbfeqp4wYULxpk2bnAFeTIsaHR1t7enpWQsAUFNTQ4wePZo9derU4o8++uilG+r0xeCCOUmSUFBQ8FJCg8bl6suaXmurp4s2XmVfr+N4dfG+aWs9hFDL0Wg0iI2NffJ///d/lkwmMyAgIIA3bdo0r6+//joXAODatWsWPB5PwOVy+ePGjfPZtm1btqurq0Imk9HeeustXw6HwxcIBHxXV1f5woULiwAAbty4wXB2dg68cOGC7YIFCzzZbLYAAKCgoMCIJMmXLkfU1dURixcvdvf29hZwuVx+TEyM7ffff59Drf/qq6/cqUfTuFwuv6am5qU2Vq9eLYmOjnYoLS2lbd++3YnNZgu4XC5/x44dzj///HMGAMChQ4ds7969a3Hs2DEHqq1bt27pPdGOwU3nWlBQACdOnIDJkyeDq6urxnL1ZQBo8rW2erpo41X29TqOV/3339b3TVvrIdQahj6da2dNgdqU48ePWz958sT0q6++eimHemcwceJErzFjxpQ39Rhcc7rU3OwkSYJEIgEXF5cXEho0LldfBoAmX2urp4s2XmVfr+N4W/NdoKb3TVvrIdQahh7M09LSjAcMGMCzsbFRYOa0ths7dqz3vXv3LLZu3Zr9/vvvlze/xYu6VDBHCCFDY+jBHHUOmGgFIYQQ6sIwmCOEEEIGDoM5QgghZOAwmCOEEEIGDoM5QgghrdRToFJldDo9lMvl8n19fQUjR47sXllZSWtcPnToUPazZ88a5jKNj48369u3L8fLy8vf09PT/4svvnBVqerzlTx48MCsR48eXBMTk5BVq1a9kAktJibGysvLy5/FYvmvWLHChSqfPXu2u7e3t4DD4fDffvttH/W+GuvsqVuTk5NNuVwun8FgBLfiV9MAgzlCCKFmqadABQAwNTVViUQi4ePHj5ONjY3Jbdu2OTYut7GxUWzZssURAKCqqooYP348e8mSJZLMzMykpKQkYVxcnAU1K5yTk5Nix44d2bNnz36q3q9CoYAFCxawLly4IBaLxcmnTp2yo3KHjxgxokIsFieLxWIhm82uWblypQtoQKPRYM+ePVnLli3zkEqlREVFBe2bb75h7t27N7up8j179mRTqVsHDhxYlZubm5icnJxy4sSJ9JycHBOA/6VuFYlEwitXrjz+/PPPPamZ7ADqU7dS870zmUwFwP9St4aFhRWrj08gENS257E/zJqGEEJdEEmSkJuba+7u7q7X+dkBAAYMGFCVkJDw0ixnffv2rabK9+/fb9+zZ8+qCRMmVADUpz398ccfs4cNG+a3fPnyIiaTqWAymYqzZ8/aqLfx559/dvP09Kzl8/l1AAATJkwoiYmJsQkNDZVQbQEA9OvXrzomJsYWtDDE1K0thcEcIYS6oNzcXPOTJ0/6vPfee088PDz0kjkNAEAul8Ply5ethg8fXqFerlAo4Pr165b/+te/ngEAJCcnm4WEhLyQBlUgENRKpVJaSUkJTT0/uLqcnBwT9TnU3d3d6+Li4iwa1/v5558dJk2aVNLceA0tdWtLYTBHCKEuyN3dXfbee+89cXd310sgr62tpVHfoffp06dy/vz5z9TL8/LyTPz9/aXjxo2rAAAgSVLjmWtzszI2Uf+FwqVLl7rQ6XRyzpw5zQZzKkWrhYWFsqnUrY3L1UVGRrLu3LljYWxsTCYlJaUA/C916/37980++OAD70mTJpUzGAwyOjo63dvbW15aWkobM2aMz+7du+3nzZtX3FS7uqDX78xJkoTKSuELCTOoZU2vddVXW7bTdRst3f/O0oa2/e/I8epTa8bbljbbemwQai+CIMDDw0Nvl9ip78ZFIpHw8OHDOWZmZqR6eWZmZmJdXR2xceNGJwAAgUAgu3fvHkO9DaFQaMJgMFS2trZNnpUDALBYrLq8vDwTajk3N9fEzc2t4YvpXbt22V++fNnm9OnTGS098zW01K0toddgXlWVAgmJ/4aqqpSXljW91lVfbdlO1220dP87Sxva9r8jx6tPrRlvW9ps67FByNDZ29srd+7cmf3DDz8419bWEh9//HHx3bt3Lc+cOWMJUH9D3Ny5c1mffvqpRFs7b775ZnVmZqaZSCQyqampIU6fPm03ceLEMoD6u9y/++47lwsXLqSpf3+dkZFh3K9fP0579+FVp25tF+qsoCU/oaGhZGuoVCqyoiKZVKlULy1ret1WbW1DF+No6X61Zd2rbEPb/nfkePWpNeNtS5ttPTaoawGAeLIVf2sb/zx8+DCTJMn4jvoRiUQJbDZbpl5mbm6ubKpu4/IhQ4aUff/99+kkScbHxcUl9+rVq9LT07PGw8OjZuHChXlKpTKeJMn4rKysh05OTnXdunVTWlhYKJycnOqKi4vvkyQZ/+uvvz729PSscXd3r1myZEku1baHh0eNs7NznZ+fn9TPz0/6/vvvF5IkGX/jxg1h//79yzXtz4IFC/JXrlyZ05LyzMzMR6NHjy5xc3Or9ff3r+7du3fFvn37npAkGf/999+n+/j4yPz8/KQ8Hq86KioqjSTJ+PLy8vt8Pr/a19dX6uPjI/vwww+fyuXyeJIk4//880+hk5NTnZmZmdLa2lrh4+PTouNKkmT88/dBk+8RTLSCEEJ6ZuiJVgwpBSoAwIYNGxw9PT3rIiIiWp2ZrKMxGIxgqVT6oKl12hKt4A1wCCGEtKLT6WRlZSWdy+XyDSEF6ooVK4qar9W5JCcnm06cONHH3t5e3nztl2EwRwghpBWbzZZLJJKEjh5HV9beSWNwBjiEEELIwGEwRwghhAwcBnOEEELIwGEwRwghhAwcBnOEEEJaNZUCVVtq0PPnz1sOGTKE3bid3r17+/n7+/Oo5b/++ovRu3dvPwAAiURC79OnD4fBYARPnz6d1Xi71qQT1SQsLMxbfQKYP/74oxuHw+FrKpfL5VBeXk6LiIhgeXh4+PN4PL5AIOBt27bNAaB+ohiBQMDjcrl8Npst2Lx5c0MbEydO9GIymQHUmG/dumUOoDnVa1VVFcHlcvnGxsYh1IQzrYF3syOEEGqWegpUKjVoeHh4cWxsbAZAfWA7efKkTXPtFBcXG504ccJq8uTJLyRmYTAY5Nq1a/MfPXpknpSU9FIGtqioqPRBgwa9kPSESie6ceNG58b1m7J79+6cfv368SIjI0udnZ0Vn332GWvXrl3ZfD6/pqlyY2NjiIiI8PLy8qrNzMxMotPpkJ+fb/TDDz84AACwWCx5fHy8yNzcnCwvL6fx+XzB5MmTy7y8vOQAAOvWrcudMWNGqfoYqFSvjTO8WVhYkM9TpQa0ZF8aw2COEEJdEEmqoLziobm1VQ8ZQej2ImxsbKzW1KDazJs37+nGjRvdGgdzKysr1YgRI6pSU1NNWzqO1qYT9fDwUMybN0/y2Wefuffs2bOaz+dLR4wYUfV8XC+VJycnmz58+LDb2bNn0+l0OgAAuLm5KdavXy8BAKDmoweoT4GqUmmcYr6BplSv7dVhl9lJkoSkSqnGpB7Uusb1NK3TRb2WjrGtfbWljdYcQ30em7aOVxfjaMnvRN/HECFDU17x0Dwp6VOf8oqHL53ltldLUoNqMnDgwCoTExNVbGysZWu2mzlzpheXy+V/8cUXri0Jmpp88cUXRWKx2GzXrl0uO3fuzNVW/vDhQzMejyelAnlT0tLSjDkcDt/b2zvws88+k1Bn5QAAa9asYXI4HP6//vUvD5lMptek8h0WzJOrZPCvpExIrnp57nn1dY3raVqni3otHWNb+2pLG605hvo8Nm0dry7G0ZLfib6PIUKGxtqqh8zff9cTa6seen8DR0ZGsvz8/Pjq34drs2LFioINGza4trT96OjodLFYLLx9+7bo1q1bFrt377Zv61jpdDp89NFHRUOGDCl3cXFRNleubunSpS5cLpfv5OQUSJWx2Wy5WCwWpqSkJB07dswhJyfHCABg+/bteenp6UmPHj1KKS0tpa9cudKlrWNuiQ4L5gILczjo7wUCi5c/NKqva1xP0zpd1GvpGNvaV1vaaM0x1Oexaet4dTGOlvxO9H0METI0BEEDG+sQnV9iB2hdatCmjB07trK2tpZ28+bNFqUF1XU60ZamQA0KCqpJSUlhKJX1sX3Tpk0SkUgkrKqqeulU3cvLS+7n5yf7/fffLQEAPD095TQaDczNzcmPPvqo+N69e4abAlUbgiDA35LRZFJ69XWN62lap4t6LR1jW/tqSxutOYb6PDZtHa8uxtGS34m+jyFC6H+aSw3aEkuXLi3YtWtXs2erbU0n6u3tLWjNeJri7+9fGxgYWD1//nymQqEAAACpVNqQoOzJkyfGVVVVBABAUVERPT4+3kIgENQAAGRlZRkD1N8sePr0aRsej6fXKyR4AxxCCKFWodFoEBsb+2Tu3LkeO3fudLGzs1MwGAzl119/3fAd9O3bt62cnZ0bLkf/8ssvT9TbmDJlSvnatWsV6mVMJjOgqqqKLpfLicuXL9tcuHBB7OvrW/fWW2/5yuVyQqVSEQMHDqxYuHBhEQDAjRs3GJMnT2ZXVFTQr127ZrN+/Xq3tLS05IKCAiOSJHXyKfzo0aOZ8+bN8/D09AywsbFRmJmZqVauXJkLAJCQkGC+dOlSd4IggCRJmDdvnqR3796y5/vnXVJSYkSSJMHn86VRUVFZAADZ2dlGvXr14ldXV9MJgiD37t3rnJKSkmRnZ9f2GwEAMAUqQgjpG6ZAfbWOHz9u/eTJE9Ovvvqq2bvrOxsmkxkQHx+f4urqqmi8DlOgIoQQajNDS4H6/vvvG1we86qqKqJnz548uVxO0Gi0Vj9Cg8EcIYSQVpgCVf+oSWPauj1O54oQQggZOAzmCCGEkIHDYI4QQggZOAzmCCGEkIHDYI4QQkirplKgMhiMYOq1thSi3333nT2Hw+FzOBy+r6+v4OjRozYAALNnz3b39vYWcDgc/ttvv+3z7NmzhlnV4uLizHv06MFls9kCDofDl0qlGp8Z37p1q8Po0aO7U8slJSU0Dw8P/y1btjRZLhKJTAAAvv76a2eqfz8/P/7MmTPda2trCQCAgQMH+vr5+fHZbLYgPDycRU0Ys3PnTntbW9sgKq3p9u3bHaj2Bw4c6GtpadmjcerXsWPHeltbW/doSYrW9jC4YE6SJCTnl7cogQZVT9s22uppWqfr9nTRV2vaaG+9jmwDIdQx1FOgNrZ79+6cXbt2ueTn5xsplUqgUohmZ2cbb9u2zfX27dupYrFYGB8fn9KzZ08pAMCIESMqxGJxslgsFrLZ7Bpq7nK5XA6RkZHeP/74Y1ZaWlryX3/9lWpiYqLxP//ChQufFRQUmJw5c8YSAGDx4sXM8PDwZ4sWLWqynMvl1m3evNnx2rVrVnfv3hWJxWLho0ePUpycnBTV1dUEAMDZs2efpKamCsVicXJxcbHxoUOHGgJxWFhYqUgkEopEIuHChQsbnv1fvHixZO/evRmNx3fu3LmMt956q6xNB70VDC6YCwsq4JOj90FYUNHietq20VZP0zpdt6eLvlrTRnvrdWQbCKGWUZEk3C2vNle9gg/B6qlFt2zZ4kilEC0oKDDu1q2bytraWgkAYG1treJyuXUAABMmTKgwNjYGAIB+/fpV5+XlmQAAnD592prH48n69esnAwBwcXFRGhlpfoqaRqPBjz/+mPXFF1+w/vrrL8bNmzct16xZ81RTOQDA9u3bXffv35/l4OCgBKhPZbphwwYJNQsb9a9cLifkcjnRkimd33333UorK6t2zeLWLtSZUEt+QkNDyY6mUqnIpLwyUqVStbietm201dO0Ttft6aKv1rTR3nod2QZChggA4slW/K1t/PPw4cNMkiTjW/Nzp6wqucffSTV3yqqSW7tt4x+RSJTAZrNl6mXm5uZK9WWFQhEfEBBQ5ebmVltQUPCAJMl4uVwe379//3IXF5faiRMnPvvll18eN9X+kCFDyn744Yd0kiTj16xZk/3uu+8W9+/fv5zH41V/+eWXOS0Z48yZMyUWFhaKc+fOpWorLykpuW9paalorr3+/fuXW1paKsaMGVMsl8vjSZKM37FjR4aDg0Odr6+vdMSIESWPHz9+pL5NbGxs6uDBg8satzVhwoRnhw4detLe38Pz90GT7xGDOzMnCAIEbtYtSqBB1dO2jbZ6mtbpuj1d9NWaNtpbryPbQAi1TKgVQ7ZP4PUk1IrxSnL4NpVC1MjICP7666/Hx44de+Lr61uzbNkyj4ULF7qpb7d06VIXOp1OzpkzpwQAQKFQEHfv3rU4efJkRlxcXOr58+dtz54922zu8wULFhQ6OTnJw8LCKrWVkyT5wt+SU6dOWXG5XD6TyQy4evVqQ2azmzdvPpZIJI/q6uposbGxVgAAkydPLsvOzk4Ui8XCoUOHVk6bNs27HYdMpwwumCOEEGoejSCgl3U3Ge0VfghuKrUojUaDIUOGSL/99lvJ0aNH08+fP29Drdu1a5f95cuXbU6fPp1Bbefu7l7Xt2/fSldXV4WlpaXq7bffLo+Pj2dAM+h0epNpTRuX29nZqczNzVXUjXATJ06sEIlEQg6HI6utrX2hAQaDQY4ZM6bst99+swGov+Rvbm5OAgAsXLiwKDk5udlxvSoYzBFCCOlFZmam8c2bNxsCXnx8PIPJZNYBAMTExFh99913LhcuXEiztLRs+K55/PjxFSkpKeaVlZU0uVwOf//9tyWVVnT8+PFe169fb3cA/fzzzwtmzZrlSd1Br1KpgArk5eXlNCp9qVwuh0uXLllzuVwZwP/SmgIAHDt2zKZ79+417R2LruDc7AghhFqtpqaGpp7i9JNPPnlqZ2f3Qqavuro6YvHixe5Pnz41NjU1Je3s7OT79+/PBgBYuHAhq66ujjZ06FAOAEBISEjVsWPHsh0dHZXz5s17GhwczCMIAoYNG1Y+derUcgCAlJQUhoeHh7y9Y1+yZEmRVCql9ezZk2diYqLq1q2bqnfv3lX9+vWTVlRU0EaPHs2uq6sjVCoV0b9//4ovvviiCABg8+bNTpcvX7ah0+mkjY2N4ueff86k2gwNDfVLT083k8lkdGdn58Ddu3dnTpw48ZXdtYspUBFCSM8wBWr7lZSU0CIiIrwuXryY3lFjaKuJEyd6jRkzpnzGjBml7WlHWwpUvMyOEEJIK/UUqB01Bjs7O5UhBvKxY8d6//PPP5ZmZmZ6fWwNL7MjhBDSClOgtt25c+demkhGH/DMHCGEEDJwGMwRQgghA4fBHCGEEDJwnT+YkyRAQUL9v42XNb3WRb2m+tbXmLTtsy7aaOnxbWtfLdXW46brvlqzHUIIGYDOH8wliQAnIuv/bbys6bUu6jXVt77GpG2fddFGS49vW/tqqbYeN1331ZrtEEJNpkDNyckxCgsL83Z3dw8QCAS8Hj16cKOiomwAAM6fP2/ZOBXozp077adPn85SL+vdu7ffX3/9xQAA+PTTT5kuLi6B6qlVAQBkMhkxevTo7iwWyz8wMJCbmppqAgAgFotNBAIBj8vl8tlstmDz5s2OoMXp06etevTowVWp6m8qVygUwOVy+YsXL3Ztqpya2nX37t12HA6Hz2azBX5+fvwpU6Y0TDYzefJkTz8/Pz6Hw+G/88473cvLy2nU/ltaWvagUqUuXrzYlRrHe++952VnZxfk6+srUB/f7Nmz3R0cHIJWrVrl3OwvpCmtmey/QxKtqFQkmf+o/t/Gy5pe66JeU33ra0za9lkXbbT0+La1r5Zq63HTdV+t2Q4hHYAOSLSiy5/GiVaUSmV8UFBQ1aZNm7KostTU1IR169ZlkxoSjuzYsSMjMjKyUL2sV69elTdu3BCSJBn/+++/p2RmZj5qnMDl22+/zXr//fcLSZKM37t375NRo0aVkCQZL5PJ7kml0nskScaXlZXdd3Nzq83IyHikbT9Gjx5dsm3btkzyeUKXKVOmFGkrP3nypJjP51enp6c/Ip8njvnPf/6T8fDhw0SSJOOLi4vvU23/61//kixfvjxX0/5TPxcuXBD93//9n7Bx4hqSJOMXLFiQv3LlSo2JZbQlWun8j6YRBIBroOZlTa91Ua/xOn2OSVv7umhDE1301VJtPW766Kul2yFkoFQkCQ+yy8yDWTY6n589NjbW0tjYmFyyZEkRVcbhcOq+/PLLwra2OWzYsOqmys+fP2/z9ddf5wMAzJgxo3Tp0qUslUoFZmZmDd+FyWQygjqz1uaHH37IGTBggN+gQYOqDhw44HT37t0UbeXffvut68aNG3O9vb3lAPWJYz7//PNiqj0qVapKpQKZTEZrSTKokSNHVlFXF3Sp819mRwgh1GoPssvM5x677/Mgu8xc120nJiaaBwYGSnXdblOePn1q4u3tXQcAYGxsDBYWFsqnT58aAQCkpaUZczgcvre3d+Bnn30m8fLy0jrVq6enp3zOnDmFgwcP5i1evLjA2dlZqa08LS3N/I033tC6n5MmTfJydHQMSktLM1u2bFnDh5kHDx5Y+Pn58QcNGuQbHx9v1t7j0BwM5ggh1AUFs2xkP4SHPAlm2eg9BWpkZCTLz8+P7+/vz9NUR9NZa3Nns2QTN6MSBEEC1E9mIxaLhSkpKUnHjh1zyMnJafZq87JlywqVSiV89tlnxS0pp9y5c8ecy+XyPTw8/Pfv329LlcfExGQ+ffr0ka+vb82hQ4dsAQDeeOON6qysrITU1FTh3LlzCydOnMhuqk1dwmCOEEJdEI0gINTTVi8pUAMCAmQJCQkN2cuOHDmS/eeff4pLS0s1BlMHBwdFWVkZXb2srKyM7uzsrNC0DQCAi4tLXUZGhglAfRazqqoqupOTk1K9jpeXl9zPz0/2+++/N5v3nE6nN/kBoqlyNpstu3XrFgMAoHfv3jKRSCQcMmRIhUwmeyF2GhkZwfvvv19y5swZW4D6y+/W1tYqAIApU6aUKxQKoqCgQK9fa2MwRwgh1CphYWGVtbW1xKZNmxruIK+qqtIaTwYMGFB97949i+zsbCMAgL/++otRV1dH8/HxqdO23ejRo8sOHTpkDwDw008/2fbr16+SRqPBkydPjKuqqggAgKKiInp8fLwFlSp17ty5TOrO+vZYsmSJZNmyZe5PnjxpSH1aU1NDANR/T56UlGRKvT579qyNr69vDQBAdna2EfUd/vXr1xkqlQqa+9DSXp3/BjiEEEKdCo1Gg9jY2Cdz58712Llzp4udnZ2CwWAov/7661yqzu3bt63UU6T+8ssvTzZt2pTzzjvv+KpUKqJbt27Ko0ePptPp9Sfrc+bMcf/tt9/sqNSqERERz7Zv354/f/78ZxMnTvRmsVj+1tbWyujo6CcAAAkJCeZLly51JwgCSJKEefPmSXr37i0DABAKhebjx48va+9+TpkypbywsNBo5MiRvkqlkrCyslJyuVzZu+++W0GSJEyfPt27qqqKRpIkwePxpD///HMWAMDRo0dtDx065ESn00kzMzNVVFRUOo1W/1knLCzM+59//rEsLS01cnZ2Dly2bFn+ggUL2p0RD1OgIoSQnmEK1FdrwIABvjdv3nzc0eNorYULF7pZWFgo165d+7Sp9ZgCFSGEUJt1hhSorWGIgXz27NnuMTExdt26dWtTqlS8zI4QQkgrTIGqf3v37s0FgNxmK2qAZ+YIIYSQgesUwZwkSRCViJp8nrAjqY9L02td1Gtpv21dp62epn713Zeu2jAkXWU/EEKdT6cI5qmlqbDg+gJILU3t6KG8QH1cml7rol5L+23rOm31NPWr77501YYh6Sr7gRDqfDrF3ewkSUJqaSr42fo1OxvQq6Q+LgBo8jX1WER76jXeZ2312rJOWz31vhuX67Ovth7fzvT+aK2ush+o9Qz9bnbUOXT6u9kJggCuHbfT/YFTH5em17qo19J+27pOWz1N/eq7L121YUi6yn6g1496ClSJREKnUns6ODgEOTk5BXK5XL6fnx8/ODiYe+LECStquwMHDtgOHDjQFwCATqeHcrlcvq+vr2Do0KFsKo0oAMDAgQN9LS0tezROmyoSiUwCAwO5np6e/qNHj+6uPmHLhx9+6MFisfw5HA7/5s2bDNBi3rx5zE8++YRJLYvFYhN3d/eAadOmsZoqf/bsGV0ul8O8efOYnp6e/tT+Ll261AUAQCqVEgEBATw/Pz8+m80WLFiwwI1qY+HChW7UMeFyufzo6GhrAACJRELv06cPh8FgBDdOBUuVU+lg26JTBHOEEEKdm4eHR61IJBK6uLgoRSKRUCQSCadPn140Z86cpyKRSJiamirct29f5rJlyzykUilRUVFB++abb5h79uzJBgAwNTVViUQi4ePHj5NtbGwUW7ZsaZg9bvHixZK9e/dmNO5z4cKF7vPmzXualZWVZG1trdixY4cDAMDJkyet09PTzTIzM5N+/PHHrH//+9+sxtuq27hxY/6lS5ds7t+/bwYA8O9//9vjyy+/zNuzZ09OU+UODg7K+fPnMwsKCoxTUlKSRSKR8Pbt2yK5XE4DADAzMyNv3ryZmpqaKkxOThZeu3bN6tq1a92o/qhjIhKJhFOmTCkHAGAwGOTatWvz1SfWocTFxYn9/f3blbgGgzlCCHVFpAogJ84cyDY9ttwmvXr1qhk+fHj5ypUrXZYsWeI2efLkYoFAUNu4Xt++favz8vIa0oC+++67lVZWVi8MVKVSwe3bty1nzJhRCgDw0UcfFcfGxtoAAJw9e9YmIiKimEajwbBhw6orKiqMsrKyjEEDCwsLcuPGjbmffPIJ68SJE1bV1dX0Tz75pERTeWVlJe3YsWOOBw4cyGYwGCQAgK2trWr79u35APUz4FFzr9fV1REKhYJo7oqblZWVasSIEVVmZmZ6+YVgMEcIoa4o9645nPzQB3Lv6jwFqjabN2/OP3XqlP0ff/xhtXbtWknj9QqFAq5fv245bty4Mm3tPH361MjS0lJpbFwfo728vOqePn1qAgBQUFBg7OXl1TCnu6ura522YA5QPzWrjY2Ncvbs2d579uzJ0lYuFApNXV1d62xtbTUGXoVCAVwul+/s7Bz05ptvVgwdOrQhH/vBgwedOBwO/7333vMqKiqia2pDlzCYI4RQV+TeSwbv/fwE3HvpPQWqOisrK9W4ceNKJk+eXGxubt5wh3VtbS2Ny+XybW1te5SVlRmNGzeuQls72lKfaljX7Njmzp1bGBgYWB0UFFTbknLKjh077LlcLt/FxSUwLS3NGKA+U5pIJBJmZ2cn3L9/v9vdu3fNAAAWLFhQmJWVlZiSkiJ0cXGR//vf//ZodmA6gMEcIYS6IoIG4NFHBsSr/zNPo9GASixCob4zz8zMTKyrqyM2btzopK0NFxcXRWVlJV0ulwMAQGZmpomTk5McAMDNzU2emZnZcJm+oKDAhMViydsyrqbK+Xx+bUFBgUlpaSkNAGD+/PnFIpFIaGlpqVQqlS98anBwcFAOGDCgMjY21hoAwMPDQ2FkZAR0Oh3mzZtX9PDhw27wCmAwRwgh9MrY29srd+7cmf3DDz8419bWajydptFo0Ldv38qffvrJFgDg0KFD9mPGjCkDABg7dmzZL7/8Yq9SqeDatWvdLC0tlZ6ennIAgH79+nEyMjK0XnJvjqWlpWrq1KnP/vWvf7GkUikBUH9ZXS6XEwAA+fn5RtTd+FVVVcSff/5pxePxagAA1C/3//rrrzZ+fn6v5MoIzs2OEELolerfv7+Mx+PJDhw4YDt37tyS0NBQv/T0dDOZTEZ3dnYO3L17d+bEiRMrtm3bljtlyhSfdevWMQUCgXT+/PnPAAAmT55c/t///tfa09PT39zcXHXgwIFMAAClUglZWVmmjo6O7c4dvmPHjrwFCxa4cblcQbdu3VRmZmaqKVOmPPP09JQ/evTI7MMPP/RWKpVAkiTx7rvvlrz//vvlAADz5893FwqF5gAA7u7udT/99FPD9/NMJjOgqqqKLpfLicuXL9tcuHBBHBoaWtPesQJ0kkljEEKoKzP0SWMMJQXq3bt3zfbu3etw4MCBNics6Si9e/f227p1a86gQYM0PqLW6SeNQQgh1HkZSgrUXr161RhiIO/Tpw8nJyfHxNjYuM2JGwwumJMkCTUpKTpLVqHenra229Jv42209dXScbRlX1raly7Gq+82dKGl7em6X4QMFZUCVSQSCTt6LF1RXFycuKCgILFfv35t/n7d4IJ5rUgEuZ9+BrUikc7b09Z2W/ptvI22vlo6jrbsS0v70sV49d2GLrS0PV33ixBC+mJw35mTJAm1IhGYcnUzx7V6ewCgse229Nt4G219tXQcbdmXlvbVmjF1VBu6/p1ra0/X/aLXl6F/Z446B23fmRtcMEcIIUODwRzpAt4AhxBCCHVhGMwRQghppZ4CFaD+MapTp05ZqddZu3at07Rp01gA9ZOqGBkZhWzZssVBvQ6TyQwoKCh4YX6To0eP2nA4HD6Xy+X7+/vzLl++bEGte/bsGf2dd97p7u3tLejevbvg999/1zibWnx8vJmXl5d/VVVVw3digwcPZi9dutSlqfJ9+/bZAgDExMRYBQQE8Ly9vQVcLpc/evTo7o8fPzYBAJg/f74bNbb+/fv7ZmZmGjc+Hlwulx8eHt6Qte3TTz9luri4BDIYjGD18a1Zs8bJ1dU1oHH6U13BYI4QQqhZVApUAID33nuv+Pjx43bq60+dOmU3bdq0EgCAqKgo26CgoOqTJ0/aN9duWFhYBZUu9ODBg5lz5szxpNZ9/PHHHsOHD6/IyMhIFgqFwh49emicYKVnz541o0aNKl2xYoUrAMCRI0ds5HI5sWnTJklT5R9//HHp3bt3zRYtWsQ6fPhwRkZGRrJIJBKGh4cXp6WlmQAArF69WiIWi4UikUg4cuTIcqoN9eMhEomEx44dy6bKx40bVxYXF5fSeHyrV68uXL58eX5zx6OtMJgjhFAXpCJV8LDwoblKDylQIyMjS69du2Ytk8kIgPoz1cLCQuPhw4dXAQCcPHnSbuvWrTkSicS4ualVra2tVdS86JWVlTTqZtOSkhJaXFyc5eeff/4MoD6HuIODg1JbW5s2bSo4d+6c3a1bt8xXrVrVkEtdU/n69etdFy5cWBASEtLwISEiIqJ85MiRVQAAdnZ2DQevurqa1pIbYYcNG1ZNTS37KmEwRwihLiihKMF88Y3FPglFCTpPgeri4qIMCgqqPnXqlDUAwOHDh+3Gjh1bSqPRIC0tzfjZs2fGQ4YMkY4dO7b08OHDds21FxUVZePt7S2YOHGi7759+zIBAEQikamdnZ3ivffe8+LxePwpU6Z4VlRUaI1ZlpaWqg0bNuQMHz6cO378+NKAgIBabeVisdisd+/eGmdcA/jfZfOYmBj7LVu2NJxZ5+bmmvB4PH6vXr38Ll26ZKGtjVcBgzlCCHVBgY6Bsq1vbn0S6Biol0QfkydPLomOjrYFADh9+rRdZGRkCcD/AjsAQGRkZElMTEyzwXz69OllGRkZyb/++mvaqlWrmAAACoWCSElJYcydO7coJSVFyGAwVCtXrnRprq3w8PByS0tLxaJFiwpbUk6RSCR0LpfL9/Ly8l+1apUzVb5r1648iUSSMGnSpOItW7Y4AQCwWCx5RkZGQkpKinD79u05H374YfeSkpIOjacYzBFCqAuiETTo4dRDRtNTCtSIiIiyv//+2+rmzZuMmpoa2oABA6QA9d+dR0dH2zOZzIAJEyawU1NTzRMTE01b0ubIkSOrsrKyTAsKCoy8vLzqnJ2d64YOHVoNADBlypTSR48eMVrSTktTnXI4nJo7d+4wAOqvNohEIuH06dOLqqqq6I23nTFjRsn58+dtAQDMzc1JFxcXJQDAwIEDpSwWqzYpKcmsJWPTFwzmCCGEWs3a2lrVt2/fypkzZ3pNmDChBADg0aNHplKplF5YWJiQl5eXmJeXlzhv3jxJVFSUxrPzpKQkU5Wq/qvpmzdvMuRyOeHs7KxgsVgKFxeXukePHpkCAFy5csXKz8+vBgBgw4YNjhs2bHBs7z6sWLFCsm3bNtf79+83BGKpVNoQF9U/hJw8edLGx8dHBlB/t75CUZ+YTSgUmmRmZpr6+fnVtnc87YEpUBFCCLXJ1KlTSz744AOf48ePpwMAHD582H7UqFGljeqUhoeHd9+yZUsBAEBQUBCfupEsLCysxNnZWREdHW1vZGREmpmZqY4cOZJOnT3v2rUrOyIiontdXR3BYrFqjx8/ngkAIBKJzPv371/V3vH37t1btnnz5pzp06d7V1dX02xtbZVMJrN2/fr1+QAAixcvdk9PTzcjCIJ0d3evO3jwYBYAwJUrVyzWrVvHpNPpJJ1OJ7/77rssZ2dnJQDAnDlz3H/77Te7mpoamrOzc2BERMSz7du36+0udsprPwMcSZLwLLcKHNzr71+gXjc1naumeprWtbReZ2rjVe5zW/rCaVWRITL0GeA6WwrUIUOGsC9evPjEzMzMoLIg7dy50z4+Pr5bVFRUdvO1X4YzwGnxLLcKLu1NhGe5VS+8bk09TetaWq8ztfEq97ktfSGEXr3OlgL1+vXraYYWyNesWeO0fft2VysrK62P17UVnpl3krPUztLGq9zntvSFZ+bIEBn6mTnqHDDRCkIIdSAM5kgX8DI7Qggh1IVhMEcIIYQMHAZzhBBCyMBhMEcIIaSVespPlUoFoaGhfidOnGhIgXrgwAHbAQMG+DZVPnDgQF8AgJycHKOwsDBvd3f3AIFAwOvRowc3KirKBgDg+vXrDCqdqJ+fH58qB6hPt+rl5eVPrc/LyzMCALh48aIFn8/nGRkZhf7000+2ze1DYGAgl8vl8l1dXQNsbW2DqPZSU1NNGqdmPX/+vOWQIUPYAPWPk1H1vb29BWvWrHGi6mkbw65du+w9PT39PT09/Xft2vVS9rgPPvjAQz1N6v79+21ZLJY/1W9r4aQxCCGEmqWeAnXPnj1ZU6ZM8RkzZoxQoVAQ33zzDfPChQuPpVIp0VS5SqWCsLAwdnh4eHFsbGwGAIBYLDY5efKkDUB9+tLExEShsbExZGVlGQcHB/Pff//9MmPj+oRrUVFR6YMGDXohIUr37t3rfvrpp8yNGzc6QwskJCSIANr2rHdYWFhpVFRUtkQiofN4PP+IiIhSNpst1zSGp0+f0jdt2uR27949IY1Gg+DgYP7UqVPLHB0dlQAAf/31F6O8vPyF+Dtr1qxSV1dXxbZt21q0P41hMEcIoS6IVKlA9uiRuXlQkIxoYp7y9ujVq1fN8OHDy1euXOlSXV1Nnzx5crFAIKgFAGiq/OzZs5bGxsbkkiVLiqg2OBxO3ZdfflkIUJ/VjCqXyWRESx5B9fPzqwOAJudg1xcXFxcli8WqzcnJMWaz2XJNYzhz5oz1oEGDKqhZ4QYNGlRx+vRp69mzZ5coFAr44osv3E+cOJHB4/FsdDU2DOYIIdQFyR49Ms/7fIEP87v/PGEEB+s8c9rmzZvzAwMD+SYmJqpHjx6laCtPTEw0DwwM1Jpq9I8//uj28ccfe+Xn55vs2bMngzorBwCYOXOmF41Gg7CwsNJNmzYV6COAv/nmmxyqXalUSvPx8alpXOfx48cmtbW1tD59+mg9nnl5ecbu7u511DKTyazLy8szBgD49ttvnUaNGlWm65znGMwRQqgLMg8KkjG/+88T86AgvaRAtbKyUo0bN67EwsJCaW5uTjZXri4yMpJ1584dC2NjYzIpKSkFAGDo0KHVaWlpyffv3zf74IMPvCdNmlTOYDDI6OjodG9vb3lpaSltzJgxPrt377afN29esa7358aNG2JXV1cFQP135uqXu2NjY23ZbLZlZmam2bZt2zIZDIbWCVqamr+FIAjIzMw0PnPmjO0///yTquvx4w1wCCHUBRE0GjCCg3V+iV1dS1ONBgQEyBISEhrSlx45ciT7zz//FJeWlr50QhkSElLDYDCU8fHx5gAA3t7ecgAAW1tb1ZQpU0ru3LnTTR/7ok1YWFhpWlpa8qVLl0SrVq3yyM7O1noi7O7uLs/NzTWhlvPy8kzc3Nzk//zzDyMrK8vMy8srgMlkBtTU1NBYLJa/LsaIwRwhA0KSJBRmpjd88ldf1vRaWz1dtPEq+2prG6hjhYWFVdbW1hKbNm1qSFtaVVXVEH9EIpGJXF5/1VksFptkZGSY+fr61snlcqDuMq+trSUuXLhg7e/v3+yVBm9vb4EedgPeeuut6gkTJhRv2rRJ601q48aNK79x44ZVUVERvaioiH7jxg2rcePGlU+dOrX82bNnj6j0sGZmZqrs7OwkXYwNgzlCBqQoKwPObdsARVkZLy1req2tni7aeJV9tbUN1LFoNBrExsY++b//+z9LJpMZEBAQwJs2bZrX119/nQsAcO3aNQsejyfgcrn8cePG+Wzbti3b1dVVIZPJaG+99ZYvh8PhCwQCvqurq3zhwoVFAAA3btxgODs7B164cMF2wYIFnmw2WwAAUFBQYESSpN6SOKxevVoSHR3tUFpaStM0BmdnZ+UXX3yRHxoaygsNDeUtWbIkn7oZTl9wbnaEDAhJklCUlQGOnt4NSWioZQBo8rW2erpo41X21dY2OjpBj6HPzd7ZUqBqc/z4cesnT56YfvXVV4UdPZbWor6rv379elpT6zHRCkIIdSBDD+ZpaWnGAwYM4NnY2CioZ82Rbu3fv99248aNbgEBAdIzZ840eTlJWzDHu9kRQghpxWaz5RKJJKGjx9GVzZo1q3TWrFmlbd0evzNHCCGEDBwGc4QQQsjAYTBHCCGEDBwGc4QQQsjAYTBHCCGklXoKVID6tKSnTp2yUq+zdu1ap2nTprEAAPLz842MjIxCtmzZ4qBep3GqUQCABw8emPXo0YNrYmISsmrVqhcmY4mJibHy8vLyZ7FY/itWrHChyg8dOmTLZrMFNBot9K+//mKAFp09ZWtycrIpl8vlq6dDbQsM5gghhJqlngL1vffeKz5+/Lid+vpTp07ZTZs2rQQAICoqyjYoKKj65MmTL+XxbszJyUmxY8eO7NmzZz9VL1coFLBgwQLWhQsXxGKxOPnUqVN29+7dMwMA6NGjh+zUqVNpPXv2rGqufRqNBnv27MlatmyZh1QqJSoqKmjffPMNc+/evdlNle/ZsyebStk6cODAqtzc3MTk5OSUEydOpOfk5JgA/C9lq0gkEl65cuXx559/7knNYPd8/9NFIpFQJBIJmUymAuB/KVvDwsJemFdeIBDU6uJxv1Y9mnbv3r1nBEFktbdTA+YAAB32rGcnhMfjf/BYvAiPx4v8XnWHJEmCJL3C3KW7lUzXk+ZERkaWbtiwgSmTyQhzc3MyNTXVpLCw0Hj48OFVAAAnT56027p1a84HH3zQPSMjw5iaX70pTCZTwWQyFWfPnrVRL//zzz+7eXp61vL5/DoAgAkTJpTExMTYhIaGSkJCQl7KaKbN65CytVXBnCRJx+ZrdV0EQcS3Z+KHrgaPx//gsXgRHo8XEQTxymfbkqRXmF/en+QzYpb/E1cfa51mTnNxcVEGBQVVnzp1ynratGllhw8fths7dmwpjUaDtLQ042fPnhkPGTJEOnbs2NLDhw/bff3110+bb/VFOTk5JkwmsyGNqLu7e11cXJxFW8fc1VK2NoaX2RFCqAty6W4lGzHL/4lLdyu9pECdPHlySXR0tC0AwOnTp+0iIyNLAACowA4AEBkZWRITE2OnrR1NNKQRbXPWHCo16+TJk4ubStnauFxdZGQky8/Pj+/v78+jyqiUrTdv3kzZsmWLq1QqJQAAoqOj08VisfD27duiW7duWezevbvZrxp0AYM5Qgh1QQRBgKuPtc4vsVMiIiLK/v77b6ubN28yampqaAMGDJAC1H93Hh0dbc9kMgMmTJjATk1NNU9MTDRtbfssFqsuLy+vIY1obm6uiZubm8bL9S3RlVO2YjBvnX0dPYBOBo/H/+CxeBEejxd1ueNhbW2t6tu3b+XMmTO9JkyYUAIA8OjRI1OpVEovLCxMoNJ8zps3TxIVFdXqs/M333yzOjMz00wkEpnU1NQQp0+ftps4cWKZtm0yMjKM+/Xrx2njLjV41SlbdQGDeSuQJNnl/kO2Bx6P/8Fj8SI8Hi/qqsdj6tSpJampqeZql9jtR40aVdqoTunp06cbgnlQUBDf2dk50NnZOXDmzJnu2dnZRs7OzoH79u1z/s9//uPq7OwcWFJSQjM2NoZt27Zlv/POOxxfX1/BuHHjSnr27FkDABAVFWXj7Owc+PDhw27jx4/3HTBgAPU4mTGdTm93AvtXmbJVV1qVNQ0hhNCr19FZ0wwlBeqGDRscPT096yIiIso7eiytxWAwgqVS6QNtdTBrGkIIoTaj0+lkZWUlncvl8jtzCtQVK1YUNV+rc0lOTjadOHGij729ffvuB9DVgLoSgiA8CIK4ThBECkEQyQRBzH9ebkcQxFWCIB4//9e2uba6EoIg6ARBPCAI4vzz5df2eBAEYUMQRAxBEKLn75N+r+vxIAhiwfP/J0kEQRwnCMLsdToWBEEcIgiikCCIJLUyjftPEMRygiDSCIJIJQhiRMeMunWoFKidOZAbKmrSmJycnKTma2uGwbxpCgBYRJIkDwD6AsBcgiD4ALAMAK6RJOkLANeeL79O5gNAitry63w8dgDAJZIkuQAQBPXH5bU7HgRBMAHgMwDoSZKkPwDQAWAqvF7H4mcAeKdRWZP7//zvyFQAEDzfZjdBEPRXN1TUVWEwbwJJkgUkSd5//roS6v9QMwHgXQA4/LzaYQAY1yED7AAEQbgDwGgAOKBW/FoeD4IgrABgEAAcBAAgSbKOJMkyeE2PB9R/XWdOEIQRADAAIB9eo2NBkuRfAFDSqFjT/r8LAL+SJFlLkmQGAKQBQO9XMU7UtWEwbwZBEF4AEAwAcQDgTJJkAUB9wAcApw4c2qv2HQAsAQCVWtnrejy6A0ARAPz0/GuHAwRBdIPX8HiQJJkHAFsBIBsACgCgnCTJK/AaHotGNO0/EwBy1OrlPi9DqF0wmGtBEIQFAJwCgM9Jkqzo6PF0FIIgxgBAIUmS9zp6LJ2EEQCEAMCPJEkGA0A1dO3LyBo9/y74XQDwBgA3AOhGEMS0jh1Vp9bUDC74SBFqNwzmGhAEYQz1gfwXkiRPPy9+ShCE6/P1rgBQ2FHje8X6A8BYgiAyAeBXABhKEMRReH2PRy4A5JIkGfd8OQbqg/vreDzeAoAMkiSLSJKUA8BpAHgDXs9joU7T/ucCgIdaPXeo/1qiU1NPgdrZU4pqU15eTouIiGB5eHj483g8vkAg4G3bts2B2kdfX18BgOa0qdS+vPfee152dnZBVH1169evd/Ly8vJns9mCOXPmuAPUTyAzYcIELw6Hw+/evbtg+fLlDelc+/Tpw2EwGMHNpXJtDgbzJjxPgXMQAFJIktyutuocAHzw/PUHAHD2VY+tI5AkuZwkSXeSJL2g/uadP0iSnAav7/GQAEAOQRBUJqxhACCE1/N4ZANAX4IgGM//3wyD+ntMXsdjoU7T/p8DgKkEQZgSBOENAL4AcKcDxtdqVArUzp5SVJuIiAgvW1tbZWZmZlJKSorw6tWrj0tKSl56RFvTPu7ZsycbAOCjjz56du7cuceNt4uNjbX873//a5OSkpKclpaWvHLlSgkAwE8//WRbV1dHE4vFwkePHqVERUU5pqammgAAxMXFif39/bUmdWkJfM68af0BIBIAEgmCePi8bAUAbASAEwRB/Avq/4i91zHD6zRe5+PxKQD8QhCECQCkA8AMqP9w/FodD5Ik4wiCiAGA+1D/FMgDqJ+61AJek2NBEMRxABgMAA4EQeQCwGrQ8H+DJMlkgiBOQP2HPwUAzCVJUqmPcZEkCQWPReauvlydz89uiClFk5OTTR8+fNjt7Nmz6XR6/QMEbm5uivXr10tau48jR46sooKxuh9//NFxyZIlBVTCFuqDB0EQIJVKaXK5HKqrqwljY2PSxsZGp793DOZNIEnyJjT93RZA/ZnHa4skyT8B4M/nr4vhNT0eJEk+BICmUny+dseDJMnVUB/A1NXCa3IsSJJ8X8OqJvefJMn1ALBefyOqV/BYZB77n40+YQuWPXHj8HQ+P7ihpRR9+PChGY/Hk1KBvD37qEl6errZjRs3LFetWsU0NTUlt27dmvPmm29KP/zww9LY2FgbJyenoJqaGto333yT4+zsjMEcIYSQdq6+XFnYgmVPXH25ekn0QaUOtbCwUDaVUrRxubrIyEjWnTt3LIyNjcmkpKQUgP+lFL1//77ZBx984D1p0qRyBoNBRkdHp3t7e8tLS0tpY8aM8dm9e7f9vHnzWnxpXZOlS5e6nD171q6kpMSosLAwoTX7qIlSqSRKS0vpDx8+FN24cYMRHh7uk5OTk3jjxg0GjUYjJRJJwrNnz+j9+/fnjho1qoLP59c112ZL4XfmCCHUBREEAW4cnt5SoAIYVkrRoKCgmpSUFIZSWX9CvGnTJolIJBJWVVVpPVXXtI9NcXFxqZs0aVIZjUaDIUOGSJ8HcKMjR47YjxgxotzU1JRkMpmKXr16Vd26dUunqVExmCOEENKrV51S1Nvb+6W7zP39/WsDAwOr58+fz1QoFAAAIJVKdZpsLCwsrOz333+3BABISEgwlcvlNBcXFwWLxaq7fv26lUqlgoqKCtr9+/e7BQQE1OisY8DL7AghhPSMSik6d+5cj507d7rY2dkpGAyGUj2l6JgxY1yNjIxIGo1GUilFKyoqaG+99ZavXC4nVCoVMXDgwAr1lKKTJ09mV1RU0K9du2azfv16t7S0tOSCggIjkiSbvBxx9OjRzHnz5nl4enoG2NjYKMzMzFQrV67MpdZnZGSYOjs7B1LL3377bU5T7YSFhXn/888/lqWlpUbOzs6By5Yty1+wYMGzzz777NmUKVO8fH19BcbGxqp9+/Zl0Gg0WLJkSeHUqVO9OByOgCRJCA8Pf9anTx+dfv2BKVARQqiTwxSoLXf8+HHrJ0+emH711VcGM7dB7969/bZu3ZozaNAgrTcJYgpU1CkQBGEP9UknAABcAEAJ9dOisgEgiiTJf+uhz88BoIQkyahWbmcCAL8DwFCSJBW6HhdChsRQUqACALz//vsGlcu8T58+nJycHBNjY+N2nVljMEevzPNH2XoAABAE8TUAVJEkuVVf/T1P/PER1M/O1iokSdYRBHENAKYAwC+6HhtCraRSqVQEjUbrkEupVArUjui7q4uLixO3pJ5KpSLgxdwYL8Ab4FCHIwhisFqO9K8JgjhMEMQVgiAyCYKYQBDEZoIgEgmCuPR8ml0gCCKUIIgbBEHcIwjiMjV1ZiNDAeA+dWZNEMSfBEFsIgjiDkEQYoIgBj4vFzwve0gQRAJBEL7Ptz8DABF6PwAINS+pqKjI+vkfdPSaUalURFFRkTUAaMx5jmfmqDPyAYAhAMAHgNsAMJEkySUEQfwGAKMJgvgvAOwCgHdJkiwiCGIK1E/C8VGjdvoDQOPkMEYkSfYmCGIU1E908hYAzAGAHSRJUjO6UY+qJAFALz3sH0KtolAoZkokkgMSicQf8CTsdaQCgCSFQjFTUwUM5qgzukiSpJwgiESoD6yXnpcnAoAXAPgBgD8AXH3+DC0d6tNvNuYK9fOEq6OS5tx73hZA/QeGL5/nbD9NkuRjAACSJJUEQdQRBGH5PK89Qh0iNDS0EADGdvQ4UOeFwRx1RrUAACRJqgiCkJP/e+RCBfXvWQIAkkmS7NdMOzIAMGuqbai/+c7oeT/HCIKIA4DRAHCZIIiZJEn+8byeKQDo9HlQhBDSNbxcgwxRKgA4EgTRD6A+XS1BEC9NEgH1Z+Xs5hojCKI7AKSTJLkT6rNaBT4vtwcAKrUnQgh1WhjMkcEhSbIOACYBwCaCIB4BwEOoz6Hd2EUAGNSCJqcAQNLzDHlcAKAeYxsCABfaO16EENI3nDQGdWnPb5pbQn0P3sptTwPAcpIkU3U/MoQQ0h08M0dd3TKovxGuVZ7f1X4GAzlCyBDgmTlCCCFk4PDMHCGEEDJwGMwRQgghA4fBHCGEEDJwGMwRQgghA4fBHCGEEDJw/w89ctRlmTlFGAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "(fig, ax, lgd) = plot_concurrence(concurrence, x_values=times)\n",
    "plt.xlabel(\"Time (ns)\")\n",
    "plt.savefig(\"concurrences.pdf\", bbox_extra_artists=(lgd,), bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This plot shows when each contact occurred. The x-axis is time. Each dot represents that a specific contact pair is present at that time. The contact pairs are separated along the vertical axis, listed in the order from `contact_list` (here, decreasing order of frequency of the contact). The contacts are also listed in that order in the legend, to the right.\n",
    "\n",
    "This trajectory shows two groups of stable contacts between the protein and the ligand; i.e. there is a change in the stable state. This allows us to visually identify the contacts involved in each state. Both states involve the ligand being in contact with Phe33, but the earlier state includes contacts with Ile28, Gly29, etc., while the later state includes contacts with Ser32 and Gly168.\n",
    "\n",
    "This change occurs around 60ns (which is also frame 60), and is evident when viewing the MD trajectory. If you have [NGLView](https://github.com/arose/nglview/) installed, visualize the trajectory with the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for visualization, we need to clean up the trajectory:\n",
    "traj.topology.create_standard_bonds()  # required for image_molecules\n",
    "traj = traj.image_molecules().superpose(traj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "81c8628c0a5740d58140b9522455fd7f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bd9f9289582044a29392e24f0fee93f3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "NGLWidget(frame=60, max_frame=99)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import nglview as nv\n",
    "view = nv.show_mdtraj(traj)\n",
    "view.remove_cartoon()\n",
    "view.add_cartoon('protein', color='#0000BB', opacity=0.3)\n",
    "view.add_ball_and_stick(\"YYG\")\n",
    "\n",
    "# update to my recommeded camera orientation\n",
    "camera_orientation = [-100,  45, -30, 0,\n",
    "                        50,  90, -45, 0,\n",
    "                         0, -45,-100, 0,\n",
    "                       -50, -45, -45, 1]\n",
    "view._set_camera_orientation(camera_orientation)\n",
    "\n",
    "# start trajectory at frame 60: you can scrub forward or backward\n",
    "view.frame = 60\n",
    "view"
   ]
  }
 ],
 "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.8.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}