\n",
"This course is derived from the course Biological Circuit Design by Michael Elowitz and Justin Bois, 2020 at Caltech. The original course material has been changed by Matthias Fuegger and Thomas Nowak.\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This lecture covers:\n",
"\n",
"**Concepts**\n",
"\n",
"- Bio chemical reaction networks (BCRNs) to describe species and reactions among them\n",
"- The chemical master equation: definition and procedural formulation\n",
"- The Gillespie algorithm to sample from the solution of the master equation\n",
"\n",
"**Techniques**\n",
"\n",
"- Stochastic simulation of reaction networks\n",
"- Profiling of code for optimization\n",
"\n",
"\n",
"\n",
"Let's start by importing libraries we will need:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The line_profiler extension is already loaded. To reload it, use:\n",
" %reload_ext line_profiler\n"
]
}
],
"source": [
"# multiprocessing & progress\n",
"import joblib\n",
"import tqdm\n",
"\n",
"# numbers\n",
"import numpy as np\n",
"import scipy.stats as st\n",
"import numba\n",
"\n",
"# plotting\n",
"import matplotlib.pyplot as plt\n",
"figsize = (5,4)\n",
"\n",
"# profiler\n",
"# -- install if necessary --\n",
"#\n",
"# by conda (uncomment):\n",
"# %conda install line_profiler\n",
"#\n",
"# or by pip (uncomment):\n",
"# %pip install line_profiler\n",
"#\n",
"# -- end install --\n",
"%load_ext line_profiler"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Sampling from distributions\n",
"\n",
"### Flipping coins\n",
"Let's start with an example of someone doing $n$ coin flips.\n",
"The output of each coin-flip is either heads $H$ or tails $T$. \n",
"Let the **sample space** $\\Sigma = \\{H,T\\}^n$ be the set of all possible sequences.\n",
"Events are described by an **event space** $F$ that, here, contains all subsets of $\\Sigma$.\n",
"For example, $E = \\{T^n, H^n\\} \\in F$ is the event that either all flips resulted in heads or\n",
"all in tails.\n",
"Finally events are assigned **probabilities** by a function $P: F \\to [0,1]$.\n",
"In general such functions are restricted to [follow three axioms](https://en.wikipedia.org/wiki/Probability_axioms).\n",
"In our case, we will assume that:\n",
"* Heads are assumed to occur **with probability** $p \\in [0,1]$ in a single round.\n",
" Let $X_i$, with $1 \\leq i \\leq n$, be the outcome of the flip in round $i$.\n",
" We can now write $P(X_i = H)$ to denote the event that round $i$ resulted in\n",
" heads.\n",
" This is short for $P(E)$ with $E = \\left\\{ s \\in \\{H,T\\}^n \\mid s_i = H \\right\\} \\subseteq F$.\n",
" We have,\n",
" $$\n",
" P(X_i = H) = p \\quad\\text{ and }\\quad P(X_i = T) = 1-p\\enspace.\n",
" $$\n",
"* Flips in rounds are independent.\n",
" Two events $E_1$ and $E_2$ are **independent** if $P(E_1 \\cap E_2) = P(E_1) \\cdot P(E_2)$.\n",
" Thus, for $A,B \\in \\Sigma$,\n",
" $$\n",
" P(X_i = A \\wedge X_j = B) = P(X_i = A) \\cdot P(X_j = B)\\enspace,\n",
" $$\n",
" where we denoted events by Boolean formulas rather than sets; the sets\n",
" are defined as the intersection of the sets $X_i = A$ and $X_j = B$.\n",
"\n",
"\n",
"### Distribution of flipping coins\n",
"\n",
"It can be shown that the probability of having $k$ heads in a sample is,\n",
"let's denote it by $P(\\mbox{heads} = h; n,p)$, is\n",
"$$\n",
"P(\\mbox{heads} = h; n,p) = P(\\mbox{heads} = h) = {n \\choose k}p^{h}(1-p)^{n-h} \\enspace.\n",
"$$\n",
"\n",
"The observable heads is called a **stochastic variable** and its\n",
"probabilities its distribution.\n",
"The distribution of the stochastic variable heads is a [Binomial distribution (see here also for the proof)](https://www.randomservices.org/random/bernoulli/Binomial.html).\n",
"\n",
"But sometimes, it is too difficult to prove what the distribution is, so we can numerically compute its properties by **sampling** out of the distribution. Sampling involved using a random number generator to *simulate* the generation of the distribution from more fundamental events (like in our case single round coin tosses, with known probabilities).\n",
"Let's demonstrate this with the Binomial distribution.\n",
"We will take $n = 25$ and $p = 0.25$ and compute\n",
"\n",
"$$\n",
"P(\\mbox{heads} = h ; n, p)\n",
"$$\n",
"\n",
"the probability of getting $h$ heads in $n$ flips, each with probability $p$ of landing heads.\n",
"We will draw 10, 30, 100, and 300 samples and plot them versus the expected Binomial distribution."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAVAAAAEWCAYAAAAw6c+oAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAgLUlEQVR4nO3de3hU9b3v8feXCJh4iVaiG7kYdCO2GPGSgk9xF8I+WLxU7LZW2XgprQ+lu2g5z9ajrbXq2XDU1t2H2qKWWhEtyj6bXS0teKkH6q21kiAQsMVDCSCXI8FKvCRGiN/zx1oJkzAJMyuzMjOZz+t5eGbWb/3Wmt/qNB/XWr9Zv5+5OyIikr4+2W6AiEi+UoCKiESkABURiUgBKiISkQJURCQiBaiISEQKUJFuMLM7zOyX2W6HZIcCVHqUmf3SzHaZ2Xtm9qaZXddh/T+a2V/MrNHMVprZSdlqq8ihKEClp90FlLv70cAlwGwzOwfAzAYAvwJuAz4FVAP/ka2GihyKAlR6lLtvcPfm1sXw3ynh8j8BG9z9P939I+AOYJSZnZZsX2b2VTPbbGbvm1mdmU0Ny08xsxVm9o6Z7TGzRWZ2TMJ2W8zsJjNbZ2YfmtkvzOwEM3s63NfzZnZsWLfczNzMppvZzvDs+V87Oz4zO9fM/mBme81srZmNP1R7JX8pQKXHmdn9ZtYI/AXYBSwPV40E1rbWc/cPgb+G5R33cQRwH3CBux8FfA5Y07qa4Ez3RODTwBCCME50GTAROBX4IvA08F1gAMHfxQ0d6lcBw4HzgVvM7L8ladMgYBkwm+AM+kbgv8ys7BDtlTylAJUe5+7/AhwF/APBJXvrGemRQEOH6g1h3WQ+AU43s2J33+XuG8L9b3L337l7s7vXAz8CxnXY9ifu/ra77wBeAv7k7q+HZ8dPAmd1qH+nu3/o7rXAAmBKkvZcBSx39+Xu/om7/47gNsSFXbVX8pcCVLLC3Vvc/WVgMPDNsPgD4OgOVY8G3k+y/YfAFcAMYJeZLWu91Dez481ssZntMLP3gF8SnFkmejvhfVOS5SM71H8r4f1WgrPbjk4CLg8v3/ea2V7gPGBgV+2V/KUAlWw7jAP3QDcAo1pXhJe9p4TlB3H3Z919IjCQ4HbAz8NVdxHcWz0j7Ky6iuCyvjuGJLwfCuxMUuct4DF3Pybh3xHufvch2it5SgEqPSY8M7zSzI40syIz+wLBpfCKsMqTBJe4l5nZ4cD3gXXu/pck+zrBzC4JQ7aZ4Oy1JVx9VLi8N7wveVMGmn+bmZWY2UhgGsl/HfBL4Itm9oXw+A43s/FmNvgQ7ZU8pQCVnuQEl+vbgXeBe4FZ7v5rgPB+5WXAnHD9GODKTvbVB/hXgjPBvxHc4/yXcN2dwNkE90+XEdxn7a4XgE3A/wHudffnOlZw97eAyQSdUfUEZ6Q3hW3tqr2Sp0wDKot0zszKgTqgr7vvz3JzJMfoDFREJCIFqIhIRLqEFxGJSGegIiIRHZbtBmTSgAEDvLy8PNvNEJFepqamZo+7l3Us71UBWl5eTnV1dbabISK9jJltTVauS3gRkYgUoCIiESlARUQiUoCKiESkABURiUgBmu/qFsFT5fB4n+C1blG2WyRSMHrVz5gKTt0ieG06tDQGy41bg2WAYZpuRyRuOgPNZ2tvPRCerVoag3IRiZ0CNJ81bkuvXCQH7N27l/vvvx+A3//+91x88cVZblF0CtB8VjI0vXKRNK2o3cHV961g0r8t4+r7VrCidke395kYoHHZv79nhm5VgOazUXOgqKR9WVFJUC7STStqdzB3WS27G5pwYHdDE3OX1XY7RG+55Rb++te/cuaZZ3LTTTfxwQcf8OUvf5nTTjuNqVOn0jpCXE1NDePGjeOcc87hC1/4Art27QJgzZo1nHvuuZxxxhl86Utf4t133wVg/PjxfPe732XcuHHMmTOHYcOGsW/fPgDee+89ysvL25YzRQGaz4ZNhdHzoeQkwILX0fPVgSQZsWDlRpr3tZ+2qXlfCwtWbuzWfu+++25OOeUU1qxZww9/+ENef/115s6dyxtvvMHmzZt55ZVX2LdvH9dffz1LliyhpqaGr33ta9x6a3Bv/5prruGee+5h3bp1VFRUcOedd7bte+/evbzwwgvcfvvtjB8/nmXLlgGwePFiLrvsMvr27duttnekXvh8N2yqAlNiUd/QlFZ5VKNHj2bw4MEAnHnmmWzZsoVjjjmG9evXM3HiRABaWloYOHAgDQ0N7N27l3HjxgFw7bXXcvnll7ft64orrmh7f9111/GDH/yASy+9lAULFvDzn2d+ElQFqIgkVVZazO4kYVlWWpzRz+nfv3/b+6KiIvbv34+7M3LkSP74xz+2q9vQ0NDlvo444oi292PHjmXLli288MILtLS0cPrpp2e03aBLeBHpxLSqEfTvW9SurH/fIqZVjejWfo866ijef//9LuuMGDGC+vr6tgDdt28fGzZsoLS0lGOPPZaXXnoJgMcee6ztbDSZa665hilTpjBt2rRutbkzOgMVkaQmVAwCgnuh9Q1NlJUWM61qRFt5VMcddxxjx47l9NNPp7i4mBNOOOGgOv369WPJkiXccMMNNDQ0sH//fmbNmsXIkSNZuHAhM2bMoLGxkZNPPpkFCxZ0+llTp07le9/7HlOmTOlWmzvTq+ZEqqysdA2oLCKtlixZwq9//Wsee+yxbu3HzGrcvbJjeaxnoGY2CfgxUAQ85O53d1g/Fbg5XPwA+Ka7rw3XbQHeB1qA/ckaLyLSmeuvv56nn36a5cuXx/YZsQWomRUB84CJwHZglZktdfc3EqrVAePc/V0zuwCYD4xJWF/l7nviaqOI9F4/+clPYv+MODuRRgOb3H2zu38MLAYmJ1Zw9z+4+7vh4qvA4BjbIyKSUXEG6CDgrYTl7WFZZ74OPJ2w7MBzZlZjZtNjaJ+ISLfEeQ/UkpQl7bEysyqCAD0voXisu+80s+OB35nZX9z9xSTbTgemAwwdqmfARaTnxHkGuh0YkrA8GNjZsZKZnQE8BEx293day919Z/i6G3iS4JbAQdx9vrtXuntlWdlB0zaLiMQmzgBdBQw3s2Fm1g+4EliaWMHMhgK/Aq529zcTyo8ws6Na3wPnA+tjbKuI9ELjx48nzp82xnYJ7+77zWwm8CzBz5gedvcNZjYjXP8g8H3gOOB+M4MDP1c6AXgyLDsMeNzdn4mrrSLSibpFwQDdjduCYRJHzdHYCwli/R2ouy8HlncoezDh/XXAdUm22wyMirNtInIIMU0Z8+GHH/KVr3yF7du309LSwm233cbGjRv5zW9+Q1NTE5/73Of42c9+hpkxfvx4zjrrLGpqaqivr+fRRx/lrrvuora2liuuuILZs2ezZcsWJk2axJgxY3j99dc59dRTefTRRykpaT/U43PPPcftt99Oc3Mzp5xyCgsWLODII4+MfBygZ+FFpDMxTRnzzDPPcOKJJ7J27VrWr1/PpEmTmDlzJqtWrWL9+vU0NTXx29/+tq1+v379ePHFF5kxYwaTJ09m3rx5rF+/nkceeYR33gm6TTZu3Mj06dNZt24dRx999EEDNu/Zs4fZs2fz/PPPs3r1aiorK/nRj37UreMABaiIdCamKWMqKip4/vnnufnmm3nppZcoLS1l5cqVjBkzhoqKClasWMGGDRva6l9yySVt240cOZKBAwfSv39/Tj75ZN56K/il5JAhQxg7diwAV111FS+//HK7z3z11Vd54403GDt2LGeeeSYLFy5k69at3ToO0GAiItKZkqHBZXuy8m449dRTqampYfny5XznO9/h/PPPZ968eVRXVzNkyBDuuOMOPvroo7b6rcPd9enTp93Qd3369GmbuiPsL2nTcdndmThxIk888US32t6RzkBFJLmYpozZuXMnJSUlXHXVVdx4442sXr0agAEDBvDBBx+wZMmStPe5bdu2tqHvnnjiCc4777x2688991xeeeUVNm3aBEBjYyNvvvnmQftJl85ARSS51o6iDPfC19bWctNNN9GnTx/69u3LAw88wFNPPUVFRQXl5eV89rOfTXufn/70p1m4cCHf+MY3GD58ON/85jfbrS8rK+ORRx5hypQpNDc3AzB79mxOPfXUbh2LhrMTkby2ZcsWLr74Ytavj++n4p0NZ6dLeBGRiBSgIpLXysvLYz377IoCVEQkIgWoiEhEClARkYgUoCIiESlARUQiUoCKiESkABURiUgBKiISkQJURCQiBaiISEQKUBGRiBSgIiIRKUBFRCJSgIqIRKQAFRGJSAEqIhKRAlREJCIFqIhIRApQEZGIFKAiIhEpQEVEIoo1QM1skpltNLNNZnZLkvVTzWxd+O8PZjYq1W0ly+oWwVPl8Hif4LVuUbZbJNLjYgtQMysC5gEXAJ8BppjZZzpUqwPGufsZwL8B89PYVrKlbhG8Nh0atwIevL42XSEqBSfOM9DRwCZ33+zuHwOLgcmJFdz9D+7+brj4KjA41W0li9beCi2N7ctaGoNykQISZ4AOAt5KWN4elnXm68DT6W5rZtPNrNrMquvr67vRXElZ47b0ykV6qTgD1JKUedKKZlUEAXpzutu6+3x3r3T3yrKyskgNlTSVDE2vXKSXijNAtwNDEpYHAzs7VjKzM4CHgMnu/k4620qWjJoDRSXty4pKgnKRAhJngK4ChpvZMDPrB1wJLE2sYGZDgV8BV7v7m+lsK1k0bCqMng8lJwEWvI6eH5SLFJDD4tqxu+83s5nAs0AR8LC7bzCzGeH6B4HvA8cB95sZwP7wcjzptnG1VSIYNlWBKQXP3JPeWsxLlZWVXl1dne1miEgvY2Y17l7ZsVxPIomIRKQAFRGJSAEqIhJRbJ1IkltW1O5gwcqN1Dc0UVZazLSqEUyo6Oq5BhE5FAVoAVhRu4O5y2pp3tcCwO6GJuYuqwVQiIp0gy7hC8CClRvbwrNV874WFqzcmKUWifQOCtACUN/QlFa5iKRGAVoAykqL0yoXkdQoQAvAtKoR9O9b1K6sf98iplWNyFKLRHoHdSIVgNaOIvXCi2SWArRATKgYpMAUyTBdwouIRKQAFRGJSAEqIhKRAlREJCIFqIhIROqFl3Y06IhI6hSgeS6TgadBR0TSo0v4PNYaeLsbmnAOBN6K2h2R9qdBR0TSowDNY5kOPA06IpIeBWgey3TgadARkfQoQPNYpgNPg46IpEcBmscyHXgTKgYx66IKji8txoDjS4uZdVGFOpBEOqFe+DwWxyhLGnREJHUK0DynwBPJHl3Ci4hEpAAVEYlIASoiElGsAWpmk8xso5ltMrNbkqw/zcz+aGbNZnZjh3VbzKzWzNaYWXWc7RQRiSK2TiQzKwLmAROB7cAqM1vq7m8kVPsbcANwaSe7qXL3PXG1UUSkO+I8Ax0NbHL3ze7+MbAYmJxYwd13u/sqYF+M7RARiUWcAToIeCtheXtYlioHnjOzGjOb3lklM5tuZtVmVl1fXx+xqSIi6YszQC1Jmaex/Vh3Pxu4APiWmX0+WSV3n+/ule5eWVZWFqWdIiKRxBmg24EhCcuDgZ2pbuzuO8PX3cCTBLcERERyRpwBugoYbmbDzKwfcCWwNJUNzewIMzuq9T1wPrA+tpaKiEQQWy+8u+83s5nAs0AR8LC7bzCzGeH6B83s74Bq4GjgEzObBXwGGAA8aWatbXzc3Z+Jq60iIlHE+iy8uy8HlncoezDh/f8juLTv6D1gVJxtExHpLj2JJCISkQJURCQiBaiISEQK0EJRtwieKofH+wSvdYuy3SKRvKcBlQtB3SJ4bTq0NAbLjVuDZYBhU7PXLpE8pzPQQrD21gPh2aqlMSgXkci6PAM1s8OBGcDfA7XAL9x9f080TDKocVt65SKSkkOdgS4EKgnC8wLg32NvkWReydD0ykUkJYcK0M+4+1Xu/jPgy8A/9ECbJNNGzYGikvZlRSVBuYhEdqgAbRunU5fueWzYVBg9H0pOAix4HT1fHUgi3XSoXvhRZvZe+N6A4nDZAHf3o2NtnWTOsKkZDcwVtTsyOh+9SD7qMkDdvainGiL5Y0XtDuYuq6V5XwsAuxuamLusFkAhKgVFP2OStC1YubEtPFs172thwcqNWWqRSHYoQCVt9Q1NaZWL9FYKUElbWWlxWuUivZUe5cxRudxJM61qRLt7oAD9+xYxrWpEFlsl0vMUoDko1ztpWtuQqwEv0lMUoDmoq06aXAmpCRWDcqYtItmie6A5SJ00IvlBAZqD1Ekjkh8UoDloWtUI+vdt/wyDOmlEco/ugeYgddKI5AcFaI5SJ41I7tMlvIhIRApQEZGIFKAiIhEpQEVEIlKAiohEFGuAmtkkM9toZpvM7JYk608zsz+aWbOZ3ZjOtpJH6hbBU+XweJ/gtW5RtlskkhGx/YzJzIqAecBEYDuwysyWuvsbCdX+BtwAXBphW8kHdYvgtekH5qVv3Bosg+ZkkrwX5xnoaGCTu29294+BxcDkxAruvtvdV5EweV2q20qeWHvrgfBs1dIYlIvkuTgDdBDwVsLy9rAso9ua2XQzqzaz6vr6+kgNlRg1bkuvXCSPxBmglqTMM72tu89390p3rywrK0u5cdJDSoamVy6SR+IM0O3AkITlwcDOHthWcsmoOVBU0r6sqCQoF8lzcQboKmC4mQ0zs37AlcDSHthWcsmwqTB6PpScBFjwOnq+OpCkV4itF97d95vZTOBZoAh42N03mNmMcP2DZvZ3QDVwNPCJmc0CPuPu7yXbNq62SsyGTVVgSq9k7qnelsx9lZWVXl1dne1miEgvY2Y17l7ZsVxPIomIRKQAFRGJSAEqIhKRAlREJCIFqIhIRApQEZGIFKAiIhFpVk6J3YraHZqiWXolBajEakXtDuYuq6V5XwsAuxuamLusFkAhKnlPl/ASqwUrN7aFZ6vmfS0sWLkxSy0SyRwFqMSqvqEprXKRfKIAlViVlRanVS6STxSgEqtpVSPo37eoXVn/vkVMqxqRpRaJZI46kSRWrR1F6oWX3kgBKrGbUDFIgSm9ki7he1qqc6Rnay71bM7hrvnjJc/oDLQnpTpHerbmUs/mHO6aP17ykM5Ae1Kqc6Rnay71bM7hrvnjJQ8pQHuQdzIX+kHl2ZpLPZtzuGv+eMlDCtAe9E5L8nnrDyrP1lzq2ZzDXfPHSx5SgPagh/ZczUef9G9X9tEn/Xloz9XtK2ZrLvVszuGu+eMlDylAe9CGoguZWz+Tt/eV8Ykbb+8rY279TDYUXdi+YrbmUs/mHO6aP17ykKY17kEdRyaC4KmcWRdV6HeSIjmss2mN9TOmHqSnckR6FwVoD9NTOSK9h+6BiohEpAAVEYlIASoiElGsAWpmk8xso5ltMrNbkqw3M7svXL/OzM5OWLfFzGrNbI2Z5W7XuogUrNg6kcysCJgHTAS2A6vMbKm7v5FQ7QJgePhvDPBA+Nqqyt33xNVGEZHuiPMMdDSwyd03u/vHwGJgcoc6k4FHPfAqcIyZDYyxTSIiGRNngA4C3kpY3h6WpVrHgefMrMbMpnf2IWY23cyqzay6vr4+A80WEUlNnAFqSco6PvbUVZ2x7n42wWX+t8zs88k+xN3nu3ulu1eWlSUfrENEJA5xBuh2YEjC8mBgZ6p13L31dTfwJMEtARGRnBHnk0irgOFmNgzYAVwJ/HOHOkuBmWa2mKDzqMHdd5nZEUAfd38/fH8+8D9jbKvkgBW1O/SYq+SV2ALU3feb2UzgWaAIeNjdN5jZjHD9g8By4EJgE9AITAs3PwF40sxa2/i4uz8TV1sl+zoOtLK7oYm5y2oBFKKSs2J9Ft7dlxOEZGLZgwnvHfhWku02A6PibJvklgUrN7YbpQqgeV8LC1ZuVIBKztKTSJIT6hua0ioXyQUKUMkJZaXFaZWL5AIFqOSEaVUj6N+3qF1Z/75FTKsakaUWiRyaAlRywoSKQcy6qILjS4sx4PjS4s5H6q9bBE+Vw+N9gte6RT3cWpGABlSWnJHSYNN1i+C16QfmkG/cGiyD5k+SHqczUMkva289EJ6tWhqDcpEepgCV/NK4Lb1ykRgpQCW/lAxNr1wkRgpQyS+j5kBRSfuyopKgXKSHqRMpQ/Qcdw9p7Shae2tw2V4yNAhPdSBJFihAM0DPcfewYVMVmJITdAmfAV09xy0ivZcCNAP0HLdIYVKAZoCe4xYpTLoHmgHTqka0uwcKeo47Tuqwk1yhAM2A1j9e/VHHTx12kksUoBmS0nPc0m0aeFlyie6BSl5Rh53kEgWo5BV12EkuUYBKXtHAy5JLdA9U8oo67CSXKEAl76jDTnKFAlR6Lf1eVOKmAO2C/gDzl34vKj1BnUidaP0D3N3QhHPgD3BF7Y5sN01SoAFepCcoQDuhP8D8pt+LSk/QJXwn9AeY38pKi9md5LtK9ntR3aqRqAoyQFP5g0nnD1ByT6oDvKRzr1RBKx3FGqBmNgn4MVAEPOTud3dYb+H6C4FG4KvuvjqVbaNaUbuD1St/zL2lj1B23B7q9w/gsZVfBb7d7o9hWtUIVq/8MVeXPkLZYWG9hq9ydtW3k++4bpGmmcghqf5eNNVn6+MK2lTrFlq9bH92qmILUDMrAuYBE4HtwCozW+rubyRUuwAYHv4bAzwAjElx20j+8so8Zn7qPg7v0wzACX3rmfmp+3j4lcOYUPG/2upNOPL3jCv7KUXe1Fbvv5f9lKIjRwEdgrFuEbw2/cB85Y1bg2VQiGZRKr8XTfVWTVxBm0rdQquX7c9OR5ydSKOBTe6+2d0/BhYDkzvUmQw86oFXgWPMbGCK20ZyWfFDbeHZ6vA+zVxW/FD7imtvbQvPVkXeFJxldrT21gPh2aqlMXldySmpPlufiaDtKNW6hVYv25+djjgDdBDwVsLy9rAslTqpbAuAmU03s2ozq66vrz9ko44/bE9q5Y3bku8gWXk6dSWnpPpsfaaDNp26hVYv25+djjgD1JKUeYp1Utk2KHSf7+6V7l5ZVlZ2yEZ91O/E1MpLhibfQbLydOpKTplQMYhZF1VwfGkxBhxfWsysiyoOuqzLdNCmU7fQ6mX7s9MRZ4BuB4YkLA8GdqZYJ5VtIymuvIcWa/8/WosVU1x5T/uKo+ZAUUn7sqKSoLyjdOpKzplQMYjHbpjAM7ddxGM3TEh6TyzTQZtO3UKrl+3PTkecvfCrgOFmNgzYAVwJ/HOHOkuBmWa2mKATqcHdd5lZfQrbRjNsKkXQrse8KFmPeetyKj3r6dSVvJVKp1Q6o0WlWrfQ6mX7s9Nh7kmvjDPCzC4E5hL8FOlhd59jZjMA3P3B8GdMPwUmEfyMaZq7V3e27aE+r7Ky0qurq+M4FBEpYGZW4+6VB5XHGaA9TQEqInHoLED1LLyISEQKUBGRiBSgIiIRKUBFRCLqVZ1I4c+ftqaxyQAg+aNJ+UfHknt6y3GAjuUkdz/oSZ1eFaDpMrPqZD1r+UjHknt6y3GAjqUzuoQXEYlIASoiElGhB+j8bDcgg3Qsuae3HAfoWJIq6HugIiLdUehnoCIikSlARUQiKtgANbNJZrbRzDaZ2S3Zbk93mNkWM6s1szVmljejqZjZw2a228zWJ5R9ysx+Z2b/N3w9NpttTFUnx3KHme0Iv5c14QhjOc/MhpjZSjP7s5ltMLNvh+V59d10cRwZ+14K8h5oOGndmyRMWgdMycSkddlgZluASnfPqx86m9nngQ8I5sU6PSz7AfA3d787/A/bse5+czbbmYpOjuUO4AN3vzebbUtXOC/ZQHdfbWZHATXApcBXyaPvpovj+AoZ+l4K9Qw0tknrJHXu/iLwtw7Fk4GF4fuFBP+Hz3mdHEtecvddrdOLu/v7wJ8J5iTLq++mi+PImEIN0JQnrcsTDjxnZjVmNj3bjemmE9x9FwR/AMDxWW5Pd800s3XhJX5OX/ImY2blwFnAn8jj76bDcUCGvpdCDdCUJ63LE2Pd/WzgAuBb4eWkZN8DwCnAmcAu4N+z2po0mdmRwH8Bs9z9vWy3J6okx5Gx76VQAzS2Seuywd13hq+7gScJblHkq7fDe1et97B2Z7k9kbn72+7e4u6fAD8nj74XM+tLEDqL3P1XYXHefTfJjiOT30uhBmjbhHdm1o9g0rqlWW5TJGZ2RHiDHDM7AjgfWN/1VjltKXBt+P5a4NdZbEu3tIZN6EvkyfcSzlX2C+DP7v6jhFV59d10dhyZ/F4Kshceok1al4vM7GSCs04IZll9PF+OxcyeAMYTDC/2NnA78BTwv4GhwDbgcnfP+c6ZTo5lPMFlogNbgG+03kPMZWZ2HvASUAt8EhZ/l+D+Yd58N10cxxQy9L0UbICKiHRXoV7Ci4h0mwJURCQiBaiISEQKUBGRiBSgIiIRKUAlJ5lZeeLIRhne9xYzG5BG/bmtT3d1tq2ZXWxmd2aynZL7FKAiXTCzTwHnhoOFdGUZcImZlfRAsyRHKEAllxWZ2c/DsRyfM7NiADM7xcyeCQdPecnMTgvLv2hmfzKz183seTM7ISw/Ltz+dTP7GeFYCOFTXMvMbK2ZrTezK5K04cvAMx3Krjez1eEYrKcBePCD6t8DF8fyv4TkJAWo5LLhwDx3HwnsBS4Ly+cD17v7OcCNwP1h+csEZ4tnEQxR+D/C8tuBl8PypQRP0gBMAna6+6hwDM+OQQkwlmAcyUR7wsFbHgg/v1U18A9RDlTy02HZboBIF+rcfU34vgYoD0fW+Rzwn8GjzgD0D18HA/8RPuvcD6gLyz8P/BOAuy8zs3fD8lrgXjO7B/itu7+UpA0DgfoOZa2Da9S07je0GzgxrSOUvKYzUMllzQnvWwj+g98H2OvuZyb8+3RY5yfAT929AvgGcHjC9gc9s+zubwLnEATpXWb2/SRtaOqwn8R2tbap1eFhfSkQClDJK+F4jnVmdjkEI+6Y2ahwdSmwI3x/bcJmLwJTw/oXAMeG708EGt39l8C9wNlJPvLPwN+n2LxTyZMRlyQzFKCSj6YCXzeztcAGDkzHcgfBpf1LQOL8UHcCnzez1QTD/W0LyyuA18xsDXArMDvJZy0jGFUpFVVhfSkQGo1J5BDM7GXgYnff20WdEwiGEvzHHmuYZJ0CVOQQzGwM0OTu67qo81lgX0KnlxQABaiISES6ByoiEpECVEQkIgWoiEhEClARkYgUoCIiEf1/fIPnrxDLH4EAAAAASUVORK5CYII=\n",
"text/plain": [
"