Heyde-et-al-2021-Python.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Python code for:\n",
    "# Alexander Heyde, David Rohde, Cameron S. McAlpine, Shuang Zhang, Friedrich F. Hoyer, Jeffrey M. Gerold, David Cheek, Yoshiko Iwamoto, Maximilian J. Schloss, Katrien Vandoorne, Oriol Iborra - Egea, Christian Muñoz - Guijosa, Antoni Bayes - Genis, Johannes G. Reiter, Morgan Craig, Filip K. Swirski, Matthias Nahrendorf, Martin A. Nowak, Kamila Naxerova. \n",
    "# Increased stem cell proliferation in atherosclerosis accelerates clonal hematopoiesis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np                    # mathematics\n",
    "import kingman                        # coalescence\n",
    "import seaborn as sns                 # visualization\n",
    "from matplotlib import pyplot as plt  # plotting\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Table S2\n",
    "# All time measurements in years\n",
    "b = 365/28     # Average baseline HSC proliferation rate\n",
    "age = 40       # Age of onset for elevated HSC proliferation rate\n",
    "T = 70         # Age at sequencing\n",
    "thresh1 = .02  # Detection frequency (baseline)\n",
    "thresh2 = .08  # Detection frequency (alternative)\n",
    "pop = 10000    # Number of HSCs\n",
    "mutrate = 1.5  # Mean number of neutral mutations acquired by each daughter cell per division\n",
    "\n",
    "M = 1000       # Number of independent realizations\n",
    "body = 10**13  # Number of cells in the body after development\n",
    "progen = 100   # Number of HSC progenitor cells\n",
    "u = 50         # Possible mutational realizations before birth\n",
    "U = 50         # Possible mutational realizations after birth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "999\n",
      "998\n",
      "997\n",
      "996\n",
      "995\n",
      "994\n",
      "993\n",
      "992\n",
      "991\n",
      "990\n",
      "989\n",
      "988\n",
      "987\n",
      "986\n",
      "985\n",
      "984\n",
      "983\n",
      "982\n",
      "981\n",
      "980\n",
      "979\n",
      "978\n",
      "977\n",
      "976\n",
      "975\n",
      "974\n",
      "973\n",
      "972\n",
      "971\n",
      "970\n",
      "969\n",
      "968\n",
      "967\n",
      "966\n",
      "965\n",
      "964\n",
      "963\n",
      "962\n",
      "961\n",
      "960\n",
      "959\n",
      "958\n",
      "957\n",
      "956\n",
      "955\n",
      "954\n",
      "953\n",
      "952\n",
      "951\n",
      "950\n",
      "949\n",
      "948\n",
      "947\n",
      "946\n",
      "945\n",
      "944\n",
      "943\n",
      "942\n",
      "941\n",
      "940\n",
      "939\n",
      "938\n",
      "937\n",
      "936\n",
      "935\n",
      "934\n",
      "933\n",
      "932\n",
      "931\n",
      "930\n",
      "929\n",
      "928\n",
      "927\n",
      "926\n",
      "925\n",
      "924\n",
      "923\n",
      "922\n",
      "921\n",
      "920\n",
      "919\n",
      "918\n",
      "917\n",
      "916\n",
      "915\n",
      "914\n",
      "913\n",
      "912\n",
      "911\n",
      "910\n",
      "909\n",
      "908\n",
      "907\n",
      "906\n",
      "905\n",
      "904\n",
      "903\n",
      "902\n",
      "901\n",
      "900\n",
      "899\n",
      "898\n",
      "897\n",
      "896\n",
      "895\n",
      "894\n",
      "893\n",
      "892\n",
      "891\n",
      "890\n",
      "889\n",
      "888\n",
      "887\n",
      "886\n",
      "885\n",
      "884\n",
      "883\n",
      "882\n",
      "881\n",
      "880\n",
      "879\n",
      "878\n",
      "877\n",
      "876\n",
      "875\n",
      "874\n",
      "873\n",
      "872\n",
      "871\n",
      "870\n",
      "869\n",
      "868\n",
      "867\n",
      "866\n",
      "865\n",
      "864\n",
      "863\n",
      "862\n",
      "861\n",
      "860\n",
      "859\n",
      "858\n",
      "857\n",
      "856\n",
      "855\n",
      "854\n",
      "853\n",
      "852\n",
      "851\n",
      "850\n",
      "849\n",
      "848\n",
      "847\n",
      "846\n",
      "845\n",
      "844\n",
      "843\n",
      "842\n",
      "841\n",
      "840\n",
      "839\n",
      "838\n",
      "837\n",
      "836\n",
      "835\n",
      "834\n",
      "833\n",
      "832\n",
      "831\n",
      "830\n",
      "829\n",
      "828\n",
      "827\n",
      "826\n",
      "825\n",
      "824\n",
      "823\n",
      "822\n",
      "821\n",
      "820\n",
      "819\n",
      "818\n",
      "817\n",
      "816\n",
      "815\n",
      "814\n",
      "813\n",
      "812\n",
      "811\n",
      "810\n",
      "809\n",
      "808\n",
      "807\n",
      "806\n",
      "805\n",
      "804\n",
      "803\n",
      "802\n",
      "801\n",
      "800\n",
      "799\n",
      "798\n",
      "797\n",
      "796\n",
      "795\n",
      "794\n",
      "793\n",
      "792\n",
      "791\n",
      "790\n",
      "789\n",
      "788\n",
      "787\n",
      "786\n",
      "785\n",
      "784\n",
      "783\n",
      "782\n",
      "781\n",
      "780\n",
      "779\n",
      "778\n",
      "777\n",
      "776\n",
      "775\n",
      "774\n",
      "773\n",
      "772\n",
      "771\n",
      "770\n",
      "769\n",
      "768\n",
      "767\n",
      "766\n",
      "765\n",
      "764\n",
      "763\n",
      "762\n",
      "761\n",
      "760\n",
      "759\n",
      "758\n",
      "757\n",
      "756\n",
      "755\n",
      "754\n",
      "753\n",
      "752\n",
      "751\n",
      "750\n",
      "749\n",
      "748\n",
      "747\n",
      "746\n",
      "745\n",
      "744\n",
      "743\n",
      "742\n",
      "741\n",
      "740\n",
      "739\n",
      "738\n",
      "737\n",
      "736\n",
      "735\n",
      "734\n",
      "733\n",
      "732\n",
      "731\n",
      "730\n",
      "729\n",
      "728\n",
      "727\n",
      "726\n",
      "725\n",
      "724\n",
      "723\n",
      "722\n",
      "721\n",
      "720\n",
      "719\n",
      "718\n",
      "717\n",
      "716\n",
      "715\n",
      "714\n",
      "713\n",
      "712\n",
      "711\n",
      "710\n",
      "709\n",
      "708\n",
      "707\n",
      "706\n",
      "705\n",
      "704\n",
      "703\n",
      "702\n",
      "701\n",
      "700\n",
      "699\n",
      "698\n",
      "697\n",
      "696\n",
      "695\n",
      "694\n",
      "693\n",
      "692\n",
      "691\n",
      "690\n",
      "689\n",
      "688\n",
      "687\n",
      "686\n",
      "685\n",
      "684\n",
      "683\n",
      "682\n",
      "681\n",
      "680\n",
      "679\n",
      "678\n",
      "677\n",
      "676\n",
      "675\n",
      "674\n",
      "673\n",
      "672\n",
      "671\n",
      "670\n",
      "669\n",
      "668\n",
      "667\n",
      "666\n",
      "665\n",
      "664\n",
      "663\n",
      "662\n",
      "661\n",
      "660\n",
      "659\n",
      "658\n",
      "657\n",
      "656\n",
      "655\n",
      "654\n",
      "653\n",
      "652\n",
      "651\n",
      "650\n",
      "649\n",
      "648\n",
      "647\n",
      "646\n",
      "645\n",
      "644\n",
      "643\n",
      "642\n",
      "641\n",
      "640\n",
      "639\n",
      "638\n",
      "637\n",
      "636\n",
      "635\n",
      "634\n",
      "633\n",
      "632\n",
      "631\n",
      "630\n",
      "629\n",
      "628\n",
      "627\n",
      "626\n",
      "625\n",
      "624\n",
      "623\n",
      "622\n",
      "621\n",
      "620\n",
      "619\n",
      "618\n",
      "617\n",
      "616\n",
      "615\n",
      "614\n",
      "613\n",
      "612\n",
      "611\n",
      "610\n",
      "609\n",
      "608\n",
      "607\n",
      "606\n",
      "605\n",
      "604\n",
      "603\n",
      "602\n",
      "601\n",
      "600\n",
      "599\n",
      "598\n",
      "597\n",
      "596\n",
      "595\n",
      "594\n",
      "593\n",
      "592\n",
      "591\n",
      "590\n",
      "589\n",
      "588\n",
      "587\n",
      "586\n",
      "585\n",
      "584\n",
      "583\n",
      "582\n",
      "581\n",
      "580\n",
      "579\n",
      "578\n",
      "577\n",
      "576\n",
      "575\n",
      "574\n",
      "573\n",
      "572\n",
      "571\n",
      "570\n",
      "569\n",
      "568\n",
      "567\n",
      "566\n",
      "565\n",
      "564\n",
      "563\n",
      "562\n",
      "561\n",
      "560\n",
      "559\n",
      "558\n",
      "557\n",
      "556\n",
      "555\n",
      "554\n",
      "553\n",
      "552\n",
      "551\n",
      "550\n",
      "549\n",
      "548\n",
      "547\n",
      "546\n",
      "545\n",
      "544\n",
      "543\n",
      "542\n",
      "541\n",
      "540\n",
      "539\n",
      "538\n",
      "537\n",
      "536\n",
      "535\n",
      "534\n",
      "533\n",
      "532\n",
      "531\n",
      "530\n",
      "529\n",
      "528\n",
      "527\n",
      "526\n",
      "525\n",
      "524\n",
      "523\n",
      "522\n",
      "521\n",
      "520\n",
      "519\n",
      "518\n",
      "517\n",
      "516\n",
      "515\n",
      "514\n",
      "513\n",
      "512\n",
      "511\n",
      "510\n",
      "509\n",
      "508\n",
      "507\n",
      "506\n",
      "505\n",
      "504\n",
      "503\n",
      "502\n",
      "501\n",
      "500\n",
      "499\n",
      "498\n",
      "497\n",
      "496\n",
      "495\n",
      "494\n",
      "493\n",
      "492\n",
      "491\n",
      "490\n",
      "489\n",
      "488\n",
      "487\n",
      "486\n",
      "485\n",
      "484\n",
      "483\n",
      "482\n",
      "481\n",
      "480\n",
      "479\n",
      "478\n",
      "477\n",
      "476\n",
      "475\n",
      "474\n",
      "473\n",
      "472\n",
      "471\n",
      "470\n",
      "469\n",
      "468\n",
      "467\n",
      "466\n",
      "465\n",
      "464\n",
      "463\n",
      "462\n",
      "461\n",
      "460\n",
      "459\n",
      "458\n",
      "457\n",
      "456\n",
      "455\n",
      "454\n",
      "453\n",
      "452\n",
      "451\n",
      "450\n",
      "449\n",
      "448\n",
      "447\n",
      "446\n",
      "445\n",
      "444\n",
      "443\n",
      "442\n",
      "441\n",
      "440\n",
      "439\n",
      "438\n",
      "437\n",
      "436\n",
      "435\n",
      "434\n",
      "433\n",
      "432\n",
      "431\n",
      "430\n",
      "429\n",
      "428\n",
      "427\n",
      "426\n",
      "425\n",
      "424\n",
      "423\n",
      "422\n",
      "421\n",
      "420\n",
      "419\n",
      "418\n",
      "417\n",
      "416\n",
      "415\n",
      "414\n",
      "413\n",
      "412\n",
      "411\n",
      "410\n",
      "409\n",
      "408\n",
      "407\n",
      "406\n",
      "405\n",
      "404\n",
      "403\n",
      "402\n",
      "401\n",
      "400\n",
      "399\n",
      "398\n",
      "397\n",
      "396\n",
      "395\n",
      "394\n",
      "393\n",
      "392\n",
      "391\n",
      "390\n",
      "389\n",
      "388\n",
      "387\n",
      "386\n",
      "385\n",
      "384\n",
      "383\n",
      "382\n",
      "381\n",
      "380\n",
      "379\n",
      "378\n",
      "377\n",
      "376\n",
      "375\n",
      "374\n",
      "373\n",
      "372\n",
      "371\n",
      "370\n",
      "369\n",
      "368\n",
      "367\n",
      "366\n",
      "365\n",
      "364\n",
      "363\n",
      "362\n",
      "361\n",
      "360\n",
      "359\n",
      "358\n",
      "357\n",
      "356\n",
      "355\n",
      "354\n",
      "353\n",
      "352\n",
      "351\n",
      "350\n",
      "349\n",
      "348\n",
      "347\n",
      "346\n",
      "345\n",
      "344\n",
      "343\n",
      "342\n",
      "341\n",
      "340\n",
      "339\n",
      "338\n",
      "337\n",
      "336\n",
      "335\n",
      "334\n",
      "333\n",
      "332\n",
      "331\n",
      "330\n",
      "329\n",
      "328\n",
      "327\n",
      "326\n",
      "325\n",
      "324\n",
      "323\n",
      "322\n",
      "321\n",
      "320\n",
      "319\n",
      "318\n",
      "317\n",
      "316\n",
      "315\n",
      "314\n",
      "313\n",
      "312\n",
      "311\n",
      "310\n",
      "309\n",
      "308\n",
      "307\n",
      "306\n",
      "305\n",
      "304\n",
      "303\n",
      "302\n",
      "301\n",
      "300\n",
      "299\n",
      "298\n",
      "297\n",
      "296\n",
      "295\n",
      "294\n",
      "293\n",
      "292\n",
      "291\n",
      "290\n",
      "289\n",
      "288\n",
      "287\n",
      "286\n",
      "285\n",
      "284\n",
      "283\n",
      "282\n",
      "281\n",
      "280\n",
      "279\n",
      "278\n",
      "277\n",
      "276\n",
      "275\n",
      "274\n",
      "273\n",
      "272\n",
      "271\n",
      "270\n",
      "269\n",
      "268\n",
      "267\n",
      "266\n",
      "265\n",
      "264\n",
      "263\n",
      "262\n",
      "261\n",
      "260\n",
      "259\n",
      "258\n",
      "257\n",
      "256\n",
      "255\n",
      "254\n",
      "253\n",
      "252\n",
      "251\n",
      "250\n",
      "249\n",
      "248\n",
      "247\n",
      "246\n",
      "245\n",
      "244\n",
      "243\n",
      "242\n",
      "241\n",
      "240\n",
      "239\n",
      "238\n",
      "237\n",
      "236\n",
      "235\n",
      "234\n",
      "233\n",
      "232\n",
      "231\n",
      "230\n",
      "229\n",
      "228\n",
      "227\n",
      "226\n",
      "225\n",
      "224\n",
      "223\n",
      "222\n",
      "221\n",
      "220\n",
      "219\n",
      "218\n",
      "217\n",
      "216\n",
      "215\n",
      "214\n",
      "213\n",
      "212\n",
      "211\n",
      "210\n",
      "209\n",
      "208\n",
      "207\n",
      "206\n",
      "205\n",
      "204\n",
      "203\n",
      "202\n",
      "201\n",
      "200\n",
      "199\n",
      "198\n",
      "197\n",
      "196\n",
      "195\n",
      "194\n",
      "193\n",
      "192\n",
      "191\n",
      "190\n",
      "189\n",
      "188\n",
      "187\n",
      "186\n",
      "185\n",
      "184\n",
      "183\n",
      "182\n",
      "181\n",
      "180\n",
      "179\n",
      "178\n",
      "177\n",
      "176\n",
      "175\n",
      "174\n",
      "173\n",
      "172\n",
      "171\n",
      "170\n",
      "169\n",
      "168\n",
      "167\n",
      "166\n",
      "165\n",
      "164\n",
      "163\n",
      "162\n",
      "161\n",
      "160\n",
      "159\n",
      "158\n",
      "157\n",
      "156\n",
      "155\n",
      "154\n",
      "153\n",
      "152\n",
      "151\n",
      "150\n",
      "149\n",
      "148\n",
      "147\n",
      "146\n",
      "145\n",
      "144\n",
      "143\n",
      "142\n",
      "141\n",
      "140\n",
      "139\n",
      "138\n",
      "137\n",
      "136\n",
      "135\n",
      "134\n",
      "133\n",
      "132\n",
      "131\n",
      "130\n",
      "129\n",
      "128\n",
      "127\n",
      "126\n",
      "125\n",
      "124\n",
      "123\n",
      "122\n",
      "121\n",
      "120\n",
      "119\n",
      "118\n",
      "117\n",
      "116\n",
      "115\n",
      "114\n",
      "113\n",
      "112\n",
      "111\n",
      "110\n",
      "109\n",
      "108\n",
      "107\n",
      "106\n",
      "105\n",
      "104\n",
      "103\n",
      "102\n",
      "101\n",
      "100\n",
      "99\n",
      "98\n",
      "97\n",
      "96\n",
      "95\n",
      "94\n",
      "93\n",
      "92\n",
      "91\n",
      "90\n",
      "89\n",
      "88\n",
      "87\n",
      "86\n",
      "85\n",
      "84\n",
      "83\n",
      "82\n",
      "81\n",
      "80\n",
      "79\n",
      "78\n",
      "77\n",
      "76\n",
      "75\n",
      "74\n",
      "73\n",
      "72\n",
      "71\n",
      "70\n",
      "69\n",
      "68\n",
      "67\n",
      "66\n",
      "65\n",
      "64\n",
      "63\n",
      "62\n",
      "61\n",
      "60\n",
      "59\n",
      "58\n",
      "57\n",
      "56\n",
      "55\n",
      "54\n",
      "53\n",
      "52\n",
      "51\n",
      "50\n",
      "49\n",
      "48\n",
      "47\n",
      "46\n",
      "45\n",
      "44\n",
      "43\n",
      "42\n",
      "41\n",
      "40\n",
      "39\n",
      "38\n",
      "37\n",
      "36\n",
      "35\n",
      "34\n",
      "33\n",
      "32\n",
      "31\n",
      "30\n",
      "29\n",
      "28\n",
      "27\n",
      "26\n",
      "25\n",
      "24\n",
      "23\n",
      "22\n",
      "21\n",
      "20\n",
      "19\n",
      "18\n",
      "17\n",
      "16\n",
      "15\n",
      "14\n",
      "13\n",
      "12\n",
      "11\n",
      "10\n",
      "9\n",
      "8\n",
      "7\n",
      "6\n",
      "5\n",
      "4\n",
      "3\n",
      "2\n",
      "1\n",
      "[  1.00010000e+00   1.00000000e-04   1.00000000e-04 ...,   7.04400000e-01\n",
      "   9.85100000e-01   1.00000000e+00]\n",
      "Wall time: 24min 19s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# Control cohort: neutral VAF simulations\n",
    "\n",
    "# Set distribution of proliferation rates for control cohort\n",
    "R = np.array([0.4883, 1.8066, 0.9861, 0.9861, 0.9151, 1.0966, 0.2832, 0.6453, 1.9959, 0.7968]) # from Fig 2C\n",
    "lifes = (min(T,age)+max(0,T-age)*R)*b\n",
    "x = np.array([])\n",
    "weightx = np.array([])\n",
    "\n",
    "# Simulate mutation and cell division in control patients\n",
    "nx0 = [0]*M\n",
    "nx1 = [0]*M\n",
    "nx2 = [0]*M\n",
    "for m in range(M):\n",
    "    print(M-m)\n",
    "    \n",
    "    #Simulate cell division backwards in time via coalescent\n",
    "    #np.random.seed(m)\n",
    "    life = np.random.choice(lifes)\n",
    "    outp = kingman.simulate(pop,m)\n",
    "    parents = outp[0]\n",
    "    times = np.array(outp[1])*2*pop\n",
    "    times[0] = 0\n",
    "    finalStage = life+np.log(pop/progen)+np.log(1+times*body/pop-life-np.log(pop/progen))\n",
    "    times = np.where(times < life, times, np.where(times < life+pop/progen-1, life+np.log(1+times-life), finalStage))\n",
    "    nodes = len(times)\n",
    "\n",
    "    #Determine coalescent tree structure and identify mutations\n",
    "    branch = [0]*nodes\n",
    "    branchearly = [0]*nodes\n",
    "    for n in range(1,nodes):\n",
    "        branch[n] = np.maximum(0,np.minimum(times[parents[n]],life)-times[n])\n",
    "        branchearly[n] = np.maximum(0,times[parents[n]]-np.maximum(times[n],life))\n",
    "    branch[nodes-1] = 0\n",
    "    branchearly[nodes-1] = max(0,np.log(pop/progen*body)+life-times[nodes-1])\n",
    "    #np.random.seed(m+M)\n",
    "    muts = np.random.poisson(np.array(branch)*U+np.array(branchearly)*u)\n",
    "\n",
    "    #Trace mutations down ancestry to calculate variant frequency\n",
    "    cells = [0]*nodes\n",
    "    for n in [i for i in range(len(times)) if times[i] == 0]:\n",
    "        while n > -1:\n",
    "            cells[n] = cells[n]+1\n",
    "            n = parents[n]\n",
    "    freq = np.array(cells)/pop\n",
    "    \n",
    "    #Count mutations that exceed a VAF threshold of 0, thresh1, or thresh2\n",
    "    x = np.concatenate([x, freq])\n",
    "    weightx = np.concatenate([weightx, muts])\n",
    "    nx0[m] = np.dot(freq>0, muts)\n",
    "    nx1[m] = np.dot(freq>thresh1, muts)\n",
    "    nx2[m] = np.dot(freq>thresh2, muts)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "999\n",
      "998\n",
      "997\n",
      "996\n",
      "995\n",
      "994\n",
      "993\n",
      "992\n",
      "991\n",
      "990\n",
      "989\n",
      "988\n",
      "987\n",
      "986\n",
      "985\n",
      "984\n",
      "983\n",
      "982\n",
      "981\n",
      "980\n",
      "979\n",
      "978\n",
      "977\n",
      "976\n",
      "975\n",
      "974\n",
      "973\n",
      "972\n",
      "971\n",
      "970\n",
      "969\n",
      "968\n",
      "967\n",
      "966\n",
      "965\n",
      "964\n",
      "963\n",
      "962\n",
      "961\n",
      "960\n",
      "959\n",
      "958\n",
      "957\n",
      "956\n",
      "955\n",
      "954\n",
      "953\n",
      "952\n",
      "951\n",
      "950\n",
      "949\n",
      "948\n",
      "947\n",
      "946\n",
      "945\n",
      "944\n",
      "943\n",
      "942\n",
      "941\n",
      "940\n",
      "939\n",
      "938\n",
      "937\n",
      "936\n",
      "935\n",
      "934\n",
      "933\n",
      "932\n",
      "931\n",
      "930\n",
      "929\n",
      "928\n",
      "927\n",
      "926\n",
      "925\n",
      "924\n",
      "923\n",
      "922\n",
      "921\n",
      "920\n",
      "919\n",
      "918\n",
      "917\n",
      "916\n",
      "915\n",
      "914\n",
      "913\n",
      "912\n",
      "911\n",
      "910\n",
      "909\n",
      "908\n",
      "907\n",
      "906\n",
      "905\n",
      "904\n",
      "903\n",
      "902\n",
      "901\n",
      "900\n",
      "899\n",
      "898\n",
      "897\n",
      "896\n",
      "895\n",
      "894\n",
      "893\n",
      "892\n",
      "891\n",
      "890\n",
      "889\n",
      "888\n",
      "887\n",
      "886\n",
      "885\n",
      "884\n",
      "883\n",
      "882\n",
      "881\n",
      "880\n",
      "879\n",
      "878\n",
      "877\n",
      "876\n",
      "875\n",
      "874\n",
      "873\n",
      "872\n",
      "871\n",
      "870\n",
      "869\n",
      "868\n",
      "867\n",
      "866\n",
      "865\n",
      "864\n",
      "863\n",
      "862\n",
      "861\n",
      "860\n",
      "859\n",
      "858\n",
      "857\n",
      "856\n",
      "855\n",
      "854\n",
      "853\n",
      "852\n",
      "851\n",
      "850\n",
      "849\n",
      "848\n",
      "847\n",
      "846\n",
      "845\n",
      "844\n",
      "843\n",
      "842\n",
      "841\n",
      "840\n",
      "839\n",
      "838\n",
      "837\n",
      "836\n",
      "835\n",
      "834\n",
      "833\n",
      "832\n",
      "831\n",
      "830\n",
      "829\n",
      "828\n",
      "827\n",
      "826\n",
      "825\n",
      "824\n",
      "823\n",
      "822\n",
      "821\n",
      "820\n",
      "819\n",
      "818\n",
      "817\n",
      "816\n",
      "815\n",
      "814\n",
      "813\n",
      "812\n",
      "811\n",
      "810\n",
      "809\n",
      "808\n",
      "807\n",
      "806\n",
      "805\n",
      "804\n",
      "803\n",
      "802\n",
      "801\n",
      "800\n",
      "799\n",
      "798\n",
      "797\n",
      "796\n",
      "795\n",
      "794\n",
      "793\n",
      "792\n",
      "791\n",
      "790\n",
      "789\n",
      "788\n",
      "787\n",
      "786\n",
      "785\n",
      "784\n",
      "783\n",
      "782\n",
      "781\n",
      "780\n",
      "779\n",
      "778\n",
      "777\n",
      "776\n",
      "775\n",
      "774\n",
      "773\n",
      "772\n",
      "771\n",
      "770\n",
      "769\n",
      "768\n",
      "767\n",
      "766\n",
      "765\n",
      "764\n",
      "763\n",
      "762\n",
      "761\n",
      "760\n",
      "759\n",
      "758\n",
      "757\n",
      "756\n",
      "755\n",
      "754\n",
      "753\n",
      "752\n",
      "751\n",
      "750\n",
      "749\n",
      "748\n",
      "747\n",
      "746\n",
      "745\n",
      "744\n",
      "743\n",
      "742\n",
      "741\n",
      "740\n",
      "739\n",
      "738\n",
      "737\n",
      "736\n",
      "735\n",
      "734\n",
      "733\n",
      "732\n",
      "731\n",
      "730\n",
      "729\n",
      "728\n",
      "727\n",
      "726\n",
      "725\n",
      "724\n",
      "723\n",
      "722\n",
      "721\n",
      "720\n",
      "719\n",
      "718\n",
      "717\n",
      "716\n",
      "715\n",
      "714\n",
      "713\n",
      "712\n",
      "711\n",
      "710\n",
      "709\n",
      "708\n",
      "707\n",
      "706\n",
      "705\n",
      "704\n",
      "703\n",
      "702\n",
      "701\n",
      "700\n",
      "699\n",
      "698\n",
      "697\n",
      "696\n",
      "695\n",
      "694\n",
      "693\n",
      "692\n",
      "691\n",
      "690\n",
      "689\n",
      "688\n",
      "687\n",
      "686\n",
      "685\n",
      "684\n",
      "683\n",
      "682\n",
      "681\n",
      "680\n",
      "679\n",
      "678\n",
      "677\n",
      "676\n",
      "675\n",
      "674\n",
      "673\n",
      "672\n",
      "671\n",
      "670\n",
      "669\n",
      "668\n",
      "667\n",
      "666\n",
      "665\n",
      "664\n",
      "663\n",
      "662\n",
      "661\n",
      "660\n",
      "659\n",
      "658\n",
      "657\n",
      "656\n",
      "655\n",
      "654\n",
      "653\n",
      "652\n",
      "651\n",
      "650\n",
      "649\n",
      "648\n",
      "647\n",
      "646\n",
      "645\n",
      "644\n",
      "643\n",
      "642\n",
      "641\n",
      "640\n",
      "639\n",
      "638\n",
      "637\n",
      "636\n",
      "635\n",
      "634\n",
      "633\n",
      "632\n",
      "631\n",
      "630\n",
      "629\n",
      "628\n",
      "627\n",
      "626\n",
      "625\n",
      "624\n",
      "623\n",
      "622\n",
      "621\n",
      "620\n",
      "619\n",
      "618\n",
      "617\n",
      "616\n",
      "615\n",
      "614\n",
      "613\n",
      "612\n",
      "611\n",
      "610\n",
      "609\n",
      "608\n",
      "607\n",
      "606\n",
      "605\n",
      "604\n",
      "603\n",
      "602\n",
      "601\n",
      "600\n",
      "599\n",
      "598\n",
      "597\n",
      "596\n",
      "595\n",
      "594\n",
      "593\n",
      "592\n",
      "591\n",
      "590\n",
      "589\n",
      "588\n",
      "587\n",
      "586\n",
      "585\n",
      "584\n",
      "583\n",
      "582\n",
      "581\n",
      "580\n",
      "579\n",
      "578\n",
      "577\n",
      "576\n",
      "575\n",
      "574\n",
      "573\n",
      "572\n",
      "571\n",
      "570\n",
      "569\n",
      "568\n",
      "567\n",
      "566\n",
      "565\n",
      "564\n",
      "563\n",
      "562\n",
      "561\n",
      "560\n",
      "559\n",
      "558\n",
      "557\n",
      "556\n",
      "555\n",
      "554\n",
      "553\n",
      "552\n",
      "551\n",
      "550\n",
      "549\n",
      "548\n",
      "547\n",
      "546\n",
      "545\n",
      "544\n",
      "543\n",
      "542\n",
      "541\n",
      "540\n",
      "539\n",
      "538\n",
      "537\n",
      "536\n",
      "535\n",
      "534\n",
      "533\n",
      "532\n",
      "531\n",
      "530\n",
      "529\n",
      "528\n",
      "527\n",
      "526\n",
      "525\n",
      "524\n",
      "523\n",
      "522\n",
      "521\n",
      "520\n",
      "519\n",
      "518\n",
      "517\n",
      "516\n",
      "515\n",
      "514\n",
      "513\n",
      "512\n",
      "511\n",
      "510\n",
      "509\n",
      "508\n",
      "507\n",
      "506\n",
      "505\n",
      "504\n",
      "503\n",
      "502\n",
      "501\n",
      "500\n",
      "499\n",
      "498\n",
      "497\n",
      "496\n",
      "495\n",
      "494\n",
      "493\n",
      "492\n",
      "491\n",
      "490\n",
      "489\n",
      "488\n",
      "487\n",
      "486\n",
      "485\n",
      "484\n",
      "483\n",
      "482\n",
      "481\n",
      "480\n",
      "479\n",
      "478\n",
      "477\n",
      "476\n",
      "475\n",
      "474\n",
      "473\n",
      "472\n",
      "471\n",
      "470\n",
      "469\n",
      "468\n",
      "467\n",
      "466\n",
      "465\n",
      "464\n",
      "463\n",
      "462\n",
      "461\n",
      "460\n",
      "459\n",
      "458\n",
      "457\n",
      "456\n",
      "455\n",
      "454\n",
      "453\n",
      "452\n",
      "451\n",
      "450\n",
      "449\n",
      "448\n",
      "447\n",
      "446\n",
      "445\n",
      "444\n",
      "443\n",
      "442\n",
      "441\n",
      "440\n",
      "439\n",
      "438\n",
      "437\n",
      "436\n",
      "435\n",
      "434\n",
      "433\n",
      "432\n",
      "431\n",
      "430\n",
      "429\n",
      "428\n",
      "427\n",
      "426\n",
      "425\n",
      "424\n",
      "423\n",
      "422\n",
      "421\n",
      "420\n",
      "419\n",
      "418\n",
      "417\n",
      "416\n",
      "415\n",
      "414\n",
      "413\n",
      "412\n",
      "411\n",
      "410\n",
      "409\n",
      "408\n",
      "407\n",
      "406\n",
      "405\n",
      "404\n",
      "403\n",
      "402\n",
      "401\n",
      "400\n",
      "399\n",
      "398\n",
      "397\n",
      "396\n",
      "395\n",
      "394\n",
      "393\n",
      "392\n",
      "391\n",
      "390\n",
      "389\n",
      "388\n",
      "387\n",
      "386\n",
      "385\n",
      "384\n",
      "383\n",
      "382\n",
      "381\n",
      "380\n",
      "379\n",
      "378\n",
      "377\n",
      "376\n",
      "375\n",
      "374\n",
      "373\n",
      "372\n",
      "371\n",
      "370\n",
      "369\n",
      "368\n",
      "367\n",
      "366\n",
      "365\n",
      "364\n",
      "363\n",
      "362\n",
      "361\n",
      "360\n",
      "359\n",
      "358\n",
      "357\n",
      "356\n",
      "355\n",
      "354\n",
      "353\n",
      "352\n",
      "351\n",
      "350\n",
      "349\n",
      "348\n",
      "347\n",
      "346\n",
      "345\n",
      "344\n",
      "343\n",
      "342\n",
      "341\n",
      "340\n",
      "339\n",
      "338\n",
      "337\n",
      "336\n",
      "335\n",
      "334\n",
      "333\n",
      "332\n",
      "331\n",
      "330\n",
      "329\n",
      "328\n",
      "327\n",
      "326\n",
      "325\n",
      "324\n",
      "323\n",
      "322\n",
      "321\n",
      "320\n",
      "319\n",
      "318\n",
      "317\n",
      "316\n",
      "315\n",
      "314\n",
      "313\n",
      "312\n",
      "311\n",
      "310\n",
      "309\n",
      "308\n",
      "307\n",
      "306\n",
      "305\n",
      "304\n",
      "303\n",
      "302\n",
      "301\n",
      "300\n",
      "299\n",
      "298\n",
      "297\n",
      "296\n",
      "295\n",
      "294\n",
      "293\n",
      "292\n",
      "291\n",
      "290\n",
      "289\n",
      "288\n",
      "287\n",
      "286\n",
      "285\n",
      "284\n",
      "283\n",
      "282\n",
      "281\n",
      "280\n",
      "279\n",
      "278\n",
      "277\n",
      "276\n",
      "275\n",
      "274\n",
      "273\n",
      "272\n",
      "271\n",
      "270\n",
      "269\n",
      "268\n",
      "267\n",
      "266\n",
      "265\n",
      "264\n",
      "263\n",
      "262\n",
      "261\n",
      "260\n",
      "259\n",
      "258\n",
      "257\n",
      "256\n",
      "255\n",
      "254\n",
      "253\n",
      "252\n",
      "251\n",
      "250\n",
      "249\n",
      "248\n",
      "247\n",
      "246\n",
      "245\n",
      "244\n",
      "243\n",
      "242\n",
      "241\n",
      "240\n",
      "239\n",
      "238\n",
      "237\n",
      "236\n",
      "235\n",
      "234\n",
      "233\n",
      "232\n",
      "231\n",
      "230\n",
      "229\n",
      "228\n",
      "227\n",
      "226\n",
      "225\n",
      "224\n",
      "223\n",
      "222\n",
      "221\n",
      "220\n",
      "219\n",
      "218\n",
      "217\n",
      "216\n",
      "215\n",
      "214\n",
      "213\n",
      "212\n",
      "211\n",
      "210\n",
      "209\n",
      "208\n",
      "207\n",
      "206\n",
      "205\n",
      "204\n",
      "203\n",
      "202\n",
      "201\n",
      "200\n",
      "199\n",
      "198\n",
      "197\n",
      "196\n",
      "195\n",
      "194\n",
      "193\n",
      "192\n",
      "191\n",
      "190\n",
      "189\n",
      "188\n",
      "187\n",
      "186\n",
      "185\n",
      "184\n",
      "183\n",
      "182\n",
      "181\n",
      "180\n",
      "179\n",
      "178\n",
      "177\n",
      "176\n",
      "175\n",
      "174\n",
      "173\n",
      "172\n",
      "171\n",
      "170\n",
      "169\n",
      "168\n",
      "167\n",
      "166\n",
      "165\n",
      "164\n",
      "163\n",
      "162\n",
      "161\n",
      "160\n",
      "159\n",
      "158\n",
      "157\n",
      "156\n",
      "155\n",
      "154\n",
      "153\n",
      "152\n",
      "151\n",
      "150\n",
      "149\n",
      "148\n",
      "147\n",
      "146\n",
      "145\n",
      "144\n",
      "143\n",
      "142\n",
      "141\n",
      "140\n",
      "139\n",
      "138\n",
      "137\n",
      "136\n",
      "135\n",
      "134\n",
      "133\n",
      "132\n",
      "131\n",
      "130\n",
      "129\n",
      "128\n",
      "127\n",
      "126\n",
      "125\n",
      "124\n",
      "123\n",
      "122\n",
      "121\n",
      "120\n",
      "119\n",
      "118\n",
      "117\n",
      "116\n",
      "115\n",
      "114\n",
      "113\n",
      "112\n",
      "111\n",
      "110\n",
      "109\n",
      "108\n",
      "107\n",
      "106\n",
      "105\n",
      "104\n",
      "103\n",
      "102\n",
      "101\n",
      "100\n",
      "99\n",
      "98\n",
      "97\n",
      "96\n",
      "95\n",
      "94\n",
      "93\n",
      "92\n",
      "91\n",
      "90\n",
      "89\n",
      "88\n",
      "87\n",
      "86\n",
      "85\n",
      "84\n",
      "83\n",
      "82\n",
      "81\n",
      "80\n",
      "79\n",
      "78\n",
      "77\n",
      "76\n",
      "75\n",
      "74\n",
      "73\n",
      "72\n",
      "71\n",
      "70\n",
      "69\n",
      "68\n",
      "67\n",
      "66\n",
      "65\n",
      "64\n",
      "63\n",
      "62\n",
      "61\n",
      "60\n",
      "59\n",
      "58\n",
      "57\n",
      "56\n",
      "55\n",
      "54\n",
      "53\n",
      "52\n",
      "51\n",
      "50\n",
      "49\n",
      "48\n",
      "47\n",
      "46\n",
      "45\n",
      "44\n",
      "43\n",
      "42\n",
      "41\n",
      "40\n",
      "39\n",
      "38\n",
      "37\n",
      "36\n",
      "35\n",
      "34\n",
      "33\n",
      "32\n",
      "31\n",
      "30\n",
      "29\n",
      "28\n",
      "27\n",
      "26\n",
      "25\n",
      "24\n",
      "23\n",
      "22\n",
      "21\n",
      "20\n",
      "19\n",
      "18\n",
      "17\n",
      "16\n",
      "15\n",
      "14\n",
      "13\n",
      "12\n",
      "11\n",
      "10\n",
      "9\n",
      "8\n",
      "7\n",
      "6\n",
      "5\n",
      "4\n",
      "3\n",
      "2\n",
      "1\n",
      "[  1.00010000e+00   1.00000000e-04   1.00000000e-04 ...,   7.04400000e-01\n",
      "   9.85100000e-01   1.00000000e+00]\n",
      "Wall time: 24min 13s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# Disease cohort: neutral VAF simulations\n",
    "\n",
    "# Set distribution of proliferation rates for disease cohort\n",
    "R = np.array([1.9328, 2.2878, 1.6961, 3.8025, 1.4595, 3.3843, 1.7908, 3.7393, 1.9170, 2.8637]) # from Fig 2C\n",
    "lifes = (min(T,age)+max(0,T-age)*R)*b\n",
    "y = np.array([])\n",
    "weighty = np.array([])\n",
    "\n",
    "# Simulate mutation and cell division in disease patients\n",
    "ny0 = [0]*M\n",
    "ny1 = [0]*M\n",
    "ny2 = [0]*M\n",
    "for m in range(M):\n",
    "    print(M-m)\n",
    "    \n",
    "    #Simulate cell division backwards in time via coalescent\n",
    "    #np.random.seed(m)\n",
    "    life = np.random.choice(lifes)\n",
    "    outp = kingman.simulate(pop,m)\n",
    "    parents = outp[0]\n",
    "    times = np.array(outp[1])*2*pop\n",
    "    times[0] = 0\n",
    "    finalStage = life+np.log(pop/progen)+np.log(1+times*body/pop-life-np.log(pop/progen))\n",
    "    times = np.where(times < life, times, np.where(times < life+pop/progen-1, life+np.log(1+times-life), finalStage))\n",
    "    nodes = len(times)\n",
    "\n",
    "    #Determine coalescent tree structure and identify mutations\n",
    "    branch = [0]*nodes\n",
    "    branchearly = [0]*nodes\n",
    "    for n in range(1,nodes):\n",
    "        branch[n] = np.maximum(0,np.minimum(times[parents[n]],life)-times[n])\n",
    "        branchearly[n] = np.maximum(0,times[parents[n]]-np.maximum(times[n],life))\n",
    "    branch[nodes-1] = 0\n",
    "    branchearly[nodes-1] = max(0,np.log(pop/progen*body)+life-times[nodes-1])\n",
    "    np.random.seed(m+M)\n",
    "    muts = np.random.poisson(np.array(branch)*U+np.array(branchearly)*u)\n",
    "\n",
    "    #Trace mutations down ancestry to calculate variant frequency\n",
    "    cells = [0]*nodes\n",
    "    for n in [i for i in range(len(times)) if times[i] == 0]:\n",
    "        while n > -1:\n",
    "            cells[n] = cells[n]+1\n",
    "            n = parents[n]\n",
    "    freq = np.array(cells)/pop\n",
    "    \n",
    "    #Count mutations that exceed a VAF threshold of 0, thresh1, or thresh2\n",
    "    y = np.concatenate([y, freq])\n",
    "    weighty = np.concatenate([weighty, muts])    \n",
    "    ny0[m] = np.dot(freq>0, muts)\n",
    "    ny1[m] = np.dot(freq>thresh1, muts)\n",
    "    ny2[m] = np.dot(freq>thresh2, muts)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1fd3457af98>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4QAAAJTCAYAAABQAtgYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XmcXGWZ9/9PV3X1lu50ku5OAsqwKZfgI4phBFkjIgwM\nCurDiMgywMMmjDCOywgo4LC4Kzg/YCYqiziPM4IwGIbFEcImy8PiGASuGAJCJJDe0umlequq3x+n\nEopUd/Xp7jrdVV3f9+vVL6rOuc+p69wp+urrnHPfpyqTySAiIiIiIiKVJzbbAYiIiIiIiMjsUEEo\nIiIiIiJSoVQQioiIiIiIVCgVhCIiIiIiIhVKBaGIiIiIiEiFUkEoIiIiIiJSoapnO4Cotbf3FuW5\nGgsXNtDdPVCMXc0Z6pN86pN86pN86pN8xeqTtramqiKEUzGKkSP1fZ4d6veZpz6fHer34iiUH3WF\nMKTq6vhsh1By1Cf51Cf51Cf51Cf51CflS/92s0P9PvPU57ND/R49FYQiIiIiIiIVSgWhiIiIiIhI\nhVJBKCIiIiIiUqFUEIqIiIiIiFQoFYQiIiIiIiIVSgWhiIiIiIhIhZrzzyEUESlX6XSajo6Oou6z\ntbWVWEznAkVEpHwpPxaXCkIRkRLV0dHB9Xf8P+Y1LSjK/vp7N3HKx/6SxYsXF2z39NNP8p//eSuX\nXnrl1mXXXvtDdtxxJ4488qOhPmvDhte4+OIL+Nd/vYHf/e5pGhubeMc73snHPnY4d9xxz7SOQ0RE\nKpvyY3GpIBQRKWHzmhbQ1Nwy22FMy5133sGHP3wY73jHO2c7FBERmSOUH4tHBWEIo6OjPPvss3R2\n9o+5frfddqO6Wl0pInPfddf9M//zP8+QTqf51Kc+wyGHHMozzzzF9devIJ1Ok0wmufjiy0gkEgC8\n8MLzPP74o6xZ8wI77bQLw8PDXHLJhbzxxus0Nzdz2WXf4u/+7gyuvPIKFixYyqOPPsIjjzzEF77w\nj7N8pBLGRPlxLMqZIjIXRZUfv/jFC9lll10jzY/6jRzCmjVruPrBm6mfPz9vXf+mHs7nRPbYY49Z\niExEJBpPPfUk5557xtb3r732Z0488RQ2bPgz1177Y4aGhjjzzFP4y7/ch5deWsfXvvZPtLa2cdNN\nP+H++/+bww47AoB3vWt39tnng3z4w4exdOlSkskBzjzzHLbbbnvOPfcM1qx5gaOOOobbbruNU045\nmzvvvIMTTzxltg5bJqlQfhyLcqaIlLuZzo93372Sz372vEjzowrCkBoXLqBxUXlflhYRCWvZsr3z\nxkgMDPTj/sLWRDg6Osrrr79GW1sbP/jBt6mvb6C9fSPvec97x93v/PnNbLfd9gC0tLQwODjIIYd8\nhDPOOIljjvkU7e0bMXtXtAcnRaX8KCKVZKbz42mnncCnP31ipPlRBaGIiIRSU1PLXnvtzZe/fCHp\ndJobbvgRb3vb2/n7vz+X//iP22lomMdll12ct11VVRWZTHrr623V19ezzz778IMffGfrmVMREZFy\nEWV+fP/79448P6ogFBEpYf29m0pmXw0NDTQ01PPZz/4fkskBDjroQzQ0zOPww4/gs589nfr6OhYu\nbKGjo/0t2+2xx//iuuv+me22e9u4+/6bv/kbPv3pT/OFL3xlWjGKiEhlqJT8+NGPfpzPfva0SPNj\nVSaTiWznpaC9vXfaB/jcc8/xk9Urx7wlpqe9ndPf97GKHA/R1tZEe3vvbIdRUtQn+dQn+cL2SSU9\nZ2nDhpf40Y+u56tf/fq09tPW1pR/ilXGNd0cWSg/jqWSc2ax6XfrzFOfz46x+r2S8uPzz/+BW275\n90jzY6RXCM3sK8DHgBrgGuAB4AYgAzwLnOPuaTM7HTgTGAUuc/eVZlYP3AwsBnqBk9293cz2Ba7K\ntr3X3S+N8hhERGZLLBab8JlIc8Gtt/47d9+9kq997fLZDkVERMpAJeXHlSv/k69//RuRfk5kZbCZ\nLQf2A/YHDgZ2AL4HXOTuBwJVwNFmthT4XLbd4cCVZlYLnA2szra9Cbgou+vrgOOBA4B9zGyvqI5B\nRESi98lPforbb7+dHXb4i9kORUREpGR88pOf4vrr/y3y/BjlddHDgdXAbcCvgJXAMoKrhAB3AYcC\nHwAecfchd+8B1gJ7EhR8d+e2NbP5QK27v+juGeCe7D5ERERERERkkqK8ZbQV2BE4CtgZuAOIZQs5\nCG4DbQbmAz052421PHfZ5m3a7hJR/CIiIiIiInNalAVhJ/CCuw8DbmaDBLeNbtEEbCIo8JomWD5R\n23EtXNhAdXV8GocBLS3zAIjH8y+oxuMxWlrm0dbWlLeuElTqcReiPsmnPsmnPsmnPhEREZl5URaE\nDwPnmdn3gO2AecBvzGy5u68CjgDuB54ALjezOqAW2J1gwplHgCOz648AHnL3zWY2bGa7AusIbkst\nOKlMd/fAtA+ks7MfgFQqnbculUrT2dlfkbNOabatfOqTfOqTfJplNF+xvicqKkVE5r5Kyo8zIbKC\nMDtT6EEEBV0MOAd4CVhhZjXA88At7p4ys6uBh7LtLnT3QTO7FrjRzB4GhgkmkgE4C/gZECeYZfTx\nqI5BRGQ2dXR08NOnfk5jc2NR9tfX08eJy46bcGa2p59+kq997SvstNPOW5ctWLCQTZu6+eIXL2DH\nHXeaVhxDQ0Pce+9dfPSjx4Rqf/HFX+Hooz/J+9+/97Q+t9SY2T7AN919ec6y44G/c/cPZt9rFm4R\nkW0oPwaKlR8jfeyEu39pjMUHj9FuBbBim2UDwLFjtH0M2LdYMYqIlLLG5kaaFjXP+OcuW7Y3l156\n5VuWnXvuGUXZd1dXJ7/61e2hE95cZGZfAk4E+nOW7QWcRjALNzmzcO8N1AEPm9mveXMW7kvM7DiC\nWbjPI5iF+5MEd9DcaWZ7ufszM3dUIiIzR/mxeCItCEVEZO7p6+vjG9/4Oj09wbxf55//RTZs+DMP\nPriKCy64GIBTT/0M3/3uD7nvvl/zwAP3k0wmWbBgAVdc8R1uuuknvPzyS1x//QqOPfbTfOMbX2dg\noI+RkRTnn/9Fdt31Hdx663+wcuXttLS00t3dPZuHG5UXgU8APwUwsxbgCuB83jxBunUWbmDIzHJn\n4f5Wts1dwFdzZ+HO7m/LLNwqCEVEZkgU+TF3X1HlRxWEIiKS56mnnnzLGc/99jtg6+ubbvoJy5Z9\ngI9//H/z6quvcMUVl/LP//yvXHPN1SSTSV5+eR3bb/82mpsX0NPTww9+cA2xWIzPf/5cnn/+D5x0\n0qm8+OJaTjnldK655mqWLfsAZ5xxCk8//QeuuOJSLr/8W/ziFz/nppt+TiwW47TTTpiNLoiUu99q\nZjsBmFkc+DHweSCZ00yzcIuIlJiZzo+5+4oqP6ogFBGRPGPdEvPb3z4MwLp1a3n66Sf5zW/uBaC3\ndzPxeJzlyz/MAw/cx7PPruajH/04sViMRCLBJZdcSH19PRs3bmR0dPQt+9yyr4ceuo/h4VF6ezfz\n5z+vZ+edd6GmpgaA3Xd/9wwc8axaBrwTuJbg1tA9zOwHwH1EOAs3TH8m7kKzcI+l0mfmLjb148xT\nn8+Obfs9nR6g7tUE9XWJoux/pDZBa2vjhP++CxY0sN9+H+T73//+W5Y/+eRjLFzYwPr1L/P73z/N\nQw/dB8DAQB9Lly7gyCOP4Omnf8vvfvc7TjjheJYsaWbBgkauvPJiGhoa6OrqoLGxhkWL5pFIxGlr\naxpzX/39XbzrXbvxtre1ALDXXu9jwYKGaX8vVRCKiMik7LjjThx22B4cdthf0d3dxa9+dTsARx11\nNN/+9hX09PTw+c9/ibVr/8iDD65ixYobGRwc3Homs6oqRiaTfsu+PvOZY1mz5k/86le38/a3/wUv\nvbSOoaFBqqsTrFnjHHbYEbN2vFFz9yeAdwNkrxr+3N3Pz44hjGwWbpj+TNyFZuEeSyXPzF1smsF5\n5qnPZ8dY/d7R0cfg0AiJwZGifMbg0AgdHX3EYg0F223aNMDQ0EhePMPDo3R3D7Dddm9n+fLD3pIf\n29t7OeSQI7bmx7PP/nseffRp7rrrnrfkx02bBmhoSDI8HOx/rH01NrbywgtrWL++nerqBL///bMc\ndNChob6XhYpGFYQiIiWsr6dvVva17S0xEMx+BnDSSafyjW/8E3fc8UsGBvo59dSg3fbbvw2AAw88\nmFgsxtvfvgP19fWcffapALS0tNLR0c673/0eRkZGueaaq7fu6+6772DTph5OPfUMFi5cyAknnMxZ\nZ53KggULqa+vL8bhlx13f12zcIuIjK1S8mPuvqLKj1WZTKYoOypV7e290z7A5557jp+sXknjopa8\ndT3t7Zz+vo+xxx57TPdjyo7OlOVTn+RTn+TTcwjzFfE5hFVFCKdiTDdHFsqPY6nknFls+t0689Tn\ns2Osfq+k/FgshfKjrhCKiJSoWCw24TORREREKo3yY3HN3TJYREREREREClJBKCIiIiIiUqFUEIqI\niIiIiFQoFYQiIiIiIiIVSgWhiIiIiIhIhVJBKCIiIiIiUqFUEIqIiIiIiFQoFYQiIiIiIiIVSgWh\niIiIiIhIhVJBKCIiIiIiUqFUEIqIiIiIiFQoFYQiIiIiIiIVSgWhiIiIiIhIhVJBKCIiIiIiUqFU\nEIqIiIiIiFQoFYQiIiIiIiIVSgWhiIiIiIhIhVJBKCIiIiIiUqFUEIqIiIiIiFQoFYQiIiIiIiIV\nSgWhiIiIiIhIhVJBKCIiIiIiUqEmXRCa2fwoAhERESl3ypEiIlJuqidqYGZHAQcC/wT8P6DNzC52\n9/8vxLZPA5uzb18CLgduADLAs8A57p42s9OBM4FR4DJ3X2lm9cDNwGKgFzjZ3dvNbF/gqmzbe939\n0skcsIiISLFMM0fuA3zT3Zeb2fuAHwIpYAg4yd3fUH4UEZGohblCeDFwPXAc8ASwE3DKRBuZWR1Q\n5e7Lsz+nAN8DLnL3A4Eq4GgzWwp8DtgfOBy40sxqgbOB1dm2NwEXZXd9HXA8cACwj5ntFfZgRURE\nimyqOfJLwI+Auuyiq4C/c/flwC+BLys/iojITAh1y6i7vwD8NXCHu/cBNSE2ey/QYGb3mtl92TOX\ny4AHsuvvAg4FPgA84u5D7t4DrAX2JEhod+e2zd6KU+vuL7p7Brgnuw8REZFZMcUc+SLwiZz3x7n7\n77Kvq4FBlB9FRGQGTHjLKPCGmf0Q2Bs4wcy+C7wSYrsB4DsEZ0DfSZC0qrKJCoLbXJqB+UBPznZj\nLc9dtnmbtrsUCmLhwgaqq+Mhwh1fS8s8AOLx/Po5Ho/R0jKPtramaX1GuarU4y5EfZJPfZJPfZKv\nTPtkSjnS3W81s51y3m8AMLP9gHOBgwiuCkaWH2H6ObJQfhxLpefMYlM/zjz1+exQv0crTEH4aeDj\nwFXu3m9m6whukZnIGmBttgBcY2adBFcIt2gCNhEksKYJlk/Udlzd3QMhQi2ss7MfgFQqnbculUrT\n2dlPe3vvtD+n3LS1NVXkcReiPsmnPsmnPslXrD6ZhT8appoj85jZp4ALgb/OjgmMND/C9HNkofw4\nlkrOmcWm3yMzT30+O9TvxVEoP4Y5pZcEOoEPmtlJBGcdP1F4EwBOBb4LYGbbE5y9vNfMlmfXHwE8\nRDDm4kAzqzOzZmB3gglnHgGOzG3r7puBYTPb1cyqCM6ePhQiFhERkShMNUe+hZmdQHBlcLm7r8su\nVn4UEZHIhblC+G/AjsDzBLODkv3vTRNs92PgBjN7ONv+VKADWGFmNdn93eLuKTO7miBxxYAL3X3Q\nzK4FbsxuP0wwUB7gLOBnQJxgFrXHwx2qiIhI0U01R25lZnHgaoJbTX9pZgAPuPvFyo8iIhK1MAXh\nnsDuOWP/QnH33CSV6+Ax2q4AVmyzbAA4doy2jwH7TiYWERGRiEwpRwK4+8u8mc8WjdNG+VFERCIV\n5pbR54GlUQciIiJShpQjRUSkrIW5QtgAuJk9SzANNgDufkhkUYmIiJQH5UgRESlrYQrCKyKPQkRE\npDwpR4qISFmb8JZRd3+A4AzoRwmm1l6QXSYiIlLRlCNFRKTcTVgQmtmXgEsIZj97CbjQzC6IOC4R\nEZGSpxwpIiLlLswtoycA+7h7EsDMVgBPodtkRERElCNFRKSshZllNLYl0WUNAqMRxSMiIlJOlCNF\nRKSshblC+BszuxW4Ifv+ZOC+yCISEREpH8qRIiJS1sIUhOcDZwEnEVxRvA/4lyiDEhERKRPKkSIi\nUtbGLQjNbKm7vw7sANyZ/dlie4IB9CIiIhVHOVJEROaKQlcIfwQcBTwAZHKWV2Xf7xJhXCIiIqVM\nOVJEROaEcQtCdz8q+3KZu3flrjOznaIMSkREpJQpR5amdDpNR0fHpLZpbW0lFgszx56IyNxU6JbR\nHQjOdP6XmR2Rfb1lm/8C3hV9eCIiIqVHObI0tbdv5Me/vYmG5sZQ7Qd6+jhtv5NYsmRpxJGJiJSu\nQreMXgp8iGAsxIM5y0eBlVEGJSIiUuKUI0tQZ2cXr/cO0VQzL1T73t4hOju7VBCKSEUrdMvoqQBm\n9mV3/+bMhSQiIlLalCNLV7w6QaKmLnRbEZFKF+axEzeY2d8DjQS3xMSBnd39pEgjExERKX3KkSIi\nUtbCjKK+FXgfcAIwD/gYkI4yKBERkTKhHCkiImUtTEHY6u4nA78CfgksB94dZVAiIiJlQjlSRETK\nWpiCsDv7Xwfe6+49gG66FxERUY4UEZEyF2YM4X1m9gvgC8C9ZvZ+YDDasERERMqCcqSIiJS1Ca8Q\nuvuFwD+6+5+A4wnOgn486sBERERKnXKkiIiUuwkLQjO71d1fBHD3p9z9+8DNkUcmIiJS4pQjRUSk\n3I17y6iZ3Qa8F3ibma3bZptXow5MRESkVClHiojIXFFoDOHJwCLgKuBzOctHgTeiDEpERKTETTtH\nmtk+wDfdfbmZvQO4AcgAzwLnuHvazE4Hzszu9zJ3X2lm9QRXIRcDvcDJ7t5uZvtm4xkF7nX3S4tw\nnCIiMseNe8uou29295eB7wI75vzsCuw3I9GJiIiUoOnmSDP7EvAjoC676HvARe5+IMED7o82s6UE\nxeb+wOHAlWZWC5wNrM62vQm4KLuP6wjGMR4A7GNmexXhUEVEZI4LM8to7hnGBLAn8BDwYCQRiYiI\nlI+p5sgXgU8AP82+XwY8kH19F3AYkAIecfchYMjM1mb3fwDwrZy2XzWz+UDtlvGMZnYPcCjwzNQP\nTUREKsGEBaG7fyj3vZntDHw/sohERETKxFRzpLvfamY75SyqcvdM9nUv0AzMB3py2oy1PHfZ5m3a\n7hL6QGZIJp2hq6uTjRs3ht6mtbWVWCzMY5NFRGQqwlwhfAt3f8nM3hVFMCIiIuVsGjkynfO6CdhE\nUOA1TbB8orYFLVzYQHV1fArhBlpa5gEQj4cr2IYHBvj1n3/DH6peCNW+r6eX8w47gyVLtgsdTzwe\nCx1PPB6jpWUebW1NEzcuQeUadzlTn88O9Xu0JiwIzex6gkHuEIxr2J1gwLuIiEhFK2KOfMbMlrv7\nKuAI4H7gCeByM6sDanP2/QhwZHb9EcBD7r7ZzIbNbFdgHcGYwwknlenuHphCqG/q7OwHIJVKT9Ay\nMDI8zOs91Yw2Dodq39vZz5o1r1Bd3Rg6nlQqHTqeVCpNZ2c/7e29odqXkra2prKMu5ypz2eH+r04\nChXVYa4Qrsp5nQF+Afz39EISERGZE1blvJ5OjvwHYIWZ1QDPA7e4e8rMriYYkxgDLnT3QTO7FrjR\nzB4GhgkmkgE4C/gZECeYZfTxqRxQ1OKJBImauokbAvHqRMTRiIhImDGEN5pZE7AwZ/FS4JWJtjWz\nxcBTwEcIpsG+AU2pLSIic8R0cmR2ltJ9s6/XAAeP0WYFsGKbZQPAsWO0fWzL/kRERMKa8CZ7M/s2\nsJ7gLOgqglnQVoXYLgH8C5DMLtKU2iIiMqdMNUeKiIiUijC3jB4DvM3d+ya57+8QFHBfyb7XlNoi\nIjLXTDVHikxbOp2mo6NjUtto1lYR2VaYgvD3BIPZQyc7M/tboN3d7zGzLQXhrEypPd0Z1KDwLGrl\nPkPZdFXqcReiPsmnPsmnPslXpn0y6RwpUiwdHR389Kmf09gcbtKdvp4+Tlx2HIsXL444MhEpJ2EK\nwp8Ca81sNcHYPQDc/ZAC25wKZMzsUOB9BLd95v72mbEptac7gxoUnkWtnGcomy7N+pRPfZJPfZJP\nfZKvWH0yC0XlVHKkSNE0NjfStKh5tsMQkTIWpiD8PnAe8KewO3X3g7a8NrNVBDOffXs2ptQWERGJ\n0KRzpEixZDJpkoODVCdrQrVPDg6SyYR7JIeIVI4wBWGPu99UhM+qmCm1RUSkYhQrR4pMWmdnF+vW\nd9KUzEzcGOjt7KJzaRdLliyNODIRKSdhCsKHzexWgsldtj5JNmwCdPflOW81pbaIiMwl08qRItMV\nr9ZzHUVkesIUhPMIxu/tn7MsQzAuUEREpJIpR4qISFkL82D6U2YiEBERkXKjHCkiIuVOD6IRERER\nERGpUCoIRUREREREKlSYMYQiIiIiFS+dTtPR0TGpbVpbW4nFdP5dRErXuAWhmb1EMDB+W1VAxt13\niSwqERGREqYcWZk6Ojr46VM/p7G5MVT7vp4+Tlx2HIsXL444MhGRqSt0hXD5TAUhIiJSZpbPdgAy\n8zKZNPHaaqrrwz0IPjYQo71946Q+Q1cURWSmjVsQuvufAMysFjgSaCQ48xkHdga+NhMBioiIlBrl\nyMo02QfBt7+8nk2dt7P9DtuHaq8riiIyG8KMIfwl0AC8A3gIOAh4NMqgREREyoRyZIWZzIPgY/Fq\nqmvjoa8oxgeryWTS0wlPRGTSwhSEBrwTuAr4CfAF4JYogxIRESkTypEyrtToCOs39jGY6AnVvrez\ni86lXSxZsjTiyERE3hTmJvU33D0DvADs6e6vAbXRhiUiIlIWlCOloHgiuKIY5idenZjtcEWkAoW5\nQvgHM/shcC3wMzPbHtBvLBEREeVIEREpc2EKws8CH3T358zsYuDDwPHRhiUiIlIWlCNlTtIzF0Uq\nR5iC8Al3fz+Au98B3BFtSCIiImVDOVLmJD1zUaRyhCkI3zCzAwmS3lDUAYmIiJQR5UiZsxqbG2la\n1DzbYYhIxMIUhHsDDwAZM4PgOUsZd49HGZiIiEgZUI4UEZGyNmFB6O5t2y7LPohXRESkoilHikyN\nxiiKlI4JC0Ize9TdP5jzPgY8CbwnysBERERKXTFzpJklgBuBnYAUcDowCtwAZIBngXPcPW1mpwNn\nZtdf5u4rzaweuBlYDPQCJ7t7+9SPTipZJpMmOThIdbImVPvk4CCZTDr0/jVGUaR0jFsQmtl9wPLs\n69z/w0fRoHkREalgEeXII4Fqd9/PzD4CXE7wCIuL3H2VmV0HHG1mjwKfI7hdtQ542Mx+DZwNrHb3\nS8zsOOAi4LwpxiKzIJPO0NXVycaNGydsm04P0NXVCZlMJLF0dnaxbn0nTclw+9/c3skfa/5IVVW4\nK3gdHR00zp+nMYoiJWDcgtDdDwEws6vcXQlFREQkK6IcuQaozl5lnA+MAPsSjFEEuAs4jODq4SPZ\nSWyGzGwtsCdwAPCtnLZfLVJcMkMG+/q586V7eXJ49YRt62sTrH1+HaM19ZHFE69OkKipC9V2ZGiY\nX//5Af4Q+2Oo9m+88jrzFy+gaToBikhRhJlU5mkzO2nbhe5+UwTxiIiIlJNi5sg+gttFXwBagaOA\ng9x9yyWaXqCZoFjsydlurOVblhW0cGED1dVTn/+mpWUeAPF4uKtCsViMeKwqdPt4PEZLyzza2sKV\nDS0t84jHY+Hjqaoik0mSTg+Eap/JJIlNIv7JHm8mnWLTcJzEcNXEjYdH6U4OU19TF0l/TrYv47EY\nC9vms3j71lDtRwaSJGoT1NclwrWvTdDa2hj6uxCV2f78SqV+j1aYgnB5zusEcCDwIKCCUEREKt3y\nnNfTzZF/D9zj7l8xsx2A+4DcAVxNwCZgc/Z1oeVblhXU3R2uEBpPZ2c/AKlUuLFj6XSaVDoTun0q\nlaazs5/29t7Q8aRS6dD77+/p5eanb2fJG4+Hav/nF19luKY+suNNp9PEY7XE4hOP24vHY1TFqiPr\nz8n2ZSqdZmgoRXJwJFT7oeFRRmNVodsPDo3Q0dFHLNYQqn0U2tqaQn8XpXjU78VRqKgOM8voKbnv\nzWwR8O/TD0tERKS8FTlHdhPcJgrQRVBgPmNmy919FXAEcD/wBHC5mdUBtcDuBBPOPEIwDvGJbNuH\nphhHxUiNjtCVrCYWcpxcR/8Q9dXlO4nsZMYoRjk+UURKS5grhNvackuLiIiIvNV0cuT3gZ+Y2UME\nVwYvIJixdIWZ1QDPA7e4e8rMriYo+GLAhe4+aGbXAjea2cPAMHD89A6lMsQT4cfJVVeHu72xVA32\n94ce5/fKiy8zOqU/E6Ohx1SIRCfMYyfuJ5juGoIH7u4C3BllUCIiIuWgmDnS3fuAvxlj1cFjtF0B\nrNhm2QBw7FQ+WypHQ8iZPeub5rG5dyj8jjMwPDRIMpkM1TyZTJKoSoVu37mxk9vbf0XLprxHf45J\nj6kQCS/MqZ9Lcl5ngA53fy6acERERMrKJTmvlSOltE2iaBseHp7UrlOjI6zf2MdgomfixsBrr26k\npqGezelwt+C2b+hm112W6DEVIhEIM4bwATPbn+Ahu9cD+0QelYiISBlQjpRyMpmi7bXXu6lvnj+p\n/U/29tvJtI9Vl87tqyJzzYQ3VpvZecBlwOeBRuBfzOwLUQcmIiJS6pQjpdxsKcIm+in38ZIiEl6Y\n0y1/S3DG83F37zSzvySYwew7UQYmIiJSBv4W5cjIZFJp1q1bG7r9unVryWhmTBGRSQlTEKbcfdjM\ntrwfBFITbWRmcYIB70YwruKs7LY3ZN8/C5zj7mkzOx04ExgFLnP3lWZWD9wMLCZ4wO7J7t5uZvsC\nV2Xb3uvul4Y9WBERkSKbUo6UcPp7eril434Wbg43LHPjS3+iefslEUclIjK3hJmL9wEz+w4wz8yO\nAe4AfhNXCJSsAAAgAElEQVRiu48CuPv+wEXA5cD3gIvc/UCC2diONrOlwOeA/YHDgSvNrBY4G1id\nbXtTdh8A1xFMpX0AsI+Z7RXqSEVERIpvqjlSQmpY0ExzW1uon4ZJjnkTEZFwBeEXgT8C/wOcBPwX\nMOH4CHe/HTgj+3ZHYBOwDHggu+wu4FDgA8Aj7j7k7j3AWmBPgoLv7ty2ZjYfqHX3F909A9yT3YeI\niMhsmFKOFBERKRVhZhlNm9n/JUhyVdnF2wOvhNh21MxuBD4O/G/gI9lCDoLbQJuB+UDudFdjLc9d\ntnmbtrtMFIeIiEgUppMjRURESkGYB9NfAPwj0Ekw9q8q+99QhZi7n2xmXwYeB+pzVjURXDXcnH1d\naPlEbce1cGED1dXxMKGOq6VlHgDxeP4F1Xg8RkvLPNramvLWVYJKPe5C1Cf51Cf51Cf5yrFPppsj\nRUREZluYSWVOA3Z19/bJ7NjMTgTe7u5XAgNAGnjSzJa7+yrgCOB+gtnYLjezOqAW2J1gwplHgCOz\n648AHnL3zWY2bGa7AusIxhwWnFSmu3tgMmGPqbOzH4BUKp23LpVK09nZT3t777Q/p9y0tTVV5HEX\noj7Jpz7Jpz7JV6w+mYWicko5UkQmKQPDQ4Mkk8lQzQcGBmhv3zipj2htbSUWCzOaSmRuCVMQvgJ0\nTWHfvwSuN7MHgQRwPvA8sMLMarKvb3H3lJldDTxEMKbxQncfNLNrgRvN7GFgmGAiGQhmK/0ZECeY\nZfTxKcQmIiJSDFPNkSIyCanREdZv7GMw0TNxY6D95fVs6ryd7XfYPlT73u7NHLnzYbS2to7bJp0e\noKOjb+t7FZAyV4QpCP8IPGxm9xNMpw2Au3+90Ebu3g/8zRirDh6j7QqCR1TkLhsAjh2j7WPAviHi\nFhERidqUcqSITF48kSBRUxeqbay6mob582ha1Byqfd+mXm73lbRsGr8grHs1weDQSNC+p48Tlx3H\n4sWLQ+1fpJSFKQj/nP2BNwfMi4iIiHKkyJzR0NxYsICsr0uQGByZwYhEZkaYWUb14HcREZExKEeK\niEi5043PIiIiIiIiFWrcgtDM5s1kICIiIuVCOVJEROaKQlcIVwGY2TUzE4qIiEjZWAXKkSIiUv4K\njSFsNLObgb/KPiPwLdz91OjCEhERKWnKkSIiMicUKggPAz4EHAg8MDPhiIiIlAXlSBERmRPGLQjd\n/VXgJjP7H+A5wLLtn3X30RmKT0REpOQoR4pUtnQ6TUdHx6S20YPspVSFeQ5hguDBu50EYw6XmNnH\n3f3xSCMTEREpfcqRIhVoYHM/t7cXfpB9Lj3IXkpZmILwKuBTW5Kbme0L/BD4QJSBiYiIlAHlSJFS\nlIHhoUGSyWSo5snBQRrrwvxZ/KaJHmQvUi7CfPMbc890uvtjYw2gFxERqUDKkSIlKDU6wvqNfQwm\nekK1b9/Qza71NRFHJVKawhSEXWZ2tLv/J4CZHUNwa4yIiEilK2qONLOvAB8DaoBrCCasuQHIAM8C\n57h72sxOB84ERoHL3H2lmdUDNwOLgV7gZHdvn/KRiZS5eCJBoibc+ZlY9eSuDorMJWFGtp4BXGBm\nHWbWCVxAkIREREQqXdFypJktB/YD9gcOBnYAvgdc5O4HAlXA0Wa2FPhctt3hwJVmVgucDazOtr0J\nuGg6ByYiIpVhwtMh7v5HYB8zmwfE3L03+rBERERKX5Fz5OHAauA2YD7wReB03nysxV0Ej7tIAY+4\n+xAwZGZrgT2BA4Bv5bT96jRiERGRChH6+ri790cZiIiISLkqUo5sBXYEjgJ2Bu4gKDIz2fW9QDNB\nsZg7MGqs5VuWiUgYYSahyYySHBwBpjYJjUip0jdZRESkNHQCL7j7MOBmNkhw2+gWTcAmYHP2daHl\nW5YVtHBhA9XV8SkH3NIyD4B4PNyz1WKxGPFYldqXQftSimUm2mfSKV7r7GekLtxF/jfe6GH3BfXU\n1yVCtR+pTdDa2khbW9PEjSWP+i1aExaEZnaWu183E8GIiIiUkyLnyIeB88zse8B2wDzgN2a23N1X\nAUcA9wNPAJdnZzOtBXYnmHDmEeDI7PojgIcm+sDu7oFpBdzZGVwYTaXSodqn02lS6YzaF6F9PB6L\ndP+ldKwz1T4eqyUWH3+m0Xg8tnV/VbEYQ0OprVcMJzI4NEJHRx+xWEOo9vKmtrYm2ts1Ym26ChXV\nYU6bnFu8UEREROaUouVId18JPENQ0P0KOAf4B+BSM3uUYObRW9z9deBqgoLvPuBCdx8ErgXebWYP\nE0x2c2mxYhMRkbkrzC2jr5rZfcDjwNYbq93965FFJSIiUh6KmiPd/UtjLD54jHYrgBXbLBsAjp3K\n54qISOUKUxA+lvO6KqpAREREypBypIiIlLUwj524NDud9q4EYxTqNeOoiIiIcqSIiJS/MJPKHAL8\nKxAneGDu783sM+5+b9TBiYiIlDLlSBEJI51O09HRMaltWltbicXCzZIqMh1hbhm9kuBht3e5+wYz\nOxj4v4CSnYiIVDrlSBGZ0MDmfm5vX0nLptZQ7ft6+jhx2XEsXrw44shEws0yGsvOaAaAuz8XYTwi\nIiLlRDlSREJpaG6kaVFzqJ/G5sbZDlcqSJgrhOvN7CggY2YLCKbBfiXasERERMqCcqSIiJS1MFcI\nzwQ+A+wArAPeR/B8IxERkUqnHCkiImUtzCyjG4FPm9l8YMTdkxNtIyIiUgmUI0VEpNyFmWX0PcCN\nwF9k378AnOzuL0Ycm4iISElTjhSpUBkYHhokmQx3Dig5OEhjXZiRWiIzL8w38zrgQne/C8DMPg78\nBDg4ysBERETKgHKkSAVKjY6wfmMfg4meUO3bN3Sza31NxFGJTE2YMYT1WxIdgLvfBsyPLiQREZGy\noRwpUqHiiQSJmrpQP7FqXR2U0jXut9PM/iL78n/M7B+BHwOjBIPnHyq0UzNLEJwh3QmoBS4DngNu\nADLAs8A57p42s9MJBuWPApe5+0ozqwduBhYDvQS337Sb2b7AVdm297r7pVM5aBERkemYTo4UEREp\nJYVOVzxAULxVAcsJirYtMsDnCmx7AtDp7iea2SLgd9mfi9x9lZldBxxtZo9m97M3UAc8bGa/Bs4G\nVrv7JWZ2HHARcB7BrTmfJJjJ7U4z28vdn5nsQYuIiEzTdHKkiEhB6XSajo6OSW3T2tpKLBbm5j+R\ntxq3IHT3naex318At2RfVxGcNV1GkEAB7gIOA1LAI+4+BAyZ2VpgT+AA4Fs5bb+ancGtdstAfTO7\nBzgUUEEoIiIzapo5UkSkoIHN/dzevpKWTa2h2vf19HHisuNYvHhxxJHJXBRmllEjeKbSwtzl7n7q\neNu4e1922yaCwvAi4Dvunsk26QWaCcZZ5I7GHWt57rLN27TdZaL4RUREojKVHCkiEkZDcyNNi5pn\nOwypAGFGuN4G/Bz4/WR2bGY7ZLe9xt3/zcy+lbO6CdhEUOA1TbB8orYFLVzYQHV1fDKh52lpmQdA\nPJ5/GT4ej9HSMo+2tqa8dZWgUo+7EPVJPvVJPvVJvjLtkynlSBERkVIRpiDc5O5fn8xOzWwJcC9w\nrrv/Jrv4GTNb7u6rgCOA+4EngMvNrI5g8pndCSaceQQ4Mrv+COAhd99sZsNmtivBGMLDgQknlenu\nHphM6GPq7OwHIJVK561LpdJ0dvbT3t477c8pN21tTRV53IWoT/KpT/KpT/IVq09moaicdI4UEREp\nJWEKwhvM7HLgNwRjAQFw9wcLbHMBwe0zXzWzr2aXnQdcbWY1wPPALe6eMrOrCWZkixE8y2nQzK4F\nbjSzh4Fh4PjsPs4CfgbECWYZfTzsgYqIiERgKjlSRESkZIQpCJcDfwnsl7MsAxwy3gbufh5BAbit\nvAf1uvsKYMU2ywaAY8do+xiwb4iYZ0wmlWbdurXjrt9tt92o1rNnRETmquVMMkeKiIiUkjCVyt7u\n/s7IIylT/T093NJxPws3P5e/blMP53Mie+yxxyxEJiIiM0A5UkREylqYh5WsNrM9I4+kjDUsaKa5\nrS3vZ94CzQwlIjLHKUeKiEhZC3OFcBeCCWE2EIznqwIy7q5HPoiISKVTjhQRkbIWpiA8JvIoRERE\nypNypIiIlLUwBWHeRDBZNxUzEBERkTKkHCkiE8vA8NAgyWQyVPPk4CCNdZqUUGZGmG/ah3JeJ4AD\ngQdRshMRESlqjjSzxcBTwEcIHmNxA8Gspc8C57h72sxOB87Mrr/M3VeaWT1wM7AY6AVOdvf2KR2R\niBRdanSE9Rv7GEz0hGrfvqGbXetrIo5KJDBhQejup+S+N7NFwL9HFpGIiEiZKGaONLME8C/AlksI\n3wMucvdVZnYdcLSZPQp8DtgbqAMeNrNfA2cDq939EjM7DriIsR//JCKzJJ5IkKipC9U2pkeWyQwK\nM8votvqAnYoch4iIyFwwnRz5HeA64LXs+2XAA9nXdwGHAh8AHnH3IXfvAdYCewIHAHdv01ZERGRC\nE55+MLP7CW5XgWD2tF2AO6MMSkREpBwUK0ea2d8C7e5+j5l9Zcv+3H3LvnuBZmA+kHvP2VjLtyyb\n0MKFDVRXxycb7lYtLfMAiMfDnV+OxWLEY1VqXwbtSymWUmq/ZX3U8cRjMWpr49TXJUK1H6lN0Nra\nSFtbU6j25WauHlepCHM9+pKc1xmgw93zn8IuIiJSeS7JeT2dHHkqkDGzQ4H3EYxBXJyzvgnYBGzO\nvi60fMuyCXV3D0wh1Dd1dvYDkEqlQ7VPp9Ok0hm1L0L7eDwW6f5L6VhLpX08Htu6Pup4Uuk0Q0Mp\nkoMjodoPJId44YWX6ejoC9UeoLW1lVhsKjcLzqy2tiba23tnO4yyV6ioHrcgNLO/yL58aax17v7K\n9EMTEREpP8XOke5+UM72q4CzgG+b2XJ3XwUcAdwPPAFcbmZ1QC2wO8GEM48AR2bXHwE8NMlDEpEy\nNrC5n9vbV9KyqTVU+76ePk5cdhyLFy+euLHMeYWuED5AcLazKmdZBtieYCa1qd9jIiIiUt5mIkf+\nA7DCzGqA54Fb3D1lZlcTFHwx4EJ3HzSza4EbzexhYBg4vgifLyJlpKG5kaZFoe4WF3mLcQtCd985\n972ZNQLfBQ4HTo84LhERkZIVZY509+U5b/Oec+juK4AV2ywbAI6dzueKiEhlCnXjsJl9GPh99u17\n3P3X0YUkIiJSPpQjRUSknBWcVMbM5hE8B+lw4HQlORERkYBypIiIzAXjXiHMnvFcnX37v5ToRERE\nAsqRIiIyVxS6QvhrYAQ4DPi9mW1ZXgVk3H2XiGMTEREpVcqRIiIyJxQqCHcusE5ERKSSKUeKiMic\nUGiW0T/NZCAiIiLlQjlSRCKVgeGhQZLJZKjmycFBGusKTg0iMi59c0RERERESkhqdIT1G/sYTPSE\nat++oZtd62sijkrmKhWEIiIiIiIlJp5IkKipC9U2Vq0/6WXqQj2HUEREREREROYeFYQiIiIiIiIV\nSgWhiIiIiIhIhdINxyIiIiIiFSSdTtPR0TGpbVpbW4nFdC1pLlJBKCIiIiJSQQY293N7+0paNrWG\nat/X08eJy45j8eLFEUcms0EFoYiIiIhIhWlobqRpUfNshyElQAVhhDLpDF1dnWzcuHHM9br0LiIi\nIiIis0kFYYQG+/q586V7eXJ4dd66gZ4+TtvvJJYsWToLkYmIiIiIiKggjFRqdISuZDWxZCZvXW/v\nEJ2dXSoIRURERERk1kRaEJrZPsA33X25mb0DuAHIAM8C57h72sxOB84ERoHL3H2lmdUDNwOLgV7g\nZHdvN7N9gauybe9190ujjL8Y4okEiZq6/OXViVmIRkRERERE5E2RDWAzsy8BPwK2VEPfAy5y9wOB\nKuBoM1sKfA7YHzgcuNLMaoGzgdXZtjcBF2X3cR1wPHAAsI+Z7RVV/CIiIiIiInNdlDOavAh8Iuf9\nMuCB7Ou7gEOBDwCPuPuQu/cAa4E9CQq+u3Pbmtl8oNbdX3T3DHBPdh8iIiIiIiIyBZHdMurut5rZ\nTjmLqrKFHAS3gTYD84GenDZjLc9dtnmbtrtMFMfChQ1UV8encghbtbTMAyAez6+fY7EY8VjVpNfF\n4zFaWubR1tY0rdhmW7nHHwX1ST71ST71ST71iYjIFGVgeGiQZDIZqnlycJDGOk0lIoGZ/Cakc143\nAZsICrymCZZP1Lag7u6BqUec1dnZD0Aqlc5bl06nSaUzk16XSqXp7Oynvb132vHNlra2prKOPwrq\nk3zqk3zqk3zF6hMVlSJSiVKjI6zf2MdgomfixkD7hm52ra+JOCopFzP5ELxnzGx59vURwEPAE8CB\nZlZnZs3A7gQTzjwCHJnb1t03A8NmtquZVRGMOXxoBuMXERERESlJWyYyDPMTq9bVQXnTTH4b/gFY\nYWY1wPPALe6eMrOrCQq7GHChuw+a2bXAjWb2MDBMMJEMwFnAz4A4wSyjj89g/CIiIpExswTwE2An\noBa4DHiOac7QPdPHISIi5SXSgtDdXwb2zb5eAxw8RpsVwIptlg0Ax47R9rEt+xMREZljTgA63f1E\nM1sE/C77c5G7rzKz6whm6H6UYIbuvQlm8n7YzH7NmzN0X2JmxxHM0H3erByJiIiUjZm8ZVRERETG\n9wvgq9nXVQRX/6Y1Q/cMxS0iImVMNxDPkkw6Q1dXJxs3bhxzfWtrK7GY6nURkUrh7n0AZtYE3EJw\nhe8705yhW0Rk2tLpNB0dHZPaRn/Llg8VhLNksK+fO1+6lyeHV+etG+jp47T9TmLJkqWzEJmIiMwW\nM9sBuA24xt3/zcy+lbN6KjN0FzTdRzMVeizTWAo9jkntS6t9KcVSSu23rC+VeKbaPh6LUVsbp74u\nEap91+AQd718N63JtlDt+zb1cvaH/pa2tiWh2k9EM0hHSwXhLEmNjtCVrCaWzOSt6+0dorOzSwWh\niEgFMbMlwL3Aue7+m+ziZ8xsubuvIph1+36CGbovN7M6gslntp2h+wnenM27oOk+mqnQY5nGUuhx\nTGo/ufbxeCzS/ZfSsZZK+3g8tnV9KcQznfapdJqhoRTJwZFQ7YeGR4nX1ZGY1xCqffXQCB0dfcRi\n4doXokc1FUeholoF4SzaMj1w3vLqcGdrRERkTrkAWAh81cy2jCU8D7h6mjN0i4iIjEsFoYiISAlw\n9/MYe1bQac3QLSIiUohGeoqIiIiIiFQoFYQiIiIiIiIVSreMioiIiIhUkgwMDw2STCZDNU8ODtJY\np7JhrtK/bAnKpNKsW7d23PW77bYb1dX6pxMRERGRyUuNjrB+Yx+DiZ6JGwPtG7rZtb4m4qhktqiq\nKEH9PT3c0nE/Czc/l79uUw/ncyJ77LHHLEQmIiIiInPBeLPdjyWmCxFzmv51S1TDgmaa28I9/FNE\nRERERGQqNKmMiIiIiIhIhdIVwjKj8YUiIiIiIlIsqhzKjMYXioiIiEgpS6fTdHR0hG7f2tpKLKYb\nF2eLCsIypPGFIiIiIlKqBjb3c3v7Slo2tU7Ytq+njxOXHcfixYtnIDIZiwrCOUS3k4qIiIhIKWho\nbqRpUfNshyEhqDqYQ3Q7qYiIiIiITIYKwjlmvNtJdfVQRERERES2pQqgQujqoYiIiIhMSQaGhwZJ\nJpOhmicHB2msU5lRLvQvVUF09VBEREREJis1OsL6jX0MJnpCtW/f0M2u9TURRyXFor/ypeDVw76u\nTZzc9dfstpvlrUunB4A6TRMsIiIiMsfFEwkSNXWh2sZ0IaGs6F9LgPGvHiY393HnS/fy5PDqvHUj\n/QPsu+ADtLSMPaXweFcWR0dHWbNmzbix6IqkiIiISGWY6JmF6fQAHR19b1mm5xYWl/7qloJSoyN0\nJauJJTN56954aQN/SN3GwiVL8tb1dnVzzLr92WWXd+StW7duLf/5p9/SuHBB3jqNZxQRERGpHBM9\ns7Du1QSDQyNb3+u5hcWnglAmNN4tAvHqGmqaase8stjX1c0tfxz7NtSNL/2J5u2XjLmdiIiIiFSW\nQs8srK9LkBgcGXOdFIcKQonMeLeh9nV1j7tNoQlu0uk0LS2LiMfH/trq9gERERERkclRQSglpdAE\nNx2vvsbbt29mydu2y9+uazMHLdlv3PGMqdQoVVUQi439lde4RREREZEimcRjKvSIitmn3peSU+jK\nYlcyPeZ4xtfWb+SPr449nhGC21QTDXWTHu9YqJAc64rlloHPqdQoXV1dVFWNfcWy2AXoRAOydfVU\nREREZspkHlOhR1TMPhWEUlbGG89YXZ0gMa9+3HGJfV3d466faLzjeIXkWFcs62sTJIdG+POLr9I5\nMDzm5021AC20rrOzgwfbH2Xegqa8df3dvRzzrr+mrW3swdcqFkVERKTYwj6mYrKPqJjoJPhY9LdO\nYSoIRSh8VbJQIbntFcv48CipVIaO/iHqm+cXtQCdaF1j2yJaaxvz123o4uaBf5/0rbZTLU63XdfS\nMo/Ozv4Jt5vOGNEorpDqqquIiEjpmWhW0m1pVtKJlV1BaGYx4BrgvcAQ8H/cfexZSEQitu3Zr3g8\nRiyVpro6UXC7qRaghdYVOhM3lVttp1Oc5q6Lx2OkUukJt5vOGNGCV0inWPRGsc8t69ramrcWyZPZ\nbqbWFSrO0+k0VVWMeTt0oXUT3Ua9//7Lxlwuk6McKSJlZ6LxhplRkjmzjCYHB2lcMP6spDJ5ZVcQ\nAscAde7+QTPbF/gucPQsxyRS0qZyq+10itPcdbkF4WSvuG4RZozoeFdIp1P0RrHPREMdrdtvt7VP\nJrPdTK0rVJy3v7KB6voaFra1TGpdoduo+zf10NIyjyVLdsxbJ5OmHCkiZWWi8Ya5f0fA5Mcc6hbT\niZVjQXgAcDeAuz9mZnvPxIf2dW/K+wMOYKBnM4mRYXpq8//Ynuvr+ns2Uz1UUxKxFFo3k5+55ZdW\nKR3/bK/L/UU+4XYNdQz09OatG+5PkmgoPA5hqH9gytvO5D7LQXp0hPVvbKIvnX+lu+OVN6huqGXB\nYH7hXnDdhk5qm+ZFEq+8xYznyPHy41gm+l2t9uHbx+OxSPdfSsdaKu3D5rOZimeutC+U/yG/IBzs\nT7Kpo4va2tpQsbS/+jo3vPAi8xeGu6I42J/kr3Y5lEWL8k9uzpY99tgj0v2XY0E4H8g9hZAys2p3\nHx2rcVtbU9V0P/Dgg/fh4IP3me5uREREojajOVL5UUSk/JXjtdDNQO6gnth4iU5ERKTCKEeKiMik\nlGNB+AhwJEB2fMTq2Q1HRESkZChHiojIpJTjLaO3AR8xs98CVcApsxyPiIhIqVCOFBGRSanKZPIH\n/4uIiIiIiMjcV463jIqIiIiIiEgRqCAUERERERGpUOU4hnDazCwGXAO8FxgC/o+7r81Z/1Hga8Ao\n8BN3XzHeNmb2DuAGIAM8C5zj7uEeyFRiitwv7wN+CKSyy09y9zdm9ICKoJh9krPN8cDfufsHZ+5I\niqfI35PFwApgIRAn+J68OKMHVAQR/L9zXbbtmuzysvudMpU+yVm3D/BNd1+efT9nfs+KiIiUmkq9\nQngMUJf9g/wfge9uWWFmCeD7wGHAwcAZZrakwDbfAy5y9wMJBvAfPWNHUXzF7JerCIqe5cAvgS/P\n1EEUWTH7BDPbCziN4LtSrorZJ98CfubuBwEXAe+asaMormL2ycXA1939AKAW+OsZO4rimkqfYGZf\nAn4E5D7NeC79nq0oZrbEzJ6c7TgqhZktM7MbzOzGLf9PSfTM7MNmtsLMfmZm753teCqJmR1iZj+a\n7TjKXaUWhAcAdwO4+2PA3jnrdgfWunu3uw8DDwMHFdhmGfBA9vVdwKGRRx+dYvbLce7+u+zramAw\n+vAjUbQ+MbMW4Arg/BmLPhrF/J7sD7zdzP4b+AywaiYOIALF7JNngEVmVkXwPLmRGTmC4ptKnwC8\nCHxim33Npd+zFSP7Hf4S8KfZjqWC1BHkmDuBsrwLpUw1AGcA3yE40SUzIHv3yF689QSiTEFF3jIK\nzAd6ct6nzKw6+/Debdf1As3jbQNUuXtmm7blqmj94u4bAMxsP+Bc3vxjr9wUq09qgR8DnweS0YYc\nuWL+/7MT0O3uh5rZ1wiuJH8tyuAjUsw++f/bu/M4Sav60P+frupm9s2ZYViibOoXNQYVFPC6IEZc\nrihi3IEEVDTGqHFN4i+i3rhGDEYvGo2AwM11AeEKRCMBRYNbBKMochABUbbp7tmH7p6Zrv79UdVj\n0V1dXd39VNVTVZ/36zUvquo5z3lOHbr7PN/nbL8C/jflHtOtdG6QPJ86IaV0aUQcPCWvbvo727Ui\n4i38Plj/PrAJuBh4W9sK1eWm1nlK6QMRcSzwduCl7StZd5uh3pcBb6JzR0TlXq16B86OiIvbWKyu\n0KsB4TbKT94nFSo3KbWOrQC2zHRORJRqpO1UmdULQES8DHg38D9TSoNNK3VzZVInlOdRPQL4NOUn\nWY+OiHNSSp3YW5jl788w8LXKZ1cAH2hOkZsuyzr5BPDUlNIvIuIvKA+1/IvmFb1p5lMnM+mmv7Nd\nK6V0DnDO5PuI+Crlv31PioiXpJS+0rbCdakadf5E4AbguZSHn7+pTUXrajXqfR3lKRDvSSltbFvB\nutzUeld2ejUgvB44EfhyRBwD3FR17JfAIyLiIcAOyj1bH6O8mEGtc34SEcellL5N+Q/wt1rzFZoi\ns3qJiFOA1wHHpZQ2te4rZC6TOkkp/Qh4DECl9+OLHRoMQra/P/8JPA+4qJL2Fy35BtnLsk42UQ6Y\nAO6hPKy2E82nTmbSTX9nO1L1Qj+zLRg0KaV0cuXciw0G524+dU659/08YBfw2ZYVtovMs94/DqwH\nPhQRl6eULmldibvDPOtdGenVgPAy4FkR8T3KCxScXln5cXlK6bMR8Vbg3yn37JyXUro7IqadU8nr\nbcDnImIfyjc5nfxHIJN6iYgi8E/AXcBXIwLgupTSWW34TguV5c9Kt8j69+dfIuLPKQ8hfGWrv0xG\nspENjBMAACAASURBVKyT1wBfjIg9lG/qXtvqL5OROddJnby66e9sx6ks9HMqsLPy0d4FgyrB/tnU\nWegnpXRK80vZXeZb5ymla4BrWlbQLrOAej+tdaXsPv6Nab++iYmJ2VNJkqSeFBEvBn4GXJRSOiYi\nPg78KKX0xcrxu1NKB7a1kF3GOm8P6709rPf269VVRiVJUgNSSpfy4NVuZ1oQSRmxztvDem8P6739\nDAglSdJc1FswSM1hnbeH9d4e1nuLGRBKkqS5uJ7yYlDUWDBIzWGdt4f13h7We4vZ/SpJkuai2xfO\nyiPrvD2s9/aw3lvMRWUkSZIkqUfZQyg1UWXPwTuAE1JKV1d9fiflPRrvnEee5wPvTSn9Zg7nTKSU\n+qrerwTuBg6vXu4/Ip4O/GNK6QmV92sr6d6dUjq7Kt0FwPGU98ybdFVK6d1z/T6SJElqH+cQSs23\nm/IeaitmTdmYZ1AeQjFvKaVtlIdkvHzKodMob2o86ZXAFcCZETH1mu9JKT2u6p/BoCRJUoexh1Bq\nvnuAqylvrHrm1IMR8dfAS4Ei5Y263wUcBHw7pXRwJc17K8lHgQOAf4uIpwI3AD8EHgc8FXgz8Ezg\nIcAQcHJK6b4ZynVepUxnV66xGHg+8PaqNKcDbwU+STkQvXaO312SJEk5Zg+h1BpvA54dEc+q/jAi\nngMcCTwReDxwIPCqmTJJKX2YcoD5vJTScOXjr6eUgvK+PYcDT04pPRK4rV5ewHXA6oiIyvuTgGtT\nSpsrZTsC2B/4LvAl4PVTzn9/RPx31b+sekAlSZLUIgaEUgtUhmi+lulDR/8YOJpyT9+NwFHAY+aY\n/Q8r17iNcuD5mog4GzgWWF6nTBPABZSHhQKcCny+KsnpwJdTSuOUA8KTImJD1fGpQ0a3z7HckiRJ\najMDQqlFUkrf5PdDRycVgXMmgyrKweEHgAkePE9woE7WIwARcSTwTcq/15dQniM421zDLwAvqwR6\nAVxTyWuAcu/iSysL4FxdKdMZs31PSZIkdQ4DQqm13gY8m/I8QCjPyTs1IpZHRD9wOfAnwBZgTUSs\nj4hFwHOq8thD7fm/T6c87/AzwM3ACZQDzhmllO4C7gLeD1xU6TUEOBEYTCntn1I6uDKX8XXAa2ss\nLiNJkqQOZUAotVDV0NGByvsrgEspD/v8OfDfwBdSSluBfwD+C/gP4EdV2VxJeVGZQ6Zk/yXgiIj4\nGeVA82fA1DS1nA+8mvLw0UmnA5+eku5fgcWUA1pJkhYkIg6OiIka8+vvrGzbNJ88z4+Ig+Z4zsSU\n9ysjYntEHDjl86dHxI1V79dGxGhEvG1Kugsi4q4p8+w/MJ/vI7WCG9NLkiSp5SpB362UF0t77ORc\n9AXu1Tvnc6fu1Vv57ELgp1P24P088JOU0qcq7/8SeBrwR5T39Z2ofH4B5RE7F8y1/FI7uO2EJEmS\n2sWtmaQ2c8ioJEmS2smtmaQ2MiCUJElS27g1k9ReBoSSJElqK7dmktrHgFCSJEl54NZMUhsYEEqS\nJKnt3JpJag+3nZAkSZKkHmUPoSRJkiT1KANCSZIkSepRBoSSJEmS1KMMCCVJkiSpRxkQSpIkSVKP\nMiCUJEmSpB5lQChJkiRJPcqAUJIkSZJ6lAGhJEmSJPUoA0JJkiRJ6lEGhJIkSZLUowwIJUmSJKlH\nGRBKkiRJUo8yIJQkSZKkHtXf7gI02+Dg9oks8lmzZimbNz+QRVZdwzqZzjqZzjqZzjqZLqs6Wb9+\nRV8GxekZWbSR/jy3h/XeetZ5e1jv2ajXPtpD2KD+/mK7i5A71sl01sl01sl01sl01knn8v9de1jv\nrWedt4f13nwGhJIkSZLUowwIJUmSJKlHGRBKkiRJUo8yIJQkSZKkHmVAKEmSJEk9yoBQkiRJknpU\n1+9DKEmdqlQqMTQ0lGme69ato1DwWaAkqXPZPmbLgFCScmpoaIjzv/ZfLFuxOpP8dm7fwukveCL7\n7rtv3XQ33vhj/t//u5T3ve9Dez/79Kc/yUEHHczznndiQ9e69957OOusv+Wzn72A//7vG1m+fAUP\nf/gjeMELns3XvvbvC/oekqTeZvuYLQPCBuzZs4ef//znDA/vrHn8kY98JP39VqWk7C1bsZoVq9a2\nuxgLctVVX+OZzzyBhz/8Ee0uijI2W/tYi22mpCzYPman4/4iR8QRwCeB24EvpJS+1exr3nrrrfzT\ndy5mycqV047t3LKVt3Aqj370o5tdDElqu8985lP89Kc/oVQq8bKXvYrjj/9jfvKTGzj//M9RKpUY\nGRnhrLP+noGBAQBuueWX/PCH3+fWW2/h4IMPZdeuXbz3ve/m/vvvY9WqVfz933+Uv/zLM/nQhz7I\n6tX78f3vX8/113+Xt7/9r9v8TdWIeu1jLbaZkrpVs9rHd7zj3Rx66GFNbR87LiAEjgbuA8aBX7Tq\nosvXrGb5Qzr7KYQkNeqGG37MG9945t7399xzN6eeejr33ns3n/705xkbG+N1rzudJz7xaO6443be\n857/xbp167nwwvP41rf+gxNOeC4Ahx/+KI4++lie+cwT2G+//RgZeYDXve4v2H//A3jjG8/k1ltv\n4fnPP4nLLruM00//c6666muceurp7framgfbR0m9pNXt4ze+cSVveMObm9o+dmJA+J/Al4ANwNuB\nd7a3OJLUfY488qhpcyQeeGAnKd2ytyHcs2cP9913D+vXr+ecc/6BJUuWMji4kcc+9ogZ8125chX7\n738AAGvXrmV0dJTjj38WZ555Gied9DIGBzcScXhzv5wkSfPU6vbx1a8+hVe84tSmto+dGBA+DrgX\n2Exnll+SOtI++yzi8Y8/ine9692USiUuuOBfOPDAP+Cv/uqNfPnLl7N06TL+/u/PmnZeX18fExOl\nva+nWrJkCUcffTTnnPOxvU9ONX8RsQG4KqV0VLvLIkm9oJnt4xOecFTT28dcBVQRcTTwkZTScRFR\nAM4FjgDGgNeklG4D7qQ8h3A38P52lVWSWmHn9i25yWvp0qUsXbqEN7zhNYyMPMDTnvYMli5dxrOf\n/Vze8IbXsmTJYtasWcvQ0OCDznv0o/+Qz3zmU+y//4Ez5v3Sl76UV7ziFbz97X+zoDL2uojoozxy\n5jftLoskNVOvtI8nnvgi3vCGVze1feybmJhoWuZzERHvBE4FdqaUjomIk4EXpJT+LCKOAf4mpfTC\nueY7OLh9wV/w5ptv5rybrqw5R2Lr4CCvfdwLenKC/Pr1Kxgc3N7uYuSKdTKddTJdo3XSS/ss3Xvv\nHfzLv5zP3/3dwp7zrV+/Yvoj1h4SEX8O/AB4W0rplNnSL7SNrNc+1tLLbWbW/NvaetZ5e9Sq915q\nH3/5y19wySVfamr7mKcewl8DJwMXVd4/BfgGQErpBxExr6Eva9Yspb+/uKCCrV27DIBicfoPSbFY\nYO3aZaxfv2JB1+hUvfq967FOprNOpmu0TjZsWNXkkrTfxRdfzCWXXMI555zjz8rCPYvyyJonRcRL\nUkpfaXeBWmnPnj3ceuutczrHbTCkzlMoFGbdM7AbXHrpl7jyyv/H+9//4aZeJzd/AVNKl0bEwVUf\nrQS2Vr0fj4j+lNKeueS7efMDCy7b5P5K4+OlacfGx0sMD+/sySdGPimbzjqZzjqZzjp5sGc/+4Wc\ncsopDA5uX3C9dHNA2ci0ipTSyZW0F/daMAjlbTDO+fZFLFvd2IMUt8GQlGcvfvHLePGLX9b06+Qm\nIKxhG1DdshfmGgxKktQNqqdVVD46CVicUjq2Mq3ibGDvtIpGhou2w8R4idtvv21O58y1B2/Z6lWs\nWr9+rkWTpJ6V54DweuBE4MuVxu6mNpdHkqR2yeW0inpTKmoZ2b6dr266jofsvKWh9Ds2b+F9a1/H\nH/7hHzZcnmKx0HB5On3ax5o1S7jllsbqctLhhx/uENkF6NSflU5nvTdXnv8iXAY8KyK+B/QB7lQs\nSepJeZ1WUW9KRS2lUonFK1c0vAjNXKdlDA/vZHy81HB55pr/fBayaNZCFevXr+D662+Y+xDZ4xwi\nO18O928P6z0b9YLqXAWEKaU7gWMqr0vA69taIElqo15aRU1z5rSKNhgc3Mjnv3chS1ctbyj9zk3b\neNqGJ7N27bqGr+EQWWl2to/ZylVAKEn6vaGhIS664Yssb/DmczY7tu7g1CNfPuvKbDfe+GPe856/\n4eCDD9n72erVa9iyZTPveMffctBBBy+oHGNjY3zzm1/nxBNPaij9WWf9DS984Yt5whPcZ72K0yra\nYHh4E/dtH2PFPssaSn/P7zbyq99expoNGxpKn6dFbpq9YmueelvVeWwfy7JqHw0IJSnHlq9azoqH\ntH7riSOPPIr3ve9DD/rsjW88M5O8N20a5oorLm+4wVNNTqvIwFwXubn99tsoFPsZ2GdxQ+n7+wcY\nWLak4R68uZRn7dpl3H77bcxlP+mJ0gSbNg2zcePGWdPeemviCz+9iuVrVjeU945NW/jTTf+TRz4y\nGko/NDTE1++8muWrG5sb1ugNu3qH7WN2DAglSXOyY8cOPvzh97N1a3kK21ve8g7uvfduvvOdb/O3\nf3sWAGec8SrOPvuTXHvt1Vx33bcYGRlh9erVfPCDH+PCC8/jzjvv4PzzP8dLXvIKPvzh9/PAAzvY\nvXuct7zlHRx22MO59NIvc+WVl7N27To2b97czq+bG06ryN7OrVu5ZOhbrNl2c0PpN97xG1Yd0Fhv\nX7PLUywWuPe2O+ZUntEdO7nqjm/y412zdyjf/evfsnj50oaD2ZFtOxrOG2DwrntZuWEN/Uv2aSh9\n4YECg4OzB7LV7FFUqzWjfazOq1ntowGhJGmaG2748YOeeD75yU/Z+/rCC8/jyCOfxIte9Cf89rd3\n8cEPvo9PfeqznHvuPzEyMsKdd97OAQccyKpVq9m6dSvnnHMuhUKBt771jfzyl7/gtNPO4Ne/vo3T\nT38t5577Txx55JM488zTufHGX/DBD76PD3zgo3zlK1/kwgu/SKFQ4NWvzuUOCmqB+fTgzaXHDGDp\nHObg7djU/IcTjZanWCywdXB4TnmP79nNppF+CiOz19HQzjGW9C9qSt4A92zZydCuEttKjV1j8M7f\nsWX4cg546AENpbdHUc3S6vaxOq9mtY8GhJKkaWoNifne9/4TKN9033jjj7nmmm8CsH37NorFIscd\n90yuu+5afv7zmzjxxBdRKBQYGBjgve99N0uWLGHjxo3s2fPgdU8m8/rud69l1649bN++jbvv/h2H\nHHIo++xT7jl41KMe04JvrDzKWw9eNygODDQ05LW/f6BpeU/mP5f0hWI//YuKTetRLJVK9PVBX9/M\nPYql0gMMDe3Y+94eyN7U6vaxOq9mtY8GhJKkOTnooIM54YRHc8IJz2Hz5k1cccXlADz/+S/kH/7h\ng2zdupW3vvWd3Hbbr/jOd77N5z73BUZHR/c+yezrKzAxUXpQXq961Uu49dbfcMUVl/MHf/Aw7rjj\ndsbGRunvH+DWWxMnnPDctn1ftVfeevDUHuN7dvO7jTsYHdg6e2Lm3qN4/1330b+on7UbZl4RdvFv\nBxgd2w3YA6namtE+VufVrPbRgFCScmzH1h2zJ2pCXlOHxEB59TOA0047gw9/+H/xta99lQce2MkZ\nZ5TTHXDAgQA89alPp1Ao8Ad/8FCWLFnCn//5GQCsXbuOoaFBHvOYx7J79x7OPfef9ub1jW98jS1b\ntnLGGWeyZs0aTjnlT3n9689g9eo1LFmyJIuvL6nDNbVHcZ8iS1ctq7tIyZLFAwyMlgNCV0ltv15p\nH6vzalb72DfXsfadZnBw+4K/4M0338x5N11ZcyPdrYODvPZxL8jFEtGt5kah01kn01kn0zVaJ720\nz1JWPyfr16/oy6A4PWOhbWS99rGWu9OtDCxbwr5/8FDTLzB9sVjgrptvaVr+efqu801f2KefdQc2\n1kM4eOdvOezQDTz0kIfNmGbJ4gFGKgHhvbf/jrGxsbo9itXsUZy/Wu1DL7WPWanXPtpDKEk5VSgU\nvHmQpHmaa4/irrFRRkZGZk40sWdvQDgyOsry1e3Z9kC2j1kzIJQkSVJPa2SOYrFYYHy8PL9r8N7N\nHNbgcFRwiKnyzYBQkiRJPW+2HsVisUChEhAW+ud2C/3Atp1cPngla7c4xFT5Y0AoSZIkzcUEsw8x\nreIQU+WZAaEkSZI0B3PeBsMhpsoxA0JJkiRpjua0aI1DTJVjBoSSJElSzixd5RBTtYb9ypIkSZLU\no+whlCRJkpppPovQLPY2Xa3hT5okSZLURC5CozwzIJQkSZKazEVolFcGhJIkSVKHcxEazZf9xJIk\nSZLUo+whlCRJkvLERWjUQv7kSJIkSTnS7EVopGoGhJIkSVLONHMRGqmaPz2SJElSD3GbClUzIJQk\nSZJ6iNtUqJoBoSRJktRj3KZCk+z3lSRJkqQeZQ+hJEmS1MncpkIL4E+CJEmS1MHcpkILYUAoSZIk\ndTi3qdB8NXUOYUSsbGb+kiR1KttISVIeZPp4ICKeDzwV+F/AfwHrI+KslNL/zvI6eTJRmmDTpmE2\nbtxY87h7tkiSoDfbSElS/mXdX3wWcCrwcuBHwF8A3wa6trEb3bGTq+74Jj/eddO0Yw9s3cGrn3wa\nGzbs14aSSZJypufaSEndwY3su1vmA4hTSrdExIeAi1NKOyKiq2esju/ZzaaRfgojE9OObd8+xvDw\nJgNCSRLQe22kpO7gRvbdLeuA8P6I+CRwFHBKRJwN3JXxNXJnpkm8xf6BNpRGkpRTPdlGSuoObmTf\nvbLux30F5XkRz0gp7QRupzw0RpKkXmcbKUnKnawDwhFgGDg2Ik4DtgMnZ3wNSZI6kW2kJCl3sh4y\n+q/AQcAvgclJdRPAhRlfR5KkTmMbKSkfJmDX2CgjIyMNJR8ZHWX5Yvcu7FZZ/5/9I+BRKaXpK6xI\nktTbbCMl5cL4nt38buMORge2NpR+8N7NHLbENbC6VdYB4S+B/YB7M85XkqROZxspKTdmWhSxlkK/\nvYPdLOv/u0uBFBE/B0YnP0wpHZ/xdSRJ6jS2kZKk3Mk6IPxgxvlJktQtbCMlSbmT6SqjKaXrKD8B\nPRF4EbC68lmmImJDRPw463wlSWqWVrSREXFkRFwQEV+IiA1Z5i1J6k6Z9hBGxDuBFwP/B+gD3h0R\nj0kpZfZUNCL6gHcCv8kqT0mSmq0VbSSwGHgLcAJwLHB5hnlLUkNKpRJDQ0MNp1+3bh2FQta74alR\nWQ8ZPQU4OqU0AhARnwNuINthMq8HLgbelmGekiQ1W9PbyJTS9RFxLPB24KVZ5StJc/HAtp1cPngl\na7esmzXtjq07OPXIl7Pvvvu2oGSqJeuAsDDZ0FWMAnsyvsazgCOAJ0XES1JKX8k4/8xMlCbYtGmY\njRs31jzu0xBJ6ilNbyMj4omUg8znAmcBb8oyf0lq1NJVy1nxkFXtLoYakHVAeE1EXApcUHn/p8C1\njZ4cEUcDH0kpHRcRBeBcysHfGPCalNJtKaWTK2kvznMwCDC6cydX330dvyj8atoxn4ZIUs9pehsJ\nrATOA3YBn82w7JKkLpV1QPgWykM6T6O8YM21wD83cmJlbsWpwM7KRycBi1NKx0bEMcDZwAsn06eU\nTsmw3M0xAf2LivTX2MizONrPxESpDYWSJLVJ09vIlNI1wDVZF1yS1L0yCQgjYr+U0n3AQ4GrKv8m\nHQDc1UA2vwZOBi6qvH8K8A2AlNIPIuKo+ZRtzZql9PcX53PqXmvXLgOgWJw+vLNQKFAs9NU8NlEa\n557hnexevH3asW1Dm+ExY6xfv2JBZWu3Ti9/M1gn01kn01kn03VrnXRzG1mvfaylXptp+nylz1NZ\n8pR+8nheytOq9MVCgUWLiixZPNBQ+kX79DOwaKCh9LsXDbBu3fK6bUC3tg95kVUP4b8AzweuAyaq\nPu+rvD90tgxSSpdGxMFVH60Etla9H4+I/pTSnOZbbN78wFyS1zQ8XH4gOz4+vUevVCoxXpqY8Vix\nsIhCcXoPYV+hyPDwTgYHpweLnWL9+hUdXf5msE6ms06ms06my6pOcnrT0LVtZL32sZZ6babp55a+\nWCw0Nf88fde8pC8WC3uP56E8rUw/XioxNjbOyOjuhtKP7drDnkJfQ+lHx3YzNLSDQmFpzeO2mdmo\n1z5mEhCmlJ5feXlkSmlT9bEpDdhcbAOqS16Ya0MnSVK72UZKkvIsqyGjD6X8pPPfIuK5ldeT+f8b\ncPg8sr2e8ua9X67Mj7gpi7JKktRKtpGSpDzLasjo+4BnUJ4L8Z2qz/cAV84zz8uAZ0XE9yg3nqcv\nqISSJLWHbaSkzjYBu8ZGGRkZmT0tMDI6yvLFWa9dqWbJasjoGQAR8a6U0kcWkM+dwDGV1yXKq7FJ\nktSxbCMldbrxPbv53cYdjA5snT0xMHjvZg6rscq+8inr0P2CiPgrYDnlJ5ZF4JCU0mkZX0eSpE5j\nGympYxUHBhjYZ3FDaQv99g52kqz/b11KeWnsY4DLgROAn2Z8ja4wUZpg06ZhNm7cWPP4unXrKBQa\nWwpYktQRbCMlSbmTdUC4LqX0lIj4GPBV4IPAf2R8ja4wunMnV999Hb8o/GrasR1bd3DqkS9n3333\nbUPJJElNYhspScqdrLugNlf+m4AjUkpbgcZ2sOxBS1cuY8VDVk37t3zV8nYXTZKUPdtISVLuZN1D\neG1EfAV4O/DNiHgCMJrxNbpDndWaRkZHmZhobKNQSVLHsI2UpClKpRJDQ0N1jj/A0NCOB33m1Kps\nZRoQppTeHRGHpZR+ExGvBJ5GebltTVFvtabtw5sY3m8TGzbs14aSSZKawTZSkqZ7YNtOLh+8krVb\n1tU8vvi3A4yO7d773qlV2cs0IIyIS1NKLwZIKd0A3BAR1wDPzPI63WKm1ZqK/Y4gkqRuYxspSbUt\nXbWcFQ9ZVfPYksUDDIzurnlM2cgkIIyIy4AjgAMj4vYp+f82i2tIktSJbCMlSXmWVQ/hnwIPAT4B\nvKnq8z3A/Rldo2e4JYUkdRXbSElSbmUSEKaUtgHbIuJs4KAphw8DvpPFdXqFW1JIUvewjZQk5VnW\nq4xWT44fAP4I+C42dnM2uSWFJKlr2EZKknIn61VGn1H9PiIOAf4xy2v0BLekkKSuYxspScqjrHsI\nHySldEdEHN7Ma3Qjt6SQpO5nGylJyoOst504H5iovO0DHgX8PMtr9Aq3pJCk7mIbKUnKo6x7CL9d\n9XoC+ArwHxlfQ5KkTvTtqte2kZK6V53pT1ONjI6yfHFTBy1qFlnPIfxCRKwA1lR9vB9wV5bX6WVu\nSSFJnck2UlKvqDf9aarBezdz2JJ9WlAqzSTrIaP/AJwJDFc+6qP8FPTQLK/Ty9ySQpI6k22kpF4y\n0/SnqQr99g62W9b/B04CDkwp7cg4X1VxSwpJ6ki2kZKk3Mk6IPwZsAiwsWuDUqnE0NDQjMcdTipJ\nbWUbKUnKnawDwouA2yLiJmDP5IcppeMzvk7vqjNJd3jjMJcPXsHaLeunHXM4qSS1nW2kJCl3sg4I\n/xF4M/CbjPNVRb1JuoP3buawQzc4nFSS8sk2UpKUO1kHhFtTShdmnKemmGmSrpNyJSnXbCMlSbmT\ndQTxnxFxKfB1YNfkhzaAkiTZRkrSQs22ZkYtrqNRX9YB4TJgG/A/qj6bAGzsWqHO/MKR0VEmJkpt\nKJQkqcI2UpIW6IFtO7l88ErWblnXUHrX0Zhd1hvTn55lfpqbevMLtw9vYni/TWzYsF8bSiZJso2U\npGwsXbXcNTMy5KSzLjPj/MJCP5s2DbNx48aa59mVLkmSJPUeA8IeMbpzJ1fffR2/KPxq2jG70iVJ\nkqTeZEDYQ5auXGb3uiRJkqS9MgkII+IOyhPjp+oDJlJKh2ZxHUmSOo1tpCQpz7LqITwuo3zUBrMt\n3+v8QklakOPaXQBJkmaSSUCYUvoNQEQsAp4HLKf85LMIHAK8J4vrqDnqLd/r/EJJWhjbSElSnmU9\nh/CrwFLg4cB3gacB38/4GpqPWfYoXL7a5XslqclsIyVJuZN1QBjAI4BPAOcBbwcuyfgamod6exQO\n3ruZw5bsU/O8esNJS6UHgMUOJ5WkxthGSpJyJ+uA8P6U0kRE3AL8UUrpwsoQGeXAjHsU9s/8Y1Bv\nOOmem0Z5yWP/xOGkktQY20hJUu5kHRD+IiI+CXwa+D8RcQAwkPE11GJLV9UeTrp7p/9rJWkObCMl\nSbmT9Vi/NwBfTindDJwF7A+8MuNrSJLUiWwjJUm5k3UP4Y9SSk8ASCl9DfhaxvmrGWZbcGZx7R8T\nt6uQpDmxjZSkqerch5aP72FkdPfet/XuTTU/mc8hjIinUm70xjLOW00y3wVndm7dweW/c7sKSWqQ\nbaQkTVHvPhSgWCwwPl7a+77evanmJ+uA8CjgOmAiIqC8z9JESqmY8XWUsfksOAMzzy+UJE1jGylJ\nNcx0HwrlgLBQFRDOdm+qucu0RlNK66d+5gpqvcnhpJL0YLaRkqQ8yjQgjIjvp5SOrXpfAH4MPDbL\n66iF6ozrHhsbo7ik9r1Mve0qHE4qqRe1oo2MiGcCLweWAh9NKf00q7wlSd0pk4AwIq4Fjqu8LlUd\n2oOT5jtavXHdQ3cNcejB0x547+VwUklqeRu5FDgTeBxwAmBAKEmqK5OAMKV0PEBEfCKl9OYs8lR+\nzDSuu98x3JI0q1a2kSmlKyJiGfAm4F3NvJYkqTtkfUd/Y0ScNvXDlNKFWV0gIo4E/pLyZPx3ppTu\nzypvtYbzCyX1qFa0keuAjwLvSSltzCpfSVL3yjogPK7q9QDwVOA7QGaNHbAYeAvloTDHApdnmLda\nwPmFknrUcVWv59xGRsTRwEdSSsdV5h+eCxwBjAGvSSndBnwcWA98KCIuTyldkmH5JUldKOtVRk+v\nfh8RDwG+lPE1ro+IY4G3Ay/NMm/NzcTExLw2tAfnF0rqPQtpIyPincCpwM7KRycBi1NKx0bEjuAv\nAwAAHmlJREFUMcDZwAtTStN6ICVJqqfZk8B2AAdnmWFEPBG4AXgucBbleRJqg/E9e+a1oX09DieV\n1EPm0kb+GjgZuKjy/inANwBSSj+IiKPmU4A1a5bS3z//bRDXrl0GlPcJa0ShUKBY6DN9B6TPU1ny\nlH7yeF7K0w3pG0lbfaxYKLBoUZEliwcaKsvuRQOsW7ec9etXNJS+F2W97cS3gInK2z7gUOCqOZzf\nyHCYlcB5wC7gsxkWX/Mw3w3tZ+JwUkndaiFtZErp0og4uOqjlUD107jxiOhPKe2ZS5k2b35gLsmn\nGR4ud1iOj5dmSVlWKpUYL02YPoP0xWKhqfnn6bvmJX2xWNh7PA/l6Zb0s6WtrneA8VKJsbFxRkZ3\nN1SW0bHdDA3toFBY2lD6blUvIM66h/C9Va8ngKGU0s2NnDiH4TDXANdkV2Q1RZ39Cx1OKqlHvbfq\n9ZzayBq2AdWte2GuwaAkSZD9HMLrIuJ/UN5k93zg6DmcnsvhMFB/SEy9bu5uPwbMeGyiNM49wzvZ\nvXj7tGP337+VR61eUrOrf9E+/QwsGqh5bGygyMTECKVS7Sfa69evz8VwUockTGedTGedTNftdbLA\nNnKq64ETgS9XHprelEERJanrzDYdqZZem6KU9ZDRN1Pu2TsQ+ArwzxHx+ZTSx2Y7N6/DYaD+kJh6\n3dzdfgyoe16xsIhCcfo8wr5CYcau/rFde9hT6Kt5bNPgFj7/uy+xdkN+h5OuX7+CwcHpQXAvs06m\ns06my6pO8hxULqSNrOEy4FkR8T3Kw09PnyW9JPWketORasnLPWUrZT1k9M8oP/H8YUppuLIAzI+A\n+TR2DofRNA4nldTB/owFtJEppTuBYyqvS8Drm1NMSeou3j/Wl3Vf6HhKaVfV+1FgfJ55XQ88D8Dh\nMJKkLpBlGylJUiay7iG8LiI+BiyLiJOAM5n/AjAOh+lWC1hwZiaNjA/vtfHgknInyzZSkqRMZB0Q\nvgN4LfBT4DTg34DPNHqyw2F6w/ie3ZnvXzjb+PDtm7fxvENOYN262scNFiW1wILaSEmSmiHrVUZL\nEfF/KTdyfZWPDwDuyvI66nxZ718I9ceH79iyncuT+xtKah/bSElSHmW9yujfAn8NDFPeY6mv8t9D\ns7yOulgThpNOckKxpHayjZQk5VHWQ0ZfDRyWUhrMOF/1iGYMJ5WknLCNlCTlTtYB4V3ApozzVI9p\nxnDSemZbkMb5hZIyYhspScqdrO+wfwX8Z0R8i/Jy2gCklN6f8XWkzNRbkMbFaCRlyDZSkpQ7WQeE\nd1f+we8nzEvZaMP8QhejkZQh20hJUu5kvcro+7LMT6rWrvmFLkYjKQu2kZKkPGrOpCypSVo9v7Ce\nWnMPS6UHGBraATicVJIkSfmXyV10RCxLKe3MIi+pU9Sae7j4twOMju12OKmkvWwjJUl5llX3xbcB\nIuLcjPKTOsLkcNLJfyvXlv+7fNXydhdNUn58G2wjJUn5lNU4u+URcTHwnIiYNp4vpXRGRteRaquz\n4AwsfNEZSVoA20hJysos93xTeQ84u6xq5wTgGcBTgesyylNqWL0FZ6D1m9q7t6GkKraRkpSR2e75\npmr1PWAnyiQgTCn9FrgwIn4K3AxEJe+fp5T2ZHENaTYzLTgDrV90pt7ehs4vlHqLbaQkZavePd9U\n7Vh4sNNkXUMDlDfeHaY8P3FDRLwopfTDjK8jzU0T9zCcyUzbVdh7KPUs20hJUu5kfRf8CeBlk41b\nRBwDfBJ4UsbXkeakXXsY1mLvodSzbCMlSbmTdUC4vPpJZ0rpB7Um0EvtkKc9DN3sXupJtpGSpNzJ\nelzapoh44eSbiDiJ8tAYSZJ6nW2kJCl3su4aORO4OCI+D/QBvwZOyfgaUrbaML9wJs4vlLqabaQk\nKXcyvdNNKf0KODoilgGFlNL2LPOXmqFT5hdu37yN5x1yAuvWTT8GBotS3tlGSpLyqCldHymlnc3I\nV2qWGecXFvtzszrpji3buTy5GI3U6WwjJUl54sYcUh156j0Et7KQJElStjINCCPi9Smlz2SZp9Ru\neVqddCZuZSHln22kJCmPsr6jfSNgYye1gb2HUu7ZRkqScifrgPC3EXEt8ENg76SrlNL7M76O1H61\nVied2MPI6O6Wr05aj72HUm7YRkqScifrO9YfVL3uyzhvKVdqzS8sFguMj5faMr+wnpl6DyW1lG2k\nJCl3st524n2V5bQPA34OLHE1NXWzqfMLi8UChfFSruYXSsoH20hJyr/ZptrU0unTb7JeVOZ44LNA\nEXgy8LOIeFVK6ZtZXkeSpE5jGylJ+Vdvqk0t3TD9JutujA8BTwG+nlK6NyKeDvxfwMZOvaXW/MKK\nPM0vrPcUrFQq0dcHfX21n3h1+tMwqQ1sIyWpA/TaVJus70oLKaX7IgKAlNLNk6+lXpK3/QtnUu8p\n2P133Uf/on7WbnAxGikjtpGSpNzJOiD8XUQ8H5iIiNXAXwB3ZXwNqSN0wv6FMPNTsB1btlNcPOBW\nFlJ2bCMlSbmT9Z3p64BPAA8FbgeuAc7M+BqS2sytLKR5sY2UJOVO1quMbgReERErgd0ppekTqKRe\n1yHzC2fTa+PrpYWyjZQk5VHWq4w+FvgC8LDK+1uAP00p/TrL60idrFPmFzaDQ03Vy2wjJUl5lHVX\nxGeAd6eUvg4QES8CzgOenvF1pI424/zCYn/H9x5ODfpKpQcYGtoBwNDQEF+/82qWr14x7TyHmqoH\n2EZKknIn67vLJZMNHUBK6bKIeE/G15C6Vjf0Hk6dX7j4twOMju0GyiuXrtx3tUNN1atsIyVJuZNJ\nQBgRD6u8/GlE/DXweWAP8Crgu1lcQ+oVnbI6aT3V8wuXLB5gYLQcEO7Ysn3GcxxOqm5lGylJyrOs\n7jCvAyaAPuA4yiupTZoA3pTRdaTe1SWL0czElUvVxWwjJUm5lckdZErpkCzykTSzbhhOOhtXLlU3\nakcbGRHHA69MKb2m1deWpFyp80C9lm54yD5XWa8yGpT3VFpT/XlK6YwsryP1qm4YTir1qla1kRHx\ncODxwPQ/FpLUY+o9UK+lWx6yz0XWd5GXAV8EfpZxvpI0I+cfqkO0pI1MKd0GnB0RFzfzOpLUKWZ6\noF5LLz5kz/obb0kpvT/jPCXNpsvnF84W8LmdhTqEbaQkKXeyvku8ICI+AFxDeQU1AFJK38n4OpKq\ndPv8wnoLzoDbWahjLLiNjIijgY+klI6LiAJwLnAEMAa8ptI7KElSw7IOCI8Dngg8ueqzCeD4jK8j\naYpun19Yb8GZZmxn4TBUNcFxLKCNjIh3AqcCOysfnQQsTikdGxHHAGcDL5xMn1I6JYMyS5K6XNZ3\nikellB6RcZ7TRMQzgZcDS4GPppR+2uxrSupM893OYmhoiItu+CLLVy2f03lSHQttI38NnAxcVHn/\nFOAbACmlH0TEUXPNcM2apfT3F+ddoLVrlwFQLDb2cKRQKFAs9Jm+A9LnqSx5Sj95PC/l6Yb0jaSt\nPtbsshcLBRYtKrJk8UBD6XcvGmDduuWsXz992kqnyDogvCki/iil1OxFZZZSXqntccAJgAGhNJMu\nn1/YiPluZ7HcbTCUrQW1kSmlSyPi4KqPVgLV48THI6I/pbSHBm3e/MB8irLX8HC5s3J8vNRQ+lKp\nxHhpwvQZpC8WC03NP0/fNS/pi8XC3uN5KE+3pJ8tbXW9t6Ls46USY2PjjIzubij96NhuhoZ2UCgs\nbSh9u9QLWLO+EzwU+ElE3AvsorwJ70RK6dAsL5JSuiIillHezPddWeYtdZtun1/YDg4n1Txl3UZu\nA6pb+MJcgkFJkiD7gPCkjPOrKSLWAR8F3pNS2tiKa0qdrNvnF85XvcBuaGgIJiZqHpvvMFT1vKzb\nyOuBE4EvV+YQ3pRx/pKkHpD13eDTZ/j8wkYzaHAFtY8D64EPRcTlKaVLFlhuST2oXmC3d+XSGc6d\n7zBU9bQFt5FTXAY8KyK+R7m38fR55iNJ6mFZB4TPqHo9ADwV+A4NNnaNrqCWUjotuyJLPcz5hTMG\ndvVWLpXmaUFtJEBK6U7gmMrrEvD6DMsnSepBmd7tpZQe9HQyIh4CfGkOWeRuBTWov4pavZWLuv0Y\nkJuyzFbOVl6zWCzk6vvXOzZRGuee4Z3sXjw9+Ln/vi0cuqQAEzWmJE2Ms3ifJTVX4Fq0Tz8DiwYe\ndGzyda1j9c5r5Fiz8p3vsUZXG+vk1ciapdvrJIM2UpKkzDX78f8O4OBGE+dxBTWov4pavZWLuv0Y\nkJuyzFbOVl1zciWsPH3/2Y4VC4soFKcvLFMqlbj9t1vYVpq+atbgnRs57NA+VtZYgWts1x72FPr2\nrs61ZPHA3tdTj9U7r9Fjzcp3vscaWW1s/foVDA7aA1ktqzrpsKByTm2kJEnNkGlAGBHforzJLpTn\nMxwKXLWALF1BTWojF6ORstOENlKSpAXL+q7uvVWvJ4ChlNLNC8jPFdSkPHLu4Ywa2ZJCPeu9Va+z\naCMlSVqwTO7aIuJhlZd31DqWUrprnlm7gpqUQ+5tOLNGtqTYsMHVSXtJE9tISZIWLKvH+NdRftrZ\nV/XZBHAA5ZXUGl7VxRXUpM7gcNKZZb0lRSO9joVC7QWVlAuZtZGSJGUtkzu3lNIh1e8jYjnlLSKe\nDbw2i2tI6hBTh5NO7Nm7+Erd4aQTE4yNjTkMtYahoSEuuuGLLF+1fNqxyV7Hfffdtw0lUyNsIyVJ\neZb5HVZEPBP4HHA18NiUkkvpST1k6nDSyZVXATbes4mHFUqsWDO992zr9u0M7hxjW2n6cNNeGIZa\nrxdwaGiI5SuXZdrrqPawjZQk5U1mAWFELAM+TuWJZ0rp6qzyltRZqoeTFosFClVbXsw09/Ce+zaz\nZNXKrh6GOhn03X//coaGdjzo2NDQEF+/82qWr56+bcL9d93Hyn1X01EbKuhBbCMlSXmV1aIy1U88\n/zCltGOWUyT1qJnmHvb31950vptMLjhz4Mj+jI49eA/DvUFfjV7AHVvm14nk3MN8sI2UJOVZVo/d\nrwZ2AycAP4uIyc/7gImU0qEZXUeSOtrSVctZuXYVA1M2tZ9v0FePcw9zwzZSkpRbWQWEh8yeRJKU\nNecedgTbSElSbmW1yuhvsshHkjQ39fY9dO5hPthGSpLyrDtWapCkHjbTvofOPZQkSbMxIJSUf1P3\nNpyi1/cpnI/ZhprOtOKpcw8lSeou3kFJyr2pextOVW9/Q4PF2hoaalqj19HeQ0mSuot3SZI6wkzb\nVUyaKWDshU3t52s+Q03rBZLbN2/jeYecwLp104+VSiX6+qCvr3awuHbtsjmUXJKkfJjtQWkteXt4\nakAoqSvMFDB2y6b2eVIvkLw8zdzr2L+on7Ubph/bsXUHf7XuNRQKS5tSXkmSmqXeg9Ja8jj1wjsl\nSVJm6gWLxcUDboEhSeo6M7V9ncKAUFJ3q7MgjfMLJUlSr/NOSFJXq7cgjfMLJUlSrzMglNT1Mp9f\nODHB2NiYvY6SJKnjedciSbXUCfq2bt/O4M4xtpWm9y7a6yhJkjqJAaEk1TAyNsbdm7fXDPruuW8z\nS1atdFVTSZLU8bxzkaQZ9PfvUzPo6+8faENpJEmSsmdAKKl31VmBdNeuXW0okCRJaqs69wa1dMPa\nAZ1deklagHorkE4OC5UkSb2j3r1BLd2wdoABoaSeNtMKpA4LlSSpN810b1BLN6wdUGh3ASRJkiRJ\n7WFAKEmSJEk9yoBQkiRJknqUAaEkSZIk9ajOnwUpSXlSZ7nqbliaWpIkdRfvTCQpQ/WWq+6Gpakl\nSVJ3MSCUpIzNtFx1NyxNLUmSuot3J5LUKlXDSUdGdz/okMNJJUlSO3j3IUktMjmcdPfiTYyPlx50\nbOM9m3hYocSKNaumnWewKEmSmsU7DElqocnhpIUpASEw49xDg0VJktQs3kVIUk7MNPcQZg4WXahG\nkiQthAGhJHUAF6qRJEnN4J2EJHWriQnGxsbcE7GHRMSTgddV3r45pbSlneWRJOWfdwOS1KVGxsa4\ne/N2tpWmDyl1qGnXOpNyQPgk4GXAP7e3OJKkvDMglKROVrWVxVS7du2iv38fh5r2lmJKaTQi7gWO\nb3dhJEn55x2BJHWwya0sai04c899m1myamXtE+sEkg4n7WgPRMQiYH/gvnYXRpKUf7b4ktThZlpw\npr9/YMZz6gWSDifNp4g4GvhISum4iCgA5wJHAGPAa1JKtwGfpTxMdIDfzyWUJGlGBoSS1KNmXLm0\n2D+/3kMXsWmaiHgncCqws/LRScDilNKxEXEMcDbwwpTSDcCftaeUkqROZOssSXqQ+fYeuohNU/0a\nOBm4qPL+KcA3AFJKP4iIo+aT6Zo1S+nvL867UGvXLgOgWCw0lL5QKFAs9Jm+A9LnqSx5Sj95PC/l\n6Yb0jaStPpansgMUCwUWLSqyZPHMo3Kq7V40wLp1y1m/fkVD6VvBgFCSNM18eg9dxKZ5UkqXRsTB\nVR+tBKoj9vGI6E8p7ZlLvps3P7Cgcg0Plzssx8dLDaUvlUqMlyZMn0H6YrHQ1Pzz9F3zkr5YLOw9\nnofydEv62dJW13veyg4wXioxNjbOyOjuhtKPju1maGgHhcLShtJnpV4AagstSWrYvBexUda2AdWt\ne2GuwaAkKQN1FmmrZWR0lImJxoLNVunYgDAijgdemVJ6TbvLIkm9ZD6L2Chz1wMnAl+uzCG8qc3l\nkaSeVO9BaS3bhzcxvN8mNmzYr8kla1xHBoQR8XDg8cD0OxJJkrrfZcCzIuJ7QB9wepvLI0k9a6YH\npTXT5vDhaUcGhJWltc+OiIvbXRZJklohpXQncEzldQl4fVsLJEnqCh0ZEEqSOkydORYjo6OUSiUK\njS3oJklSx5oYL3H77bfN6ZxHPvKR9DdxcbbcBYQNbrwrSeog9eZYbB/exPDwMBs2LG9DySRJap2d\nW7dyydC3WLPt5sbSb9nKWziVRz/60U0rU64CwkY33p1Mn1I6pfWllCTNx0xzLPI4n0KSpGZZunoV\nq9avb3cx9spVQEgTNt5d6Ka7UH/j3XqbV3b7MSA3ZZmtnK28ZrFYyNX3z8OxRjby7ZT/v1kcg+l/\nT/JYzlYdm/wsT5v0SpLUK3IVEDZj492FbroL9Tferbd5ZbcfA3JTltnK2aprTm6emqfv3+5jjW7k\n2wn/f7M6BtP/nuSxnK06NvnZ4OD2acfmyqBSkqS5yfsUfjfelSRJkqQmyXtAeD3wPAA33pUkSZKk\nbOVqyGgNbrwrSZIkSU2Su4DQjXclSZIkqTXyPmRUkiRJktQkBoSSJEmS1KMMCCVJkiSpRxkQSpIk\nSVKP6puYmGh3GSRJkiRJbWAPoSRJkiT1KANCSZIkSepRBoSSJEmS1KMMCCVJkiSpRxkQSpIkSVKP\nMiCUJEmSpB7V3+4CtFtEFIBzgSOAMeA1KaXbqo6fCLwH2AOcl1L63GzndLp51skAcB5wMLAI+PuU\n0tdaXfZmmk+9VB3bF7gBeFZK6ZaWFryJ5lsnEfE3wAuAfYBzU0qfb3XZm2UBvz9foPz7Mw68tpd+\nTipplgJXA69OKd3S7X9nJUnKC3sI4SRgcUrpWOCvgbMnD1Ru0v4ROAF4OnBmRGyod06XmE+dnAIM\np5SeCjwH+FTLS91886mXyWP/DIy0vMTNN+c6iYjjgCcD/6Py+UNbXegmm8/PyfOA/pTSk4H3Ax9o\neambq+7fzIg4CvgOcFij56gzVH7nf9zucvSKiDgyIi6IiC9MtkFqvoh4ZkR8LiL+T0Qc0e7y9JKI\nOD4i/qXd5eh0BoTwFOAbACmlHwBHVR17FHBbSmlzSmkX8J/A02Y5pxvMp06+AvxdJU0f5d6PbjOf\negH4GPAZ4J4WlrVV5lMnzwZuAi4DrgCubGmJm28+dXIr0F/pFVsJ7G5tkZtutr+Zi4AXAbfM4Rzl\nXET0Ae8EftPusvSQxcBbgKuAY9tcll6yFDiTcnt/QpvL0jMi4uHA4yn/3GsBen7IKOWbr61V78cj\noj+ltKfGse3AqlnO6QZzrpOU0g6AiFgBXAL8f60qbAvNuV4i4s+AwZTSv1eGSXab+fz+rAMOAp4P\nHAJ8LSIOTylNtKjMzTafOtlBebjoLZTr5/mtKWrL1P2bmVK6HiAiGj5H+RMRbwH+uPL2+8Am4GLg\nbW0rVJebWucppQ9ExLHA24GXtq9k3W2Gel8GvAl4V/tK1t1q1TtwdkRc3MZidQUDQtgGrKh6X6i6\n4Zh6bAWwZZZzusF86oSIeCjlXp9zU0r/2oqCtth86uVNwERE/DHwOODCiHhBSum+VhS4BeZTJ8PA\nLZUeshQRo8B6YGMLytsK86mTvwL+PaX0N5Xfo2sj4rEppdGWlLj55vM3s9v/znadlNI5wDmT7yPi\nq5TngD4pIl6SUvpK2wrXpWrU+RMpz1d/LnAW5TZIGatR7+uAjwLvSSl1S1uWO1PrXdkxIITrgROB\nL0fEMZSHsk36JfCIiHgI5Sf4T6M8HGCizjndYM51Upmr8E3gjSmla1pd4BaZc72klC6ZTBAR3wZe\n30XBIMzv92cUeHNEfBzYH1hGOUjsFvOpk0fx+2Gim4ABoNiyEjdfvTrJ8hw1SUQcDXwkpXRcowv+\npJROrpx7scHg3M2nzin3rJ8H7AI+27LCdpF51vvHKT/Y/FBEXF7d9qsx86x3ZcSAsNyj9ayI+B7l\nuW+nR8QrgeUppc9GxFuBf6c83/K8lNLdETHtnHYVvknmUyefANYAfxcRk3MJn5tS6qaFVOZcL20s\na6vMp07ujoinAT+qfP4XKaXxNpW/Gebz+/OPwHkR8V3KK6/+bUppZ7u+QBPUrZNGz2lNUTVVRLwT\nOBWY/Jncu+BPJVg/G3jhTOenlE5pfim7y3zrvPJAtlsfyjbdAur9tNaVsvv4N6b9+iYmumXajiRJ\nylpEvBj4GXBRSumYSu/+j1JKX6wcvzuldGBbC9llrPP2sN7bw3pvP1cZlSRJM0opXcqDV76tueBP\na0vV3azz9rDe28N6bz8DQkmSNBcu+NN61nl7WO/tYb23mAGhJEmai+uB5wG44E/LWOftYb23h/Xe\nYna/SpKkuXDBn9azztvDem8P673FXFRGkiRJknqUPYRSE0XEwcAdwAkppaurPr8TOC6ldOc88jwf\neG9K6TdzOGcipdRX9X4lcDdwePX2GBHxdOAfU0pPqLxfW0n37pTS2VXpLgCOp7xn3qSrUkrvnuv3\nkSRJUvs4h1Bqvt3A5yJixawpG/MMykMo5i2ltI3ykIyXTzl0GuVNjSe9ErgCODMipl7zPSmlx1X9\nMxiUJEnqMPYQSs13D3A15Y1Vz5x6MCL+GngpUKS8Yfm7gIOAb6eUDq6keW8l+ShwAPBvEfFU4Abg\nh8DjgKcCbwaeCTwEGAJOTindN0O5zquU6ezKNRYDzwfeXpXmdOCtwCcpB6LXzvG7S5IkKcfsIZRa\n423AsyPiWdUfRsRzgCOBJwKPBw4EXjVTJimlD1MOMJ+XUhqufPz1lFJQ3rfncODJKaVHArfVywu4\nDlgdEVF5fxJwbUppc6VsRwD7A98FvgS8fsr574+I/676l1UPqCRJklrEgFBqgcoQzdcyfejoHwNH\nU+7puxE4CnjMHLP/YeUat1EOPF8TEWcDxwLL65RpAriA8rBQgFOBz1clOR34ckppnHJAeFJEbKg6\nPnXI6PY5lluSJEltZkAotUhK6Zv8fujopCJwzmRQRTk4/AAwwYPnCQ7UyXoEICKOBL5J+ff6Espz\nBGeba/gF4GWVQC+Aayp5DVDuXXxpZQGcqytlOmO27ylJkqTOYUAotdbbgGdTngcI5Tl5p0bE8ojo\nBy4H/gTYAqyJiPURsQh4TlUee6g9//fplOcdfga4GTiBcsA5o5TSXcBdwPuBiyq9hgAnAoMppf1T\nSgdX5jK+DnhtjcVlJEmS1KEMCKUWqho6OlB5fwVwKeVhnz8H/hv4QkppK/APwH8B/wH8qCqbKykv\nKnPIlOy/BBwRET+jHGj+DJiappbzgVdTHj466XTg01PS/SuwmHJAK0nSgkTEwRExUWN+/Z2VbZvm\nk+f5EXHQHM+ZmPJ+ZURsj4gDp3z+9Ii4ser92ogYjYi3TUl3QUTcNWWe/Qfm832kVnBjekmSJLVc\nJei7lfJiaY+dnIu+wL1653zu1L16K59dCPx0yh68nwd+klL6VOX9XwJPA/6I8r6+E5XPL6A8YueC\nuZZfage3nZAkSVK7uDWT1GYOGZUkSVI7uTWT1EYGhJIkSWobt2aS2suAUJIkSW3l1kxS+xgQSpIk\nKQ/cmklqAwNCSZIktZ1bM0nt4bYTkiRJktSj7CGUJEmSpB5lQChJkiRJPcqAUJIkSZJ6lAGhJEmS\nJPUoA0JJkiRJ6lEGhJIkSZLUowwIJUmSJKlHGRBKkiRJUo/6/wEkdYU31sAANwAAAABJRU5ErkJg\ngg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1fd35dd9dd8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot VAF spectrum (fn various axis scalings\n",
    "\n",
    "fig, axs = plt.subplots(2,2, figsize=(15,10))\n",
    "\n",
    "# Linear-linear\n",
    "binx = np.linspace(0, .01, 50)\n",
    "axs[0,0].hist(x/2, edgecolor='black', bins=binx, weights=weightx/(M*u/mutrate), log=False, alpha=.6)\n",
    "axs[0,0].hist(y/2, edgecolor='black', bins=binx, weights=weighty/(M*u/mutrate), log=False, alpha=.6)\n",
    "sns.despine()\n",
    "axs[0,0].set_xlabel(\"Neutral VAF\")\n",
    "axs[0,0].set_ylabel(\"Number of neutral mutations\")\n",
    "axs[0,0].legend(['Healthy','Elevated'])\n",
    "\n",
    "# Log-linear\n",
    "binx = np.power(10,np.linspace(-np.log10(pop), np.log10(.5), 30))\n",
    "xhist = axs[0,1].hist(x/2, edgecolor='black', bins=binx, weights=weightx/(M*u/mutrate), log=False, alpha=.6)\n",
    "yhist = axs[0,1].hist(y/2, edgecolor='black', bins=binx, weights=weighty/(M*u/mutrate), log=False, alpha=.6)\n",
    "sns.despine()\n",
    "axs[0,1].set_xscale('log')\n",
    "axs[0,1].set_xlabel(\"Neutral VAF\")\n",
    "axs[0,1].set_ylabel(\"Number of neutral mutations\")\n",
    "axs[0,1].legend(['Healthy','Elevated'])\n",
    "\n",
    "# Linear-Log\n",
    "binx = np.linspace(0, .1, 50)\n",
    "axs[1,0].hist(x/2, edgecolor='black', bins=binx, weights=weightx/(M*u/mutrate), log=True, alpha=.6)\n",
    "axs[1,0].hist(y/2, edgecolor='black', bins=binx, weights=weighty/(M*u/mutrate), log=True, alpha=.6)\n",
    "sns.despine()\n",
    "axs[1,0].set_xlabel(\"Neutral VAF\")\n",
    "axs[1,0].set_ylabel(\"Number of neutral mutations\")\n",
    "axs[1,0].legend(['Healthy','Elevated'])\n",
    "\n",
    "# Log-Log\n",
    "binx = np.power(10,np.linspace(-np.log10(pop), np.log10(.5), 30))\n",
    "axs[1,1].hist(x/2, edgecolor='black', bins=binx, weights=weightx/(M*u/mutrate), log=True, alpha=.6)\n",
    "axs[1,1].hist(y/2, edgecolor='black', bins=binx, weights=weighty/(M*u/mutrate), log=True, alpha=.6)\n",
    "sns.despine()\n",
    "axs[1,1].set_xscale('log')\n",
    "axs[1,1].set_xlabel(\"Neutral VAF\")\n",
    "axs[1,1].set_ylabel(\"Number of neutral mutations\")\n",
    "axs[1,1].legend(['Healthy','Elevated'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdYAAAFBCAYAAADKeY6hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEL9JREFUeJzt3XuspHddx/HPXko3Da0uZqOUKAsiPwyGcrMX5FIXClgl\n1EZBkRK7lIuXBGO1bLjZNJIsagEhAVRYiggJaLMiKDddq8hNKRZogC8WKAiUsNTSVsqll/WPZ46d\nbs/uztn+ZubMOa9XsunMnGfm+Z2nk7zP75lnnmfDgQMHAgD0sXHeAwCAtURYAaAjYQWAjoQVADoS\nVgDoSFgBoCNhBYCOVlVYW2s7Wmuvn/c4AOBorZqwttbul+QhSbbMeywAcLRWTVir6qqqunje4wCA\nu2LVhBUA1oLNs1hJa+2UJC+rqtNbaxuTvCbJSUm+l+S8qrpqJa+3f/+NUz3B8datx+W6626a5ioW\ngu0wsB0GtsPtbIvBet8O27Ydv2G5x6c+Y22tXZDk9bn9s9OzkmypqtOS7Epyh92/VfX0aY/pSDZv\n3jTvIawKtsPAdhjYDrezLQa2w/JmMWP9fJKzk7x5dP+RSd6TJFX1kdbaw1f6glu3Hjf1/6Hbth0/\n1ddfFLbDwHYY2A63sy0GtsOdTT2sVXVpa2372EMnJLl+7P6trbXNVXXLpK857V0P27Ydn/37b5zq\nOhaB7TCwHQa2w+1si8F63w6H+qNiHgcv3ZBkfDQbVxJVAFjN5hHWDyY5M0laa6cm+dQcxgAAUzGT\no4IPsjfJGa21DyXZkOTcOYwBAKZiJmGtqquTnDq6fVuS585ivQAwa04QAQAdCSsAdCSsANCRsAJA\nR/M4KhiAdWDn7n0TLbdn144pj2S2zFgBoCNhBYCOhBUAOhJWAOhIWAGgI2EFgI6EFQA6ElYA6EhY\nAaAjYQWAjoQVADoSVgDoSFgBoCNhBYCOhBUAOhJWAOhIWAGgI2EFgI6EFQA6ElYA6EhYAaAjYQWA\njoQVADoSVgDoSFgBoCNhBYCOhBUAOhJWAOhIWAGgI2EFgI6EFQA62jzvAQCwuu3cvW/iZffs2jHF\nkSwGYQVYR0Ry+uwKBoCOhBUAOhJWAOhIWAGgI2EFgI6EFQA6ElYA6EhYAaAjYQWAjoQVADoSVgDo\nSFgBoCNhBYCOhBUAOhJWAOhIWAGgI2EFgI6EFQA6ElYA6EhYAaAjYQWAjoQVADoSVgDoSFgBoCNh\nBYCOhBUAOhJWAOhIWAGgI2EFgI6EFQA6ElYA6EhYAaAjYQWAjoQVADoSVgDoSFgBoCNhBYCOhBUA\nOhJWAOhIWAGgI2EFgI6EFQA6ElYA6EhYAaAjYQWAjoQVADoSVgDoSFgBoCNhBYCOhBUAOhJWAOhI\nWAGgI2EFgI6EFQA62jzvAQBwdHbu3jfRcnt27ZjySBhnxgoAHQkrAHQkrADQkbACQEcOXgJg1VgL\nB2SZsQJAR8IKAB0JKwB0JKwA0JGwAkBHwgoAHQkrAHQkrADQkbACQEfCCgAdCSsAdCSsANCRsAJA\nR8IKAB0JKwB0JKwA0JGwAkBHwgoAHQkrAHS0edIFW2s/leT00XMuq6orpjUoAFhUE81YW2vnJHlH\nkvsmuXeSva21ndMcGAAsoklnrOcnObmqrk2S1tpLk1yWZM+UxgUAC2nSz1g3LUU1Sarqm0lum86Q\nAGBxTTpj/URr7ZVJ3jC6/8wkn5jOkABgcU06Y31Wku9n2PV7SZKbk/zmlMYEAAtrohlrVX0nyQVT\nHgsALLzDhrW19vGqemhr7bYkB8Z+tCHJgaraNNXRAcCCOWxYq+qho//eaZdxa+3YaQ0KABbVpN9j\n/fBB9zcm+dhURgQAC+xIu4L3ZTjbUka7g5fckuTvpjcsAFhMR9oVvCNJWmt/WlXPm82QAGBxTfo9\n1ue31n4xyd0zHLi0Kcl9quolUxsZACygScN6aZLjktwvyQeSPDrJhw/7DAAmtnP3vomW27Nrx5RH\nwl016QkiWpIdSfYm+aMkJye517QGBQCLatKwfqOqDiT5bJIHVdXXkvi6DQAcZNJdwVe21l6d5LVJ\n3tJaO3EFzwWAdeOIM9bWWktyYZK3V9Wnk7wkyX2SfGO6QwOAxXPYsLbWLkxyeZLPJdnSWtuU5IFJ\nfjUuGwcAd3Kk3bnPSPITSU5MclGGE/H/SJJfrqr3TnlsALBwjhTWG6vqmiTXtNZOTvKXSZ5YVbdO\nf2gAsHiOFNbx3b3frKrzpzkYAFh0Rzp4afxScd+Z5kAAYC040oz1ga21L4xu32vs9tL1WO87vaEB\nwOI5UljvP5NRAMAacaSr23xpVgMBgLVg0lMaAgATEFYA6EhYAaAjYQWAjoQVADoSVgDoSFgBoCNh\nBYCOhBUAOjrSKQ0BWIGdu/dNvOyeXTumOBLmxYwVADoSVgDoSFgBoCNhBYCOhBUAOhJWAOhIWAGg\nI2EFgI6EFQA6ElYA6EhYAaAjYQWAjoQVADoSVgDoSFgBoCNhBYCOhBUAOhJWAOhIWAGgI2EFgI6E\nFQA6ElYA6EhYAaAjYQWAjjbPewAAq9XO3fsmXnbPrh1THAmLxIwVADoSVgDoSFgBoCNhBYCOhBUA\nOhJWAOhIWAGgI2EFgI6EFQA6ElYA6EhYAaAjYQWAjoQVADoSVgDoSFgBoCNhBYCOhBUAOto87wEA\nzMLO3fsmXnbPrh1THAlrnRkrAHQkrADQkbACQEfCCgAdCSsAdCSsANCRsAJAR8IKAB0JKwB0JKwA\n0JGwAkBHwgoAHTkJP7BwJj2hvpPpMw9mrADQkbACQEfCCgAdCSsAdCSsANCRsAJAR8IKAB0JKwB0\nJKwA0JGwAkBHwgoAHTlXMDBXzvvLWmPGCgAdCSsAdCSsANCRsAJAR8IKAB0JKwB05Os2QBeTfm0m\n8dUZ1jYzVgDoSFgBoCNhBYCOhBUAOnLwEnAnDkSCo2fGCgAdCSsAdCSsANCRsAJAR8IKAB05KhjW\nOEf4wmyZsQJAR2assEAmnX2aecL8mLECQEfCCgAd2RUMc2K3LqxNqyKsrbVHJHnO6O7zqupb8xwP\nABytVRHWJM/OENaTkzw1yZ/NdzgwOV9nAcatls9YN1XVd5Nck+Se8x4MAByt1TJjvam1dmyGqH59\n3oNhbTiazzDNPoG7auphba2dkuRlVXV6a21jktckOSnJ95KcV1VXJfnzDLt/j8ntn7XC/3OgD7Ao\nphrW1toFSc5J8u3RQ2cl2VJVp7XWTk1ycZInV9XlSX59mmPh8I5mpna0szuRBNayac9YP5/k7CRv\nHt1/ZJL3JElVfaS19vCjedGtW4/L5s2bugzwSee/Y+Jl33nxk2f6nJU8764+ZyW2bTt+3T9nluvy\nO832ObNc12p+zizXNcvfaRamGtaqurS1tn3soROSXD92/9bW2uaqumUlr3vddTf1GN6K7d9/47p/\nzizXtZqfM8t1+Z1m+5xZrms1P2eW67qrz5nXXrBDxX3WRwXfkGR8JBtXGlUAWM1mHdYPJjkzSUaf\nsX5qxusHgKma9ddt9iY5o7X2oSQbkpw74/UDwFRNPaxVdXWSU0e3b0vy3GmvEwDmZbWceQkA1gRh\nBYCOhBUAOhJWAOhIWAGgI2EFgI6EFQA6ElYA6GjDgQMH5j0GAFgzzFgBoCNhBYCOhBUAOhJWAOhI\nWAGgI2EFgI5mfaHzVa21tjHJa5KclOR7Sc6rqqvmO6r5aK19PMkNo7tfrKp1dVH61topSV5WVae3\n1u6X5JIkB5JcmeS3RtcWXvMO2g4PSfKuJP81+vFrq+pt8xvd9LXWjkmyJ8n2JMcm+cMkn846ez8c\nYjv8d9bZ+2FSwnpHZyXZUlWntdZOTXJxkifPeUwz11rbkmRDVZ0+77HMQ2vtgiTnJPn26KGXJ3lR\nVV3WWntdhvfE3nmNb1aW2Q4PS/Lyqrp4fqOauacnubaqzmmt3SPJFaN/6+39sNx2uCjr7/0wEbuC\n7+iRSd6TJFX1kSQPn+9w5uakJMe11t7XWts3+iNjPfl8krPH7j8syb+Mbr87yeNmPqL5WG47/Hxr\n7V9ba29orR0/p3HN0l8nefHo9oYkt2R9vh8OtR3W2/thIsJ6RyckuX7s/q2ttfU4q78pyZ8keUKS\n5yZ5y3raDlV1aZKbxx7aUFVLpyi7MckPzH5Us7fMdvj3JL9fVY9O8oUkfzCXgc1QVf1vVd04isbf\nJHlR1uH74RDbYd29HyYlrHd0Q5Lxv7o2VtUt8xrMHH0uyV9V1YGq+lySa5Pcc85jmqfxz8+OT/Kt\neQ1kzvZW1eVLt5M8ZJ6DmZXW2o8m+eckb66qt2advh+W2Q7r8v0wCWG9ow8mOTNJRrs/PzXf4czN\nzgyfL6e1dmKGmfw1cx3RfP1na+300e2fS/KBOY5lnt7bWjt5dPuxSS4/3MJrQWvth5O8L8nzq2rP\n6OF19344xHZYd++HSa2b3XsT2pvkjNbahzJ8jrCujoQd84Ykl7TW/i3DkY871+nMfcn5Sf6itXa3\nJJ/JsCtsPfqNJK9urd2c5OtJnj3n8czCC5JsTfLi1trSZ4zPS/KqdfZ+WG47/G6SV6yz98NEXN0G\nADqyKxgAOhJWAOhIWAGgI2EFgI6EFQA68nUbmJLW2vYkX0zy+Kp6/9jjVyc5vaquPorXfGOSC6vq\nSyt4zoGq2jB2/4QkX03ygKr66tjjj0nyiqp66Oj+D42We+H4+WBba5ck2ZHkf8ZW8/dV9cKxZbZn\nONHIpzN8de1uSd6f5Heq6rbW2uuTvK6qPjbp7wGLwowVpuvmDN+B7XUe1Z/NEKqjVlU3ZPjO9q8c\n9KNnZLiCyZKnJXlnkme31g5e50uq6sFj/16YO/va6GcnJXlQklOTPHE0hvNElbXKjBWm62sZZmoX\nZ5kv0LfWdiV5SpJNSd6b5PlJ7p3ksqraPlrmwtHi301yYpJ/aK09KsOZbj6a5MFJHpXhxAWPTXKP\nJN9McnZVff0Q49ozGtPSGba2JPmFJL83tsy5GU4C8OoMQd+3wt993JYMs9b9o/VdlmTp93pBhvNT\n/2SGs509raq+fxfWBXNlxgrTd36SJ7TWzhh/sLX2xAxXCPnpDOdZvVeSXzvUi1TV7gyhPrOqrh09\n/O6qahlOO/mAJI+oqvsnuepwr5Xh6iw/2Fpro/tnJdlXVdeNxnZShvNDfyDJ2zJcjGHcRa21K8b+\nLTcjP3H0s08k+UqSbyT55DLLPSLJb2cI649luPgDLCxhhSkb7Xp9Vu68S/hxSU7JMPP8eIbLFD5w\nhS//0dE6rsoQ8PNaaxcnOS3J3Q8zpgMZLtb9tNFD52Q4leWSc5O8vapuzRDWs0bni11y8K7gG5dZ\nzfiu4G0ZZquvWma5K6vqK6OLhX8mw4wbFpawwgxU1fty+y7hJZuSvHIpThki+9IM52ce/0zzmMO8\n9HeSpLX2sAwnSd+Y4dy1e3Pkz2LflOSpo2C2JP80eq1jMsx2nzI60Or9ozHtPNLveSijXbtvT/Iz\ny/z4u2O3D/7dYeEIK8zO+Rl2c544ur8vyTmttbuPrnf7t0l+KcNlyLa21ra11o7N6ICfkVuy/LER\nj8nwuezrMhyJ+/gM4T6kqvpyki8nuSjDpcCWThz+pCT7q+qeVbV99Fnvc5I8a5mDmFZiR4aZOaxp\nwgozMrZL+JjR/XcmuTTD7twrk1yR5E1VdX2SP07yH0n+McMFpZe8K8PBS/c56OXfluSk1tonMwT7\nk0kOXmY5b0zyzAy7hZecm+S1By331gwHIK3k888Txz6D/WySH09ywQqeDwvJ1W0AoCMzVgDoSFgB\noCNhBYCOhBUAOhJWAOhIWAGgI2EFgI6EFQA6+j+zIST1Vi84rAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1fd37f6dc18>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot ratio of neutral variants between disease and control cohorts (Figure 4, bottom row)\n",
    "yx_ratio = yhist[0] / xhist[0]\n",
    "fig, ax = plt.subplots(figsize=(7.5,5))\n",
    "ax.bar(range(len(yx_ratio)), yx_ratio)\n",
    "sns.despine()\n",
    "ax.set_xlabel(\"Neutral VAF Bin\")\n",
    "ax.set_ylabel(\"Ratio\")\n",
    "ax.set_yscale(\"log\")"
   ]
  }
 ],
 "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}