{ "cells": [ { "cell_type": "markdown", "id": "above-final", "metadata": { "slideshow": { "slide_type": "slide" }, "tags": [] }, "source": [ "# Persister Cells" ] }, { "cell_type": "markdown", "id": "actual-public", "metadata": { "slideshow": { "slide_type": "subslide" }, "tags": [] }, "source": [ "There are several ways that bacterial cells can survive antibiotic treatment. The better known scenario is bacterial resistance, in which cells mutate and become resistant to antibiotics. Another scenario, which is also very common, is bacterial persistence. The difference is that persistence is not caused by genetic mutations, and hence is not fully heritable. What happens is that a small fraction of the cell population enters a dormant state, such that they stop growing and, in the mean time, are insensitive to antibiotics. There is obviously a disadvantage for the cell population to have these cells (called \"persisters\"), as they do not grow or divide, hence do not contribute to the growth of the population. However, when antibiotics strike, these cells will survive, and will later exit the dormant state and grow to recover the cell population. Thus, persister cells act like an insurance policy for the bacterial population --- a certain amount of asset is set aside in case an unexpected disaster happens and everything else is lost. In layman's language, \"don't put all eggs in one basket\". In the finance jargon, it is called a \"bet-hedging\" strategy, as if you are betting on horses or stocks and you don't want to bet all your money on one option." ] }, { "cell_type": "markdown", "id": "ebfed11d-df85-4a47-8f41-551f97f3b8b1", "metadata": { "slideshow": { "slide_type": "subslide" }, "tags": [] }, "source": [ "It is important to note that the persister cells are not a special subset of the cell population that are designated to be dormant. Rather, every cell has a certain probability of entering dormancy (at least for type-II persisters; type-I persiters are induced by stress). The question is, how many persisters should there be in a population? It is not always good (for the bacteria) to have more persisters. That is because, when there is no antibiotics, too many persisters means the population is not growing as fast as it can. On the other hand, too few persisters will not be able to protect the population from antibiotics. The right balance depends on how uncertain the environment is, i.e., how likely the bacterial population would encounter antibiotics. We will use a simple model to illustrate this point." ] }, { "cell_type": "markdown", "id": "4bbffb7f-9b4f-44ec-8dab-a0d0c9e01c47", "metadata": { "slideshow": { "slide_type": "slide" }, "tags": [] }, "source": [ "## Population growth in fluctuating environments" ] }, { "cell_type": "markdown", "id": "cutting-edward", "metadata": { "slideshow": { "slide_type": "subslide" }, "tags": [] }, "source": [ "For simplicity, we will consider discrete time steps. In each time step, the environment can be in one of two conditions, either with or without antibiotics. Let us denote them by $\\varepsilon = 0$ (no antibiotics) and $1$ (with antibiotics), respectively. Suppose that at every time step, the cell population always has the same fraction of persisters. Denote the cell types by $\\phi = 0$ (normal) and $1$ (persister), and the proportions of normal and persister cells by $q_0$ and $q_1$. The growth of cells is determined by a fitness matrix, $f_{\\phi\\varepsilon}$, which represents how many cells will be produced by a cell type $\\phi$ under environment $\\varepsilon$ within a time step. Thus, for example, a normal cell in an antibiotics-free environment will give rise to $f_{00}$ cells in the next time step. Let us assume that normal cells all die under antibiotics, so $f_{01} = 0$. Furthermore, persister cells can survive antibiotics but do not grow, so $f_{10} = f_{11} = 1$." ] }, { "cell_type": "markdown", "id": "3c0fdafe-b1b4-4d86-b5f8-ce9d5203da00", "metadata": { "slideshow": { "slide_type": "fragment" }, "tags": [] }, "source": [ "Now suppose the current population size is $N_t$. Then there are $q_0 N_t$ normal cells and $q_1 N_t$ persisters. If the environment is $\\varepsilon_t$, then at the next time step the population size will be:\n", "\\begin{equation}\n", "N_{t+1} = N_t q_0 f_{0\\varepsilon_t} + N_t q_1 f_{1\\varepsilon_t} = N_t \\sum_{\\phi} q_\\phi \\, f_{\\phi,\\varepsilon_t}\n", "\\end{equation}\n", "We will use this formula to simulate the population growth in randomly changing environments. For simplicity, we assume that at every time step, the environment is randomly chosen from $\\varepsilon = 0$ and $1$ with probability $p_0$ and $p_1$, respectively." ] }, { "cell_type": "markdown", "id": "coupled-young", "metadata": { "slideshow": { "slide_type": "subslide" }, "tags": [] }, "source": [ "Let us first write a python class to simulate the population growth in changing environments." ] }, { "cell_type": "code", "execution_count": 1, "id": "loving-distinction", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "id": "current-kinase", "metadata": { "slideshow": { "slide_type": "subslide" }, "tags": [] }, "outputs": [], "source": [ "class BetHedging:\n", " \"\"\"\n", " simulate bet-hedging strategy in changing environments.\n", " \"\"\"\n", " \n", " def __init__(self, phe_dist, fit_mat, N0=1., record=True):\n", " \"\"\"\n", " initialize population.\n", " inputs:\n", " phe_dist: 1-d array-like, phenotype distribution q[i] (should sum to 1).\n", " fit_mat: 2-d array-like, fitness matrix f[i,j] = number of offspring for phenotype i in environment j.\n", " N0: real, initial population size (use real instead of integer).\n", " record: bool, whether to record history of environment and population size.\n", " \"\"\"\n", " self.phe_dist = np.asarray(phe_dist) # phenotype distribution q_i\n", " self.fit_mat = np.asarray(fit_mat) # fitness matrix f_{ij}\n", " self.pop_size = float(N0) # current population size\n", " self.time = 0 # current number of generations since the beginning of simulation\n", " self.record = record\n", " if record:\n", " self.env_hist = []\n", " self.pop_hist = []\n", " \n", " def grow(self, env_seq):\n", " \"\"\"\n", " simulation population growth in changing environments.\n", " inputs:\n", " env_seq: 1-d array-like, environment sequence, list of environment indices.\n", " \"\"\"\n", " T = len(env_seq)\n", " for t in range(T):\n", " env = env_seq[t] # environmental condition at each time step\n", " factor = np.dot(phe_dist, fit_mat[:,env]) # growth factor according to formula\n", " new_pop = self.pop_size * factor # new population size at next time step\n", " if self.record: # record history\n", " self.env_hist.append(env)\n", " self.pop_hist.append(self.pop_size)\n", " self.pop_size = new_pop # update population size" ] }, { "cell_type": "markdown", "id": "illegal-influence", "metadata": { "slideshow": { "slide_type": "subslide" }, "tags": [] }, "source": [ "Let us choose some parameter values. We already have $f_{01} = 0$ and $f_{10} = f_{11} = 1$; let us choose $f_{00} = 3$. Let the environment distribution be $p_0 = p_1 = 0.5$, which means both environments are equally likely. Then let us set the phenotype distribution to be $q_0 = 0.3$ and $q_1 = 0.7$, so that there are 70% persister cells. (These numbers are not realistic and only for illustration.)" ] }, { "cell_type": "code", "execution_count": 3, "id": "driving-journal", "metadata": { "slideshow": { "slide_type": "fragment" }, "tags": [] }, "outputs": [], "source": [ "env_dist = np.array([0.5, 0.5]) # environmental distribution\n", "phe_dist = np.array([0.3, 0.7]) # phenotype distribution\n", "fit_mat = np.array([[3., 0.],\n", " [1., 1.]]) # fitness matrix" ] }, { "cell_type": "markdown", "id": "latter-interim", "metadata": { "slideshow": { "slide_type": "fragment" }, "tags": [] }, "source": [ "Now we generate a random sequence of environments and let the population grow." ] }, { "cell_type": "code", "execution_count": 4, "id": "based-broadcast", "metadata": { "slideshow": { "slide_type": "fragment" }, "tags": [] }, "outputs": [], "source": [ "T = 1000 # number of time steps\n", "env_seq = np.random.choice(2, size=T, p=env_dist) # choose random environment sequence according to distribution\n", "\n", "N0 = 1. # initial population size\n", "bh1 = BetHedging(phe_dist, fit_mat, N0=N0, record=True) # create bet-hedging population\n", "bh1.grow(env_seq) # grow population under given environmental sequence" ] }, { "cell_type": "markdown", "id": "injured-portsmouth", "metadata": { "slideshow": { "slide_type": "subslide" }, "tags": [] }, "source": [ "Let us plot the population size on log scale." ] }, { "cell_type": "code", "execution_count": 5, "id": "demographic-designer", "metadata": { "slideshow": { "slide_type": "fragment" }, "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEJCAYAAABhbdtlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXhU9fX48ffJToCwhLAkrEEWwyYSARUFNwSXYl2hlVpLXaq2fmt/tW6tVrSutdVqqyhoq3Wrda0gKKggIhBQloQ9CIQECASSANlzfn/cyTBJJskEMpmZ5LyeJw9zP3OXc0MyJ/eziqpijDHGeBMW6ACMMcYEL0sSxhhj6mRJwhhjTJ0sSRhjjKmTJQljjDF1siRhjDGmTpYkjDHG1MmShDHGmDpFBDqAKiJyGXAx0BV4TlUXiEgYMBOIA9JU9Z/1naNLly7at29fv8dqjDEtyapVq/araoK39/yaJERkDnAJsE9Vh3qUTwKeBsKBl1T1UVV9H3hfRDoBTwILgClAEpAHZDV0vb59+5KWltb0N2KMMS2YiOyo6z1/Vze9AkyqEUw48BwwGUgBpolIiscu97neBxgELFPVO4Bf+DlWY4wxNfg1SajqYpynAE+jga2qmqmqpcCbwBRxPAbMU9XVrn2zgIOu1xXeriEiN4pImoik5ebm+uEujDGm9QpEw3USsMtjO8tV9kvgfOBKEbnZ9d67wIUi8jdgsbeTqeosVU1V1dSEBK9VasYYY45TIBquxUuZquozwDM1Co8CM5olKmOMMbUE4kkiC+jlsd0TyA5AHMYYYxoQiCSxEhggIv1EJAqYCnwYgDiMMcY0wK9JQkTeAJYBg0QkS0RmqGo5cBswH9gAvK2q6f6MwxhjzPHxa5uEqk6ro3wuMNef1zbGmNbiozXZHDpayvTT+zb5uW1aDmOMCXFvrtzJf1Y1ON74uFiSMMaYEJaRXcDSrQdI6RHnl/NbkjDGmBC1aU8hFz2zBID+Ce38cg1LEsYYE2JyC0sAWLz52CwTI3t39Mu1gmYWWGOMMQ1btu0A0178hstPTeLd1bsBeGH6KFL7dvbL9exJwhhjQkja9850eFUJAuDCId39dj1LEsYYE0K+P3C02vbs61L9ej1LEsYYE0IycgqqbZ93cje/Xs/aJIwxJkSUlFewZW8ht0zoz1kDEigp97qCQpOyJGGMMSFiy97DlFcqKYlxnN4/vlmuadVNxhgTIjKynaqmIYkdmu2aliSMMSZEpGfn0zYqnD6dY5vtmpYkjDEmBGTmHuafy3bQL6EtYWHe1m7zD0sSxhgTAh76eAMA+wtLm/W6liSMMSYElFcqAI9dObxZr2tJwhhjQsD2/Ye5eHgPxg9MaNbrWpIwxpggtC33MNNnLyft+zzSs/PZlVfkt+nA62PjJIwxJgj9beEWlmzZz5It+91lQxKbP0nYk4QxxgShsgqtVdac4yOqWJIwxpgAWbUjz702BMCBwyXk5BcBsHFP9TmaEtpHk9A+ulnjA0sSxhgTEKXllVzxj2Vc+fzXVFYq6dn5jHroM05/ZBFHS8vJ3H+EX503gDV/mEiYBKaqCYKsTUJELgMuBroCz6nqAm9lgYzRGGOawrurswDYceAoj8zbwItLtrvfu/3N71CFoYlxdIiN5JfnDmB4z+avagIQ1dr1Xk16AZE5wCXAPlUd6lE+CXgaCAdeUtVHPd7rBDypqjPqK6spNTVV09LS/HAXxhjTdFSVfnfPbXC/pXedS1LHNn6PR0RWqarXhSmao7rpFWBSjYDCgeeAyUAKME1EUjx2uc/1Pg2UGWNMUCstr6SotAJVZVvuYa5+fhnvfbu74QOBxA4xfo6uYX6vblLVxSLSt0bxaGCrqmYCiMibwBQR2QA8CsxT1dWu96RmmTHGhIoHPkpn9Y6DbNxT6C5b4VqC9PErhvPw3A3kF5UB8Py1p3Lza87H3Ae3nonz8RdYgWq4TgJ2eWxnucp+CZwPXCkiN7ve81bmJiI3ikiaiKTl5ub6OWxjjGmc5ZkHqiUIT5OHdWdMv84AvDZjDJOG9uCeiwYDMKh7+2aLsT6Barj2lh5VVZ8BnqlRWKusxvuzgFngtEk0ZZDGGHMiqnopefPbCwfRPiaScwd3ZcmW/QxNcnov3XBWMjeclRwUTxEQuCeJLKCXx3ZPIDtAsRhjTL1Kyit4fflOyisqG3XchpxCPPsG/XhMb3p1dhqiLxneA4BrTuvF13edS8fYKABEJGgSBATuSWIlMEBE+gG7ganAjwIUizHG1OvVZTt46OMNVKpy7dg+de63/3AJv37rO2495yTi20aRkeMMiDulV0eSOrXh4R8Oo6JSycw9TJ/4toCTFDq1jWqW+zgefk8SIvIGMAHoIiJZwP2qOltEbgPm43SBnaOq6f6OxRhj6pKenU9Kj7haf8VXVCpPfboZqD0KuqalW/fXmm8J4L1bznCfNzxMGNAtONobfOH36iZVnaaqPVQ1UlV7qupsV/lcVR2oqv1V9WF/x2GMMd7kHSnlN2+v4eJnvuLVb3YA8PnGfZSWO1VLn23Yy9HSCld5LvlHy7yeR1VJz/aeRIKp+qixgmrEtTHGNLd73l3HJ+l7APhiUy5Pf7aFA0dKmTGuH/ddfDJPzt/k3nf3oSJGPLiAlB5xvHz9aXSLc8YxfL11Pze+uorEjrXHNfz1mlOa50b8xOZuMsa0arsPFblf7zhwhANHnOVBl207wIacQrbsOwzAqzNGu/fLyClgzJ8Wku069qut+zlcUs7mvYdrnf+ykUn+DN/v7EnCGNOqVQ1kA9iWe6y76p6C4mptEGP6xdc69uO1OdxwdnK1aqYHLk1hZO9OrNpxkHYxof8Ra08SxphW6+nPtrAz7yiXjkjk6anHqoV+ODKJvCOlfPCd0zN/3QMTiYoIY9NDk7hyVE+iI5yPzg2u3kueSWJoUgdG9OrIz8b14+pUz57+oSn005wxxvjo47U5rPw+j56d2pCRXcC7rjmUJg/tzsSUbtz+5ncA3HrOSbz37W6+3JzLiJ4daB8TCUB0RDhPXjWCRy8fxs2vrebdb3czoldH9h8uYUSvjuw8cISTA7DEqD9ZkjDGtBq3vu59+rezByYQER7G9kcuQkRQVXp0iCEnv5gUL6vBRYSHMX5QAp9t2Mv9Hzq99++ePJjRfTsTFha6PZm8seomY0yrsHiz97ndnrxqBO2inb+Xq7qqighjk502iLoW+zmtb6dq2ymJcS0uQYAlCWNMC7cr7yjXvLCMn8xZUa28auhCXUng9P5OkhiW5H2xn8Hd47j9vAEARIQJca4qqZbGqpuMMS3W+9/u5v/e+q5a2ae/PpuTurbjnvfW8d63uzmpazuvx14+Momu7aPrXRHu1xcM5IpTe7rnY2qJLEkYY1qsmgni9OR495QY/2/iIKae1pvIcO8VKhHhYUwY1LXBa/SOjz3xQIOYJQljTIsVHiZUVB6bhvWei052v45vF018u+hAhBVSLEkYY1qkg0dKqahUbjw7GQF+M3EQURHWDNtYliSMMS1S1TTdZw9IYNyALgGOJnRZWjXGtEjp2flA3b2XjG8sSRhjmlRZRSVHSsoDHQbp2QUkdogJ6gV9QoElCWNMkykoLmPAvfMYcv98VAO35PxTn27mg++ySbGniBNmbRLGmBOmqsz+ajt5rmm2wZlRta4xCP60r7CYZxZuAfA6pYZpHHuSMMacsJz8Yh76eAOzFme6y77JPFDP/kX8b212k12/slK55711rM06VG1G1pQWNtleINiThDHmuOQdKaW0vJLuHWLcH8zllcqAru04cKSU+95fT/uYCBLaR3NG/2O9iwqKyzj9kUWAk0geumwYAMVlFfx7+U6mj+3T6K6q2w8c4fXlO6moUJI6HRv9bI3WJ86eJIwxjTbnq+2cOvNTxj6yEIAMj7/ehyTGEedabOf2N7/jRy8ur9aQvcFj39e+2UluYQkA89P3MPN/GSzcsLfR8VQlqbfSdvHUp5sBWHLnOfTq3LJHQzcHSxLGmEapqFQe/F+Ge3tfQTHp2flEuaa3SEmMY+ro3tWOWbXjoPu1Z3UQwPLtTrVUVaJZvj0PgJ0HjnLdnBXsPlTEba+vJuvg0TpjyqhxTsASRBOxJGGMqVNlpfL0Z1vcH9CqWusD+ZvteWTkFHDBkG7MvGwoV6f24qazk5kxrt+xfTzaJzJyCujSLor//uKMau9VDX6r2v7Hl1v5cnMuZz66iP+tzeFPczfUGWd6dn61Kqr3bjnjRG7beAiaJCEiySIyW0Te8SjrLSIfisgcEbkrkPEZE0wOHS0lM/ew36+zcU8hf/lsM+Me+5yV3+dx7/vrufTZrwD498/H0D46gvnpe8g6WMTQxA5MH9uHjrFRiAg/GJHoPo9nkkjPLiAlsQOj+nRiwqAEvtiUy9HSctKzC4gIEzbuKSTvSClFpRXVYtmTX+w1xqrENWVEIq/OGM3GmZMY2buT131N4/k1Sbg+3PeJyPoa5ZNEZJOIbK368FfVTFWdUeMUA4GPVfVnQIo/YzUmVOwrKOaUBz/l3D9/WW3yOn/YuOfYU8NVzy/j9eU73duj+3XmtH6d+XhtDkCtMQkjenVk4W/Gc8uE/qzeeYgb/5VGSXkFW/YWunsdjU2OJ+tgESl/mE/ekVIuGtYDgBXbD7BxT2G1823IKaSsorJWjPsKSzhwpJQhiXGcNSCBmMjwprl5A/j/SeIVYJJngYiEA88Bk3E++KeJSF0J4FtgqogsAj73Y5zGhIz5Gccadr3VxTel+s4fGR5GqsfqbN56EvVPaOdevGdBxl4ysgsor1T3vlWrv1WZOroXADe/tpqNewq5eXx/rji1J8kJbSkqq2BtVn6ta1RNv2FjIvzDr0lCVRcDeTWKRwNbXU8OpcCbwJQ6TnE9cL+qngtc7G0HEblRRNJEJC031/vyhMa0JFl5xxpw6xuLUJ+P1+ZUG/hWl/TsAgZ3b8+IXh2rlT9+xXAARvft7C7rUse026P6HEskb6xwnkSqksTQGoml5ipww5I68OerR/Cfm04HvN9vVSI7uUf7Bu/HNF4g2iSSgF0e21lAkojEi8jzwEgRudv13ifAr1zl33s7marOUtVUVU1NSEjwZ9zGBIWMnAKGJsWR3KXtcSWJ7ENF3Pr6ak6d+Sl7C5x6/spK5XCN+ZZUlYycAkb27sS7vziDjTMnMfOyoXz1u3O4+jTnL/7Uvp35+FfjyPzTRXVeLzYqgqV3nQvA22lZAPSNbws4C/tsnDmJqPAw+sbH0j4mkh+POdYzqiqZxLeLZlC39l7vNz27wH2saXqBGEznbaVwVdUDwM01CtcDVzZLVMaEAFUlPbuAC07uRliY8L812VRUKuFh3n6taqusVPc4AoAxf1rI2zedzpyvtvNJ+h4yHryQyPAw/rZoK93jYsgvKmNIYhzhYUJ4WDjTx/apdc4hPlTzJHWsvrxnmEe8MZHhXDoikYT2zpPIAz8YQm5hCQsy9tLboxvr2OTOvJ2WRVlFpXs1ubVZh5i3fg+Th3b36f5N4wUiSWQBvTy2ewJNNz7fmBbsH19uI+9IKSmJcXSMjeSNFTvJyC5gWD3rMHv6zX/W8N63u6uVXf3CMvfr1TsO8cWmfbz01XZ3WVONWh7RswNrsvKrPSlU+fPVI9yvI8PDeGH6KCq1ejIZmxzPP5ftYG1WvrsK69qXlgPYYkJ+FIgksRIYICL9gN3AVOBHAYjDmJDz+CebAKcnUdVf2d9kHvA5SdRMEDV9k3mAuety3NthAoO7N02SeOX60WzeW8iYGo3V3ogI4TUejkb36+yOcVSfTmTmHqag2Kki++kZfZskRlObv7vAvgEsAwaJSJaIzFDVcuA2YD6wAXhbVdP9GYcxoU5V+XLzsY4ZKT3i6BYX06h2icpKpeoP8+X3nMeGByfRJ776qORnP99Ktsd4hOSEdrSJapoupZ3aRvmUIOpSs13i9x84Pev7J7S1cRF+5O/eTdNUtYeqRqpqT1Wd7Sqfq6oDVbW/qj7szxiMCXUHDpfw0docrpuzAoDnrx1F22inEmBMcjwrtufVGi+hqvz9i63syjvKhpwCzn/qS1bvPEilwiOXD6NbXAxtosLdVT+Xjkhk/MBjHT9uGp9MVERY0E2QNza5M2nfHyTr4FGWbnWSxX0X2xAqf7JZYI0JYnlHShn10GfVyjw/uMcmd3a3SxSXV/DYvI08esUwps9eQU5+MY9/somrU3uydd9h/v7FtlrH/3xcMqf17czI3p34ZH2O+2lleFJHZk0fRR9XL6RgUdUuccu/V7vLzhncNYARtXwNJgkRGQj8A+imqkNFZDjwA1V9yO/RGdNKbcgpYPLTS+javvbYg54eU2FXDUa79Nmv6No+mn2FJfzitdXkeFQZVXU7XbRxHwADux0bTxAWJu6qmklDezDzsqH8/v31jOrTie4dYpr+xk5QVbvEut3OALqnPBq8jX/4Ut30InA3UAagqmtxGpuNMX7yw78vBZwpJzz9/cenInKsRbdb3LEP8qp9tzYwp1N901ZMH9uH7x+9OCgTBDjtEgO7tUMVxp3UhctP7RnokFo8X5JErKquqFEW+FXOjWnBKmtMUTTupC4sv+c899xGnk6usfpa1dLSXdtHM2GQ084QE+n8qv/q3JOaPthmVvX0ZOtXNw9fksR+EekPKICIXAnk1H+IMeZ4lVVUgsC5rrr2p6eewms/H1PtqcHTGzeMcY8TOM01l9KMcf345u7z+NmZznTd//jxKFbcex53TBzUDHfgX1VJItga1VsqXxqubwVmAYNFZDewHbjWr1EZ00rlFpZw5ztrKC2v5NIRPZjz09MaPKZjbBRr759IaUUl3+08xE/mrGBYUgfCwoSzBnTh9RvGcHpyfLVqqlB23slduXvyYCam2Cjr5tBgklDVTOB8EWkLhKlqYUPHGGOOz1OfbuLzTU4PI1+mu6gSExlOTGQ4407qwtNTT2HyUKdaSkSqrS/dEkRHhHPT+P6BDqPVaLC6SUQqRORR4GhVghCR1Q0cZoxppMLiMt5YcWzuy+Quje9+GhYmTDklyaapME3Gl5+kdNd+C0Skal7glvHcakwQ8VzQB5wZUo0JNF9+CstV9U6crrBLRGQUrkZsY0zT2X/Y6cJ6TWovVt57foCjMcbhS8O1AKjq2yKSDrwB1J7G0RhzQjbkFDI0KY7Hrhwe6FCMcfPlSeLnVS9cE/GNA37lt4iMaYWqFvhJ6WHdOk1wqfNJQkTOVdVFQB8RqbnSSP1DOo0xjbKnoJi8I6WN6tFkTHOor7ppPLAIuNTLewq865eIjGmBVLXaOIXt+49w3/vruGvSyXSMjWTzXqdnuY0iNsGmziShqve7/r2++cIxpuX53TtreSttF9ERYZSUV/LmjWOZOusbwJmYD+DX5w8Eak+xYUyg+TJO4nYRiRPHSyKyWkQmNkdwxoQ6VeWtNGfsQ0m5MyHTTa+uqrXfB9/tpm98LO2ibfZ+E1x8abj+maoWABOBrsD1wKN+jcqYELZ5byH7CospLa+k391za72fX1QGQITH+s2Z+49Ye4QJSr4kiaqf5IuAl1V1DTaYzrRw+UVlqDZ+OJCqMvEvixn98ELeWZVV7b2rRlWf1vrJq6qvhWDtESYY+ZIkVonIApwkMV9E2gOVDRxjTMjKLSxh9MOf8dHa+ic7rqxUPlyTTd6RUnfZN5l57tf3vLfO/XrdAxN54qoRXD4yCYDXbxjDlFMSefEnqVx+qlNmScIEI18qQGcApwCZqnpUROJxqpyMaVEysgv4dtdBurWPoaS8kk/W51BUWs6Vo3oRHlb94TnvSCmnzvzUvb38nvPIyCng+pdXVtsvvm0Uafed7+7Z9MRVI/jT5cPcC/9ckNINAT5N38uInh39e4PGHAc5nkfqYJWamqppaWmBDsOEqGEPzKewuPZ6Wh1jI1ly5znM/mo7N5yVzJebc6utsey536GjZdVenzWgC6/OGNPgtSsrlbAwq8U1gSEiq1Q11dt7QdOVQkSSgXuBDqp6pUd5W2AxcL+q/i9Q8ZmWz1uCADh0tIzb3/yORRv3UVmpPPv51jr3A/jbtJEM6t6eiX9Z7HMVkiUIE6z8Os2kiMwRkX0isr5G+SQR2SQiW0XkLnDWrVDVGV5O8zvgbX/GaVq3V5d9z8z/ZdS7z6KN+wB4ZtFWKj0evtf8oXpv8B+N6c2lIxIZ0LUdd08ezI9H15yswJjQ4lOSEJFwEUkUkd5VXz6e/xVgUs1zAc8Bk4EUYJqIpNRx3fOBDGCvj9czplFOeXABv/8gndlfbQfglgnOYja3neOsBZ3UsU2dx/7q3JPoEBtJxoMXckb/6ktqigg3je9P7/hYf4ZvjN81WN0kIr8E7sf5oK7q1aRAg1NVqupiEelbo3g0sNW14h0i8iYwBScZ1HQO0BYnmRSJyFxVtZ5V5oQUl1UQGR5GQVGZu4qoyjWn9eIXE/rTNiqCGeP60SYqnGtfWk7ajoPV9nv+2lFMGuosnxkbFcHU0b35etsBa3w2LY4vbRK3A4NU9UATXTMJ2OWxnQWMcfWaehgYKSJ3q+ojqnovgIj8FNjvLUGIyI3AjQC9e9sM5qZ+qsrg339Cp9hIDtZIEO1jIujdOdbdE6lT2ygAXpg+iq+3HWDCoASGPbAAoNZsrZcM60GvTm0YmmQD4kzL4kuS2AXkN+E1vbXQqSsJ3eztAFV9pa6TqeosYBY4vZuaIkDTcmUdLAKoliAGdmvHaz8fQ0K76GqT8FWJbxfNpSMSAdg4cxJrdh2qVY0UFiaM7N3Jj5EbExi+JIlM4AsR+RgoqSpU1aeO85pZQC+P7Z5A9nGey5hGSc+u/vfO2OTOvHnj6T4fHxMZzpjk+KYOy5ig5UuS2On6inJ9naiVwAAR6QfsBqYCP2qC8xrToPTsAgB6dW7DM1NH2l//xjSgwSShqn8EcE3Hoarq84JDIvIGMAHoIiJZOGMdZovIbcB8IByY41rxzhi/y8guYGC3diz49fhAh2JMSPCld9NQ4FWgs2t7P/ATXz7YVXVaHeVzgdrTYxrjRws37GXhxn1cMrxHoEMxJmT4Ut00C7hDVT8HEJEJwIvAGX6My5gmsye/mIfnbuCjNU7TV3REeIAjMiZ0+DKYrm1VggBQ1S9wxi4YExJmLc50JwiA288bEMBojAktPvVuEpHf41Q5AVwLbPdfSMY0rcMlx7q7DkmMs1HQxjSCTyvTAQnAu8B7rtc2VbgJGRv3FLpfP35lgxMFGGM8+NK76SDwq2aIxZgmV15RycY9hfzk9D7cNXkwsVFBM/GxMSGhzt8YEfmrqv6fiHyEM1dTNar6A79GZkwT2JZ7hNLySk7t3ckShDHHob7fmqo2iCebIxBj/KFqhLUtDWrM8akzSajqKtfLU1T1ac/3ROR24Et/BmZMU8jILiA6IozkLtYhz5jj4UvD9XVeyn7axHEY4xfp2QUM7t6eiHC/rq9lTItVX5vENJw5lfqJyIceb7UHmmracGP8RlXJyCngomE2wtqY41Vfm8TXQA7QBfizR3khsNafQRnTFHYfKiK/qMy9WpwxpvHqa5PYAewAfJ9H2ZggkuGa8dUarY05fg1W1IrIWBFZKSKHRaRURCpEpKA5gjPmRKRnFxAmcHJ3SxLGHC9fWvOeBaYBW4A2wM+Bv/kzKGNO1L7CYp5euIW+8W1pE2UT+hlzvHwaXaSqW0UkXFUrgJdF5Gs/x2XMcausVEY/vBCAfYUlDextjKmPL0niqIhEAd+JyOM4jdnW6dwErZ15R92vrzujTwAjMSb0+VLdNB1nBbnbgCM461Nf4c+gjGkMVeWT9TmUV1QCkJHjNJlNH9uHOy4YFMjQjAl5DSYJVd2hqkWqWqCqf1TVO1R1a3MEZ4wvvt52gJtfW83c9XuoqFR+947TQ/vei08mPEwCHJ0xoa2+wXTr8DKxXxVVtTmXTVBYt9uZn2nZtgO0iw6nsKQcgJhIa7A25kTV1yZxSbNFYcwJSHeNh1i8OZeendoAMLh7+0CGZEyL0dBgOmOC2lMLNrmXJt19qIgn5m8iNiqcD247M8CRGdMy+DKYrlBEClxfxTaYzgSLvQXFPLPIaR7r3DbKXT6wW3uiI6yqyZim4EvDdXtVjXN9xeD0bHq2qQMRkWQRmS0i73iUtRWRf4rIiyLy46a+pgltVWtFAPzz+tHu1xHWWG1Mk2n0/Mmq+j5wri/7isgcEdknIutrlE8SkU0islVE7nKdN1NVZ9Q4xeXAO6p6A2Ar4ZlqquZmWvzbcxjWswObH5rM5acm8dAPhwY4MmNajgYH04nI5R6bYUAq9fR6quEVnKeOf3mcLxx4DrgAyAJWisiHqprh5fiewDrX6wofr2laifTsAvrGx9I7PhaAqIgwnrr6lABHZUzL4suTxKUeXxfiTBU+xZeTq+piIK9G8Whgq+vJoRR4s57zZeEkijpjFZEbRSRNRNJyc3N9CcuEkG25h7n51VUcLS2vVl5UWsEXm3Jthldj/KzBJwlVvb6Jr5kE7PLYzgLGiEg88DAwUkTuVtVHgHeBZ0XkYuCjOuKbBcwCSE1N9fUJx4SIeety+CR9D1dt60mlQnJCW7IPFTF99goABtsMr8b4lS/VTcnA08BYnGqmZcCvVTXzOK/prVVRVfUAcHONwiNAUycpE0KqxkAs2bKfV77+vtb7VeMijDH+4Ut10+vA20APIBH4D/DGCVwzC2f+pyo9gewTOJ9poUrKK/h6m7NSrrcEMaJnB6acktTMURnTuviSJERVX1XVctfXa/jecO3NSmCAiPRzzS47FfiwgWNMK7M88wCD7vuE/KIy2kfXfuB9cMoQPrhtnM3NZIyf+ZIkPheRu0Skr4j0EZE7gY9FpLOIdK7vQBF5A6d6apCIZInIDFUtx5lRdj6wAXhbVdNP9EZMy3LNrG/cr8cPSnC/HtTNmW7jgpRuzR6TMa2RL+tJXOP696Ya5T/DeaJIrutAVZ1WR/lcYK4vAZrWZ23WIffrfl3a8tTVp5CSGMdFQ3vQu3Msuw8V0aODtUUY0xx86d3UrzkCMabKD55dCsA9Fw3mxrP7A3DLhJPc7/fqHAoGH7sAABRLSURBVBuQuIxpjXzp3RQJ/AI421X0BfCCqpb5MS7TSlVWHmvuumR4YgAjMcaAb9VN/wAigb+7tqe7yn7ur6BM63HoaCn7D5dyUtd2vL1yF8XlzsD6n4/rR2JHq1IyJtB8SRKnqeoIj+1FIrLGXwGZ1mXai8vZkFN7UmHr2mpMcPCld1OFiPSv2nANrrN5lEyT8JYgAAZ0a9fMkRhjvPHlSeK3ON1gM3FGS/fBRkGbJlBaXum1vFNspC09akyQ8KV300IRGQAMwkkSG1W1xO+RmaCiqqhCmMfgtayDR/nrZ1uYOWUobaIa/6G+ZV8hAGP6daZDm0gWZOzlo9vGMaxnhyaL2xhzYnzp3RQD3AKMwxkXsUREnlfVYn8HZ4LH2U98Tv+Edrxy/Wjyi8o4+/HPyS9yOrideVI8PxzZs4Ez1LZ6x0EA/nT5MPonWPWSMcHIlzaJfwFDgL/hrA2RArzqz6BMcCkuq2BXXhFfbMrllaXbGfHHBe4EAbB+d+NXs91XUMzvP3AG2veNb9tksRpjmpYvbRKDavRu+tx6N7UeOflFnP7IIvf2Ax/VXhvqm8wDjT7viu+PLTNi8y8ZE7x8eZL4VkTGVm2IyBhgqf9CMsHk8U82eS2/alRPBndvT9f20WTkFJB/1PexlcVlFdz2+rcAvPST1CaJ0xjjH748SYwBfiIiO13bvYENIrIOZx2I4X6LzgRc1YpwA7q2Y8u+wwB8cOuZjOjVEXBma71m1jes+D7P50n3Nu0pdL8+3ybqMyao+fIkMQnoB4x3ffUDLgIuwVnS1LRgmblHOP/krnx6x3j6JzhtB8M9eh9VJYsb/pXGl5tzueCpL/l80756z1m1kNDMKUP8FLUxpqn40gV2R3MEYoJPcVkF23IPM3lodwDeu/VMCovLETnWhuA5nuHOd9awt6CEj9Zkc86grgBs33+E6bOXc2rvTvz56hFEhoeRnp1P++gIrh3bp3lvyBjTaL48SZhWavZX26lUSEl0nhziYiJJ8jKf0q/OGwDA3gJn+MzyzDxUlZ0HjnLOk1+QdbCID9dkc+c7awHIyCng5MS4asnGGBOcLEm0UuUVlby1cmedDc6qyhPznUbrIYlx9Z7rjgsGctPZx5YV2X2oiKyDRcxbn1Ntv/e+3c3zX27j252HGjynMSY4WJJohfbkF3PSvfP43X/X8Zv/rEFVef7LbWTmHqag2EkaWQeL3Pv37NTwbKxjk+MBSOwQA8CnGXt5ZN7GWvs96ipL6WFJwphQYEmilXlr5U7GPrLQvf3Zhr0s2bKfR+dt5Nw/f8nwBxawcU+Bu3F5zk9TfaoWSu3biajwMC4bmUTntlE89PGx8RS/uWAg//zZ6Gr7D0m0qTeMCQW+dIE1Lcjv/ruuVtlP5qyotj3pr0vcr09P7uLTedvHRPLhL8+kZ6dYtuUeZn76XgBe/ulpnDO4K0dLy2kTGU5RmTOB8EldbRoOY0KBPUm0cCXlFTwxfyOHjpbWeu+xK4Y1eHxjJu4b3D2OdtERJHU8trzohEEJAMRGRbBh5iRW3Xc+791yBlER9qNnTCiw39QW7NudBxl03yc89/k2npi/qVqiuO/ikxtc2OeWCf3rfb8uN41Ppm98LC9MH1Wrqiq+XTQje3c6rvMaY5qfqGrDewWIiIQBM4E4IE1V/1nf/qmpqZqWltYssYWCi59Z4m5b8PTqjNGcNcD5C7+sopIPvsvmByMSeSttF1NOSWT4AwsA2P7IRdZN1ZhWQERWqarXOXKavU1CRObgjNbep6pDPconAU8D4cBLqvooMAVIAvKArOaONdRFhHt/UPTsWRQZHsaVo5xpvqe7BrfNvGwoBUVlliCMMQFpuH4FZ8rxf1UViEg48BxwAU4yWCkiH+IsdLRMVV8QkXeAhbVPZ7xxBrMd8fpefLvoeo+dbiOhjTEuzZ4kVHWxiPStUTwa2KqqmQAi8ibOU8QuoKoi3dbVboSc/GIOHi3jwSlDuDq1F2+t3MXK7/PoGBsZ6NCMMSEkWLrAJuEkhCpZOLPPPg38TUTOAhZ7O1BEbgRuBOjdu7efwwwdv33HWfJjSGIcMZHhXHdGX647o29ggzLGhJxgSRLeKr9VVY8CM+o7UFVnAbPAabj2Q2whpbS8kvLKSpZudRYCGtzdRjYbY45fsCSJLKCXx3ZPIDtAsYSkQ0dLmbP0e55ZuIWzBhwbANc2Olj+i40xoShYPkFWAgNEpB+wG5gK/CiwIYWOvQXFjPnTsTb9JVv2A/DV784JVEjGmBai2QfTicgbwDJgkIhkicgMVS0HbgPmAxuAt1U1vbljC1ULMvZ6Lfc2rbcxxjRGIHo3TaujfC4wt5nDCXkFxWX8/v31tco/vO1MG+dgjDlhNi1HiFu946D79czLhrrHOAzs1j5QIRljWpBgaZMwx2nz3kIAfn9JSrUR08YY0xTsSSKElVdU8s+vd5DYIYYZ4/oFOhxjTAtkTxIh5nBJOaf8cQF3TR7MQx9vAKBtI6bzNsaYxrAniRCTvjuf8kp1JwjAPaOrMcY0NUsSISYjp/rU351iI3ny6hEBisYY09JZkggxGTXWh/jyznNoZ6OqjTF+Yp8uISY9u4CzBnTh5vH9Wbp1P3ExNqurMcZ/7EkiRKgqt/x7FRk5BaQkxnHmSV24c9LgQIdljGnhLEmEiAUZe5m7bg8APW26DWNMM7EkESJeWpLpfn3h0O4BjMQY05pYm0QISM/OZ+X3zvQbax+YaO0QxphmY08SIeCFL489RViCMMY0J0sSQS63sIQP1zjrL6Xdd36AozHGtDaWJIKc5+C5Lu2iAxiJMaY1siQR5NKz8wFY+JvxAY7EGNMaWZIIQh+vzSH7UBHgDJ7r2akN/RPaBTgqY0xrZEkiyOQfLePW11fz7OdbAdiQXUBKj7gAR2WMaa0sSQSZ9Byneun15TtZteMg2w8cYUhihwBHZYxprWycRJDxnMDvin98DUBKoj1JGGMCw54kgkzNWV4BhliSMMYEiCWJIHKkpJyP1mYzfmACpyfHA9AxNpIeHWICHJkxprUK6iQhIpeJyIsi8oGITAx0PP5282urKKtQ+ie04+XrTyMqIowhiXGISKBDM8a0Us2eJERkjojsE5H1NconicgmEdkqIncBqOr7qnoD8FPgmuaOtbkt2bIfgL5dYomJDOfei07mhrOSAxyVMaY1C8STxCvAJM8CEQkHngMmAynANBFJ8djlPtf7Qe+jNdnkHSl1b6/LyueHf1/Kbte4h7qoKh3aRJLQPpofj+kDwHVn9GXCoK5+jdcYY+rT7ElCVRcDeTWKRwNbVTVTVUuBN4Ep4ngMmKeqq72dT0RuFJE0EUnLzc31b/AN2H2oiF++8S0vLsmkqLSCg0dKufTZr/h25yHOfHQRu/KOVtu/oLiM619ewbx1OWTnF5NfVMavzj2J8DCrXjLGBIdg6QKbBOzy2M4CxgC/BM4HOojISar6fM0DVXUWMAsgNTVVmyHWOqXvdsY4vPDlNv67Kot9hSXV3n956fdMP70P/bq0pbC4jOEPLADg8025XJDSDbDursaY4BIsScLbn86qqs8AzzR3MI21K+8oURFhpLu6r1YqtRIEwJyl25mzdDsr7j2PrXsPV3vv04y9iMDg7pYkjDHBI1iSRBbQy2O7J5AdoFh8VlJewS9eW82ijfvq3W/Rb8bzxPxNzFvvLD/67urdPDpvIwCxUeEcLa0AoF+XtrSNDpb/EmOMCZ4usCuBASLST0SigKnAhwGOCXCqjm7596pa5XlHShl03ydeE0R7jw/6uJgI+nVpy3knd3OXVSUIgIwHJ/H4lcMBbI4mY0zQafY/W0XkDWAC0EVEsoD7VXW2iNwGzAfCgTmqmt7csXkzd/0e1uw6xMEjpXRqG+Uu/++qLK/7v3z9aZx1Uhfe/XY32/cf4WhJOSLCmH6da+07fazTi2liSjfujwwntU8n/9yEMcYcp2ZPEqo6rY7yucDcZg6nTgXFZcz6MpM1uw4BsHx7HkOT4vjtf9Zy+/kDWL79AAD9E9ry7i1nMntJJlel9qJX51gArk7tVe18vTrH8sr1p7G3oJjf/XcdAH+41Onl2zE2ii9+O4F4jyRkjDHBwCrAvdi8t5CJf1lcreybzAOkfZ/HsswDLJvlJIjhPTvw4W3jALhj4qAGz1s15qFbXAyb9xYSGX6stq9bnE29YYwJPpYkPOw+VMSZjy6qVd6jQwzz1uewt6B6j6Wkjm2O6zoTBnW1QXLGmJAQLA3XQcFbggCn6qhmggD7698Y0/LZk4RLkasbqqc1f5hIdn4RR0rKeXrhFgC6tItm/+ESfnZmP+6YOLC5wzTGmGZlScJl095CAKaN7k32oSKmje5Nh9hIOsRGUlpe6d5v+T3nsXlvISdbd1VjTCtgScIlPduZUuOWCf3dPZSqREWEces5/UloF014mFiCMMa0GpYkXDKyC4iLiaBnJ++N0b+9cHAzR2SMMYFnDdcu6dkFpNgCP8YYU40lCaCiUtm4p4CUHh0CHYoxxgQVSxLA4ZJyJg/twRn94wMdijHGBBVrkwA6tInkL9ecEugwjDEm6NiThDHGmDpZkjDGGFMnSxLGGGPqZEnCGGNMnSxJGGOMqZMlCWOMMXWyJGGMMaZOliSMMcbUSVQ10DE0GRHJBXacwCm6APubKJxQYffcOrS2e25t9wsnds99VDXB2xstKkmcKBFJU9XUQMfRnOyeW4fWds+t7X7Bf/ds1U3GGGPqZEnCGGNMnSxJVDcr0AEEgN1z69Da7rm13S/46Z6tTcIYY0yd7EnCGGNMnSxJGGOMqZMlCUBEJonIJhHZKiJ3BTqepiIivUTkcxHZICLpInK7q7yziHwqIltc/3byOOZu1/dhk4hcGLjoj5+IhIvItyLyP9d2i75fABHpKCLviMhG1//36S39vkXk166f6/Ui8oaIxLS0exaROSKyT0TWe5Q1+h5FZJSIrHO994yIiM9BqGqr/gLCgW1AMhAFrAFSAh1XE91bD+BU1+v2wGYgBXgcuMtVfhfwmOt1iuv+o4F+ru9LeKDv4zju+w7gdeB/ru0Wfb+ue/kn8HPX6yigY0u+byAJ2A60cW2/Dfy0pd0zcDZwKrDeo6zR9wisAE4HBJgHTPY1BnuSgNHAVlXNVNVS4E1gSoBjahKqmqOqq12vC4ENOL9cU3A+VHD9e5nr9RTgTVUtUdXtwFac70/IEJGewMXASx7FLfZ+AUQkDufDZDaAqpaq6iFa+H3jLL/cRkQigFggmxZ2z6q6GMirUdyoexSRHkCcqi5TJ2P8y+OYBlmScD40d3lsZ7nKWhQR6QuMBJYD3VQ1B5xEAnR17dYSvhd/Be4EKj3KWvL9gvMUnAu87Kpme0lE2tKC71tVdwNPAjuBHCBfVRfQgu/ZQ2PvMcn1uma5TyxJOI9fNbWofsEi0g74L/B/qlpQ365eykLmeyEilwD7VHWVr4d4KQuZ+/UQgVMl8Q9VHQkcwamGqEvI37erHn4KTrVKItBWRK6t7xAvZSF1zz6o6x5P6N4tSThZtZfHdk+cx9YWQUQicRLEv1X1XVfxXtcjKK5/97nKQ/17cSbwAxH5Hqfa8FwReY2We79VsoAsVV3u2n4HJ2m05Ps+H9iuqrmqWga8C5xBy77nKo29xyzX65rlPrEkASuBASLST0SigKnAhwGOqUm4ejDMBjao6lMeb30IXOd6fR3wgUf5VBGJFpF+wACcBq+QoKp3q2pPVe2L8/+4SFWvpYXebxVV3QPsEpFBrqLzgAxa9n3vBMaKSKzr5/w8nDa3lnzPVRp1j64qqUIRGev6Xv3E45iGBbr1Phi+gItwev5sA+4NdDxNeF/jcB4r1wLfub4uAuKBhcAW17+dPY651/V92EQjekAE2xcwgWO9m1rD/Z4CpLn+r98HOrX0+wb+CGwE1gOv4vTqaVH3DLyB0+ZShvNEMON47hFIdX2ftgHP4pptw5cvm5bDGGNMnay6yRhjTJ0sSRhjjKmTJQljjDF1siRhjDGmTpYkjDHG1MmShGn1XDOo3uKxnSgi7zTTtfuKyI+a41rGHA9LEsY4M6a6k4SqZqvqlc107b6AJQkTtCxJGAOPAv1F5DsRecL11/16ABH5qYi8LyIfich2EblNRO5wTaT3jYh0du3XX0Q+EZFVIrJERAbXvIiIjHdd4zvX8e1d1z7LVfZrcdbCeEJEVorIWhG5yXXsBBFZLCLviUiGiDwvImGu/V9xramwTkR+3YzfN9MKRAQ6AGOCwF3AUFU9Bdwz5noaijODbgzO9Mu/U9WRIvIXnCkO/oqzCP3NqrpFRMYAfwfOrXGe/wfcqqpLXZMuFruu/f9U9RLXtW/EmdH0NBGJBpaKyALX8aNx1gzYAXwCXI6zpkKSqg51Hd+xKb4hxlSxJGFMwz5XZz2OQhHJBz5yla8Dhrs+8M8A/uOx4Fe0l/MsBZ4SkX8D76pqlpcFwia6zllV3dUBZw6eUpx5eDIBROQNnGlXFgLJIvI34GNgQc0TGnMiLEkY07ASj9eVHtuVOL9DYcChqieRuqjqoyLyMc78Wd+IyPledhPgl6o6v1qhyARqT++sqnpQREYAFwK3AlcDP/PprozxgbVJGAOFOMu7Hhd11ujYLiJXgTP7ruuDuxoR6a+q61T1MZzJ+AZ7ufZ84BeuKd4RkYGuBYTAWWWsn4iEAdcAX4lIFyBMVf8L/B5ninBjmow9SZhWT1UPiMhSV2P1POC54zjNj4F/iMh9QCTOehZrauzzfyJyDlCBM5X3PJynkXIRWQO8AjyN0+NptWta51yOLTW5DKehexiwGHjP9fplV+IAuPs4YjemTjYLrDEhwFXd5G7gNqa5WHWTMcaYOtmThDHGmDrZk4Qxxpg6WZIwxhhTJ0sSxhhj6mRJwhhjTJ0sSRhjjKnT/wfCFs38l2/WRAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure()\n", "plt.plot(bh1.pop_hist)\n", "plt.yscale('log')\n", "plt.xlabel('generations')\n", "plt.ylabel('population size')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "correct-victory", "metadata": { "slideshow": { "slide_type": "subslide" }, "tags": [] }, "source": [ "We see that, even though the population size fluctuates, it roughly increases linearly in log scale, i.e., exponentially with time. (Again, the numbers here are not realistic.) Let us find the slope of this trend, which gives the average growth rate:\n", "\\begin{equation}\n", "\\Lambda_T = \\frac{1}{T} \\, \\log \\frac{N_T}{N_0}\n", "\\end{equation}" ] }, { "cell_type": "code", "execution_count": 6, "id": "specialized-cookie", "metadata": { "slideshow": { "slide_type": "fragment" }, "tags": [] }, "outputs": [], "source": [ "lam_avg = np.log(bh1.pop_size / N0) / T # average growth rate \n", "pop_avg = N0 * np.exp(lam_avg * np.arange(T)) # trend line" ] }, { "cell_type": "markdown", "id": "careful-backup", "metadata": { "slideshow": { "slide_type": "subslide" }, "tags": [] }, "source": [ "Theoretically, in the long term, the asymptotic growth rate is given by:\n", "\\begin{equation}\n", "\\Lambda = \\lim_{T\\to\\infty} \\Lambda_T = \\sum_\\varepsilon p_\\varepsilon \\log \\Big( \\sum_{\\phi} q_\\phi \\, f_{\\phi\\varepsilon} \\Big)\n", "\\end{equation}" ] }, { "cell_type": "code", "execution_count": 7, "id": "spatial-samba", "metadata": { "slideshow": { "slide_type": "fragment" }, "tags": [] }, "outputs": [], "source": [ "lam_theory = np.dot(np.log(np.dot(phe_dist, fit_mat)), env_dist) # asymptotic growth rate\n", "pop_theory = N0 * np.exp(lam_theory * np.arange(T)) # theoretical trend line" ] }, { "cell_type": "markdown", "id": "refined-flour", "metadata": { "slideshow": { "slide_type": "subslide" }, "tags": [] }, "source": [ "Let us overlay these trend lines onto our simulation result." ] }, { "cell_type": "code", "execution_count": 8, "id": "searching-stadium", "metadata": { "slideshow": { "slide_type": "fragment" }, "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEJCAYAAABhbdtlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3yN1x/A8c/JEiKICEUQUTMkQRCborSlFKXaaqtGUd0/VTooWm11GS3V2m1pq/betWLETmLvSEgkZEfG/f7+uHElkXCRm0hy3q9XXu59xnnOE8n95jnje5SIoGmapmlZscrrCmiapmmPLh0kNE3TtGzpIKFpmqZlSwcJTdM0LVs6SGiapmnZ0kFC0zRNy5YOEpqmaVq2dJDQNE3TsmWT1xW4RSnVDXgGKAv8JCLrlVJWwDigBOAvInPvVkaZMmXEzc3N4nXVNE0rSPbv339NRFyy2mfRIKGUmgV0BsJEpG667Z2ASYA18JuIfCUiS4GlSikn4FtgPdAVqAhEAsH3up6bmxv+/v45fyOapmkFmFLqQnb7LN3cNAfolKky1sBPwFNAHaCPUqpOukM+SdsPUBPwE5H3gSEWrqumaZqWiUWDhIhsw/gUkF5j4LSInBWRJGAh0FUZfQ2sEZEDaccGA9fTXqdmdQ2l1CCllL9Syj88PNwCd6FpmlZ45UXHdUXgUrr3wWnb3gLaAz2VUoPT9i0GOiqlpgDbsipMRGaIiI+I+Li4ZNmkpmmapj2gvOi4VllsExGZDEzOtDEe6P8wF0tOTiY4OJjExMSHKUbT7mBvb4+rqyu2trZ5XRVNs5i8CBLBQKV0712BEItdLDgYR0dH3NzcUCqr+KRp909EiIiIIDg4mKpVq+Z1dTTNYvKiuWkfUF0pVVUpZQe8ACy31MUSExNxdnbWAULLUUopnJ2d9ROqVuBZNEgopRYAfkBNpVSwUqq/iKQAw4B1wDHgbxEJtHA9LFm8VkjpnyutMLBoc5OI9Mlm+2pgtSWvrWmaVlisOBzCjfgk+jZ1y/GydVoOLddt3bqVXbt2md6/9tprLFq0KMfK//HHH4mPj8+x8jTtUbdo72m27d1nkbJ1kMinUlOznDaSa1JSUh743MxB4n6JCAaDIdv9Okhohcm5A5v55NJgvogfC6kP/nuZHR0kckG3bt1o2LAhHh4ezJgxA4Bp06bx4Ycfmo6ZM2cOb731FgC///47jRs3xtvbmzfeeMMUEIoXL85nn31GkyZN8PPzY+zYsTRq1Ii6desyaNAgRASAffv24enpSdOmTRk+fDh16xozoqSmpjJ8+HAaNWqEp6cnv/zyS5b1HTduHLVq1aJDhw706dOHb7/9FoA2bdowatQoWrduzaRJk7hw4QLt2rXD09OTdu3acfHiRVJTU3F3d0dEuHHjBlZWVmzbZpzi0rJlS06fPs306dP54Ycf8Pb2Zvv27QBs27aNZs2a4e7unuVTxfnz56lduzZDhw6lQYMGXLp0iSFDhuDj44OHhwejR48GYPLkyYSEhNC2bVvatm0LwPr162natCkNGjTg+eefJzY29iH+NzXtEXEzhuuL3qHKsu4UVTfZX+tDsLZAD4KIFJivhg0bSmZBQUG336weITLr6Zz9Wj3ijmtmFhERISIi8fHx4uHhIdeuXZOwsDCpVq2a6ZhOnTrJ9u3bJSgoSDp37ixJSUkiIjJkyBCZO3euiDECyF9//XVHuSIiL7/8sixfvlxERDw8PGTnzp0iIjJixAjx8PAQEZFffvlFxo0bJyIiiYmJ0rBhQzl79myGuu7bt0+8vLwkPj5eoqOj5fHHH5eJEyeKiEjr1q1lyJAhpmM7d+4sc+bMERGRmTNnSteuXUVEpGPHjhIQECArVqwQHx8fGT9+vCQmJoqbm5uIiIwePdpUpojIq6++Kj179pTU1FQJDAzM8H255dy5c6KUEj8/vzvuPyUlRVq3bi2HDx8WEZEqVapIeHi4iIiEh4dLy5YtJTY2VkREvvrqK/n888/vKP9BZfj50rRccuPwKpHv6ohhdEmZ9fHzUvvDv+W/gAsPXB7GBKpZfq7qJ4lcMHnyZLy8vPD19eXSpUucOnUKFxcX3N3d2b17NxEREZw4cYLmzZuzadMm9u/fT6NGjfD29mbTpk2cPXsWAGtra3r06GEqd8uWLTRp0oR69eqxefNmAgMDuXHjBjExMTRr1gyAF1980XT8+vXrmTdvHt7e3jRp0oSIiAhOnTqVoa47duyga9euFC1aFEdHR7p06ZJhf+/evU2v/fz8TOX37duXHTt2AMYnhm3btrFt2zZGjhzJjh072LdvH40aNcr2e9StWzesrKyoU6cOV69ezfKYKlWq4Ovra3r/999/06BBA+rXr09gYCBBQUF3nLN7926CgoJo3rw53t7ezJ07lwsXss1lpmmPtrgIwue+QsnFfbiSaE2Pm6P5JKw9tlt+ZOXcqRa55COTKjxXPPVVrl9y69atbNy4ET8/P4oVK0abNm1MY+t79+7N33//Ta1atXjuuedQSiEivPrqq0yYMOGOsuzt7bG2tgaM8z+GDh2Kv78/lSpVYsyYMSQmJpqanLIiIkyZMoWOHTve9Zi7cXBwyHbfrSGhLVu2ZPr06YSEhDB27FgmTpzI1q1badWqVbbnFilS5J51SH/tc+fO8e2337Jv3z6cnJx47bXXspyzICJ06NCBBQsW3PW+NO2RJgJHF8HaEZROiGJSSnemxj3Ntb0rid63lOFffm5qrs5p+knCwqKionBycqJYsWIcP36c3bt3m/Z1796dpUuXsmDBAtNf6O3atWPRokWEhYUBEBkZmeVfvrc+EMuUKUNsbKypHd/JyQlHR0fTdRYuXGg6p2PHjkybNo3k5GQATp48SVxcXIZyW7RowYoVK0hMTCQ2NpZVq1Zle2/NmjUzlf/HH3/QokULAJo0acKuXbuwsrLC3t4eb29vfvnlF1q2bAmAo6MjMTEx5n4LsxQdHY2DgwMlS5bk6tWrrFmzxrQvffm+vr7s3LmT06dPAxAfH8/Jkycf6tqalquiguHPXrB4ADi58X3V3/jqYj0uzBlOYnAg85dv4t133zX9AZnTCteTRB7o1KkT06dPx9PTk5o1a2ZoLnFycqJOnToEBQXRuHFjAOrUqcP48eN58sknMRgM2Nra8tNPP1GlSpUM5ZYqVYqBAwdSr1493NzcMjTlzJw5k4EDB+Lg4ECbNm0oWbIkAAMGDOD8+fM0aNAAEcHFxYWlS5dmKLdRo0Y8++yzeHl5UaVKFXx8fEznZzZ58mRef/11Jk6ciIuLC7NnzwaMTwWVKlUy3WvLli1ZsGAB9erVA6BLly707NmTZcuWMWXKlAf6vnp5eVG/fn08PDxwd3enefPmpn2DBg3iqaeeonz58mzZsoU5c+bQp08fbt68CcD48eOpUaPGA11X03KNwQD+M2HjGBADdJxATJ0XmdplAOEHNuP0xACK1W7Fi+19LFoNda/mhfzEx8dHMi86dOzYMWrXrp1HNcobsbGxFC9eHICvvvqK0NBQJk2adN/nx8fH06pVK2bMmEGDBg0sVd18rTD+fGm5IPwkLH8LLu0G97bQ5UdW7DjK0DffJMqpFu9/MpaODapzMyWVNjXLPvTllFL7RSTLaKOfJAqgVatWMWHCBFJSUqhSpQpz5sy5r/MHDRpEUFAQiYmJvPrqqzpAaFpuSU2GnT/Cf9+AbTHoNo0r5dry9hvvcPDgQT7/9ifGHrDCp2ZlmlZzzpUq6SBRAPXu3TvDKKT79eeff+ZgbTRNM8vlA8anh6sBUKcbhk5fM/OvFXz8sRcDBgxg7ty5rAi4BgeO4FEh6yZgS9BBQtM0LS8lxcOWL2D3z1C8HLzwJyfU4wzq0pvExEQ2btyIp6cnAIEhUTjYWVOldLFcq54e3aRpmpZXzm6FaU3Bbyo0eIWkgdsZ988BmjdvTo8ePdi1a5cpQJwNj2Wu3wWqujhgZZV7GYj1k4SmaVpuS7gO6z+Bg79DaXd4bRV+odYMbNaWqlWrcuDAASpXrpzhlPGrjgFwLSYpV6uqg4SmaVpuCloGq4dD3DVo/i7RDYYy8rOxLFmyhB9//JHnn38+y7VKUgzGkahf9/TM1erq5iYLu3HjBj///LPp/datW+ncuXOu1iFz1tXp06czb968Byrr1tDa3KbTi2v5XswVWPgS/P2Kse9h4GaWxTfAw9uHpKQkAgMD6dWrV7aLWZ27FssznuVpXcMlV6utg4SFZQ4SlnK31N2ZP2AHDx7MK6+8YvE6ZabTi2uFkgjsnwtTG8PpjdB+DCHP/E7Pt8fz4YcfMn/+fH799VecnJwynHYmPJa+M/fgfz6SwJAoLkUmUKd8iVyvvg4SFvbRRx9x5swZvL29GT58OGCcrNazZ09q1arFSy+9ZMpVtH//flq3bk3Dhg3p2LEjoaGhABw6dAhfX188PT157rnnuH79OnBn6u7w8HB69OhBo0aNaNSoETt37uT8+fN3pOYeM2aMKf336dOnad++PV5eXjRo0IAzZ84QGxtLu3btaNCgAfXq1WPZsmX3vE+dXlzTshBxBuZ2gRVvw2P1MLyxg18Ci+HVwIfatWtz+PBh2rRpk+WpUzadYvupa/Sc7sczk43JMz0q5H6QyPP03jn5dc9U4Xng3LlzplTdIiJbtmyREiVKyKVLlyQ1NVV8fX1l+/btkpSUJE2bNpWwsDAREVm4cKH069dPRETq1asnW7duFRGRTz/9VN555x0RuTN1d58+fWT79u0iInLhwgWpVauWiNyZmjv9+8aNG8vixYtFRCQhIUHi4uIkOTlZoqKiRMSYartatWpiMBhERMTBweGOeyzM6cXz+udLe0SlJIvsmCQyrpzIl64i+2ZJUECAtGjRQnx9feXo0aP3LGLo7/ulyoiVGb7CohMtUl3ukiq80HVcW2LxernP1CaNGzfG1dUVAG9vb86fP0+pUqUICAigQ4cOgHGBoPLlyxMVFcWNGzdo3bo1AK+++irPP/+8qaz0k+Y2btyYIV12dHT0XRPpxcTEcPnyZZ577jnAmGUWIDk5mVGjRrFt2zasrKy4fPkyV69e5bHHHsuynPTpxYF7phdfvHgxYEwvfmvhpVvpxc+dO8fIkSP59ddfad26tUXSi8+YMYOUlBRCQ0MJCgoyDTG8JX16cYCkpCSaNm2abT00LYMrR2HZMAg9BDWf5mb7L/nq53lMmTKcMWPGMGTIEFMyvv0XIqlc2gEXR2MW5IjYmySlGihfsijHr0RnKNbFsYjpuNxU6ILE/X6gW0L6tNjW1takpKQgInh4eODn55fh2KioqLuWlT59tsFgwM/Pz/RhfS/ZfS/++OMPwsPD2b9/P7a2tri5uWWZhvte5WRVx8x0enGtwEhOhG3fwM5JUNQJes5m542yDGz9DNWrV+fgwYNUqlTJdHhSioEe0/yo4lyMLR+04diVaFOzUtDYjpy9Fsfb7arTv3lV6o9bnzdNTTxifRJKqW5KqV+VUsuUUk9mty0/MTctds2aNQkPDzcFieTkZAIDAylZsiROTk6mdvj58+ebnioye/LJJ5k69fbCI4cOHbprHUqUKIGrq6spE+zNmzeJj48nKiqKsmXLYmtry5YtW+65SI9OL64Vehd2wfQWsP07qNeLqJc3MuSnDfTq3ZuxY8eydOnSDAECYPGBYOOpEfFMWHPMFCAA3ll4CBGoW6EEJYvZ8tYT1enrmzETdG6xeJBQSs1SSoUppQIybe+klDqhlDqtlPoIQESWishA4DWgd3bb8hNnZ2eaN29O3bp1TR3XWbGzs2PRokWMGDECLy8vvL29TaN55s6dy/Dhw/H09OTQoUN89tlnWZYxefJk/P398fT0pE6dOkyfPh0wNv8sWbIkQ6fvLfPnz2fy5Ml4enrSrFkzrly5wksvvYS/vz8+Pj788ccf1KpV6673mD69ePfu3e+ZXnz27Nl4enoyf/58U3barNKLx8TEZEgvnt09mCt9evHXX389y/Tibdu2xcXFxZRe3NPTE19fX44fP/5A19QKuMRoWPk+zH4KUm/Cy4tZwpN4NGqBwWAgMDCQnj173tHMLSJ8tPio6f2v289l2L8hyNh86lHR+Hv0XocatKtdzsI3kzWLpwpXSrUCYoF5IlI3bZs1cBLoAAQD+4A+IhKUtv874A8ROZCunDu2ZaZTheedwppeXP98FWIn1sKq9yE6BHyHcLlGP4Z9MIJjx44xY8YMU1NpUoqBVINgb2vF2WtxjPz3KC80rsT7fx++5yXOTXjaIv2omeVpqnAR2aaUcsu0uTFwWkTOAiilFgJdlVLHgK+ANbeCgTJ+hzJs0x49Or24VmjEhsPaERDwL7jUxvD6HH5ZtZ/P+jZn6NChLFy4MEN/2ZgVgRy4cJ3jV243l+49HwnANz08+WL1MaISjKtFTn+5AYN/N37MLXuzea4EiHvJq47risCldO+DgSbAW0B7oKRS6nERmZ7NNhOl1CBgEHBHrhMt9+j04lqBJwJH/oK1I+FmDLQZRVDpTgzsMxQwTvj08PC447Q9ZyM4Ex53x3aAp+o9xsZjV1kfdJXf+zehRfUyjHq6Fl+uPk7NxxwtejvmyquO66zCo4jIZBFpKCKDbwWDrLZlOmmGiPiIiI+LS9bT1R+FEU1awaN/rgqRGxfhj56w5A1wrsbNfpsYvSWe1u3a8/LLL7N9+/YsA0R8Ugpnr2UdIIZ3rImjvS1P1CpLUVtr6lY0jl4a2NKdcxOext7WMmtW36+8epIIBtJ39bsCIZa4kL29PRERETg7Oz8Sj25awSAiREREmOaWaAWUIRX2/YZs/JwUgwGrjl+xK6k2Azv0pnbt2hw6dIiKFStme/qx0BjS/y3xUpPKbDsVzqXIBDp7lgegd6NKdPR4jFLF7ADLzOV6GHkVJPYB1ZVSVYHLwAvAi5a4kKurK8HBwYSHh1uieK0Qs7e3N02K1AqgsOPGleKC9xJcuhk9z3bB4evtnPT/gilTppgmoaZ3LfYm7/11iDfbPo6zgx1BocYJcd6VSlHRqShfPFePVINwNjyWKs7GOTxKKZwc7HL11u6HxYOEUmoB0AYoo5QKBkaLyEyl1DBgHWANzBKRQEtc39bWlqpVq1qiaE3TCpDAkCjqlC+BSk2GHT/A9m/BzoHUrtNo+kUgV9Z9ileL9qb5S1nZefoa208Zv9JbMrSZ6QnB2kpRvdyj0d9gjtwY3dQnm+2rgdWWvr6madrdRMYl8cWqY/x7IJifW6fy9LkJEBZEap3uhHq9y/NDRhB2JIgyXT/Cvm5DsM166VARITAkOst9j1oT0v0odGk5NE3T0hu1+Cj/BZ7nU5t/6LRnLVdxYkTi+1zbnMSeN9tTokFnyr82GGVjy+UbCXiNXU+d8iWY3a8R5UoY+6R2nb7GoPn7qVDqzj6qH3t75/Yt5ahHKi2HpmlabnMJ28V6uxH0t1nDCttONL80jAXzF7Brw0pmLVpNav2eKBtb5vdvbDonKDSaJl9uIuRGAgA7Tl8j9mYKJ6/emVK+W/3sO7bzA/0koWla4RQfCes+ZlzMn5yhPN3jRrJp5yFiD42lVKu+VG7ahWTH8kAYAE2qOt9RxKojoQxs5Z6hmWlMlzrUr+zE/gvXKW6f/z9i8/8daJqm3Q8RCFoKq4eTGh/JtJSu/Icve/75jpQij1G+3xSeb+XJkoOXWXbIODL/6JgnsbOx4sT4Tny8JIAVh0O4mWLgWNropfRBom7FknhVKoVXpVJ5cns5TQcJTdMKjY17DuG66xNqRW3nkn1N+t0YzJ7NG7C/Op5fpv3E+3uMQ1HfbPs4Sw5e5r+T4Xi5lsTR3haAIjbWfPu8F191r8fg3w+w+OBlvCqV4lrsTbwqleJiRBy182CJUUvSQULTtILPYIADc2m8eiS2pDI+uQ9TDpYgYvMEitVoyuGjAVRwKc1z3QSlFCJC+ZL2hEYlUqfCncNdbaytaF3ThY3HrjJ6uXH0/sinatHYrTRWVvl3JFNWdJDQNK1gizgDy9+GCzsIMNTh3YjuBKz/h5Soq5TpOpKp771ABZfSwO2hqkopfN2dWXLwcraL/TRyc8rwvk6FEgUuQIAe3aRpWkGVmgzbv0d+bkrcxYP872Z/OvvVY/+c8dhVqEmFfpOwd62dbRBoWs3YUV2vYtYT52o9VoJ32lUHwMZKUSKtSaqg0U8SmqYVPCGHjCk1rhxhbWoj/hfSjtNr5qBs7Fi9YQtPNqvPqCVHWXLwMo+XLZ5lEd3rV6SsYxE8XbMOEmBcDKhHA1cqlTZvyeD8SAcJTdMKjuQE2PoV7JoCDmUYEPcm/+w4TeyRrynVsi8duvWhY3PjWif/e7ImLzSqjK111g0qNtZWtKlZ9p6XrOyc9QzsgkIHCU3TCobzO4x9D5FnoP7LbCnSkXlfDsG2rDvl+03BpnhpPu58O523c/EiOBcvcpcCNdBBQtO0/C4xCjZ8BvvngJMbkc/O539Tl7J+wzBKPdGfd1/vgwI+eLImdja6G/Z+6e+Ypmn51/FV8FMTODAP8X2Thc7v4/HMGxQvXpzZK7ZR7PEmtKruwsina+sA8YD0k4SmaflPbBisHm6cOV2uLhdafMuQz3/m0qUVLFmyBF9fX2ZsOwOQ7eglzTw6tGqalqOSUw3E3UyxTOEicPAPmNoITqwmtfUofkzoTsNnXqNFixYcOHAAX19fwJgqo0JJ+0d6QZ/8QD9JaJqWY6ITk/Ecsx6AcxOeztl1FK6fhxXvwtktUMmXw+5vMuDDcRQvXpxdu3ZRo0YN06HfbzjJskMhtK9979FJ2t3pIKFp2kMTEWbuOEdkXJJp25nwuGznINwXQyrsmQ6bx4OyIuGJL/l81UVmjRzAV199Rb9+/TIEo7CYRCZvOgWQZUoN7f7oIKFp2kMLjUpk/Kpj2KRLS7H7bES2QSI0KoH9F67T2bPC3Qu+GmScFHfZH6o/yUaH5xj8+mf4+Phw9OhRypUrB4DBIHyyLIAXGlUiIl2gqlPAku3lBR0kNE17IJFxSSSlGHispL0pVXaKQahetjgRcUl8sjQAR3sbXByL0KxaGdN50YnJNJ2wGTAGkvHd6gGQmJzKH3su0te3CnYkw/bvYPv3YF+CiCd+4INZ29iyZSQ///wzzzzzTIa6nIuI4889F0lNFSo63Z79rDutH54OEpqm3bdZO84xdmUQAOe/eoagdOspeFQowaFLN4iMS+KdhYcACPy8Iw5FjB83x9Id+/vui7zTrgYujkVYF3iFcSuDqJNyjKYBY+DaCaTu8yyIacQHvT+md+/eBAYGUrz4nU8nt4LUX/6XTNu2f9iWSqUL9mzo3KBHN2madl9SDWIKEABh0YkEhkRhl5beok6FErzQuHKGc/ZfuG56nX6BHoA95yIAOHXpCqNt5uK79UVIisO/4Xe4fh7E+O9+osXQr/nf6C+zDBBAhiB1iw4QOUMHCU3TsmUwCJM2niL4ejxg7KDO/IG8+1wkQaHRdPAox7hudenlU4k3WrnTv0XV28ecjTC9DgqNpkxxO/4d0uz2vlMbee3QC7xqvZ4l1k/xXfJLtHppOPHONYh7ahz74pz5cvWxbOsZGBKVYbLckqHNcuT+tUeouUkp5Q58DJQUkZ5p2yoDU4FrwEkR+SoPq6hpj4wb8UlExiXh7pIDo4fu4viVGH7YeJIfNp7kn8FNWXLwMn/uuQjAHwOaMHj+ftYFXiH4egIvNalCX98qpnOf9arAzB3ngIxBIjAkmjoVStKwihOdH7ejxdGP4dB/xFCRF4IHsH31MprUvEKvsXPYeuX2R9SVqMQs63grcHX1qsCz3hVo5FYae1trS3w7CiWLPkkopWYppcKUUgGZtndSSp1QSp1WSn0EICJnRaR/piJqAKtE5HWgjiXrqmn5RVh0It5jN/DEd/+RahCLXuv4ldtPDc9P9zMFCIDGVUvTqGppVh0JBYzNTOl5VSrFpg9aM7RNNQ5cvMGgef7cTEnl1NUY6jzmCEcXMfHqQNql7GBiQhcarKnBrkWzcWzQmZFT/uAKGRf1ORYaQ3Kq4Y46hsXcJCIuCY8KJWhZ3UUHiBxm6eamOUCn9BuUUtbAT8BTGD/4+yilsgsAB4EXlFKbgS0WrKem5Rvrgq6aXmfVFp+T7la+rbUVPulWZ8tqJFE1l+KmxXvWB10lKCQaF8M1Xrv4EfzbH0PJyjQ98TKjZmwiKeYG/2zYQfF67Rnyx0GOX4lhcOtq9GjgiruLAwnJqRwJjrrjGoEhxm16ToRlWDRIiMg2IDLT5sbA6bQnhyRgIdA1myL6AaNF5AngmawOUEoNUkr5K6X8w8PDc6rqmvbICo6MN71O34xzP1YdCc0w8S07gSHR1HrMEa9KpTJs/6aHJwCN3UqbtpXJJu12wyrGQKIwcG7NZNYX+ZCy13Zzrcko3thdmUNr/qJ0+8G4PDucVp6PZzi3XsWSfNfLi3/eaApkfb+3Alnt8o73vB/t/uVFx3VF4FK698FARaWUs1JqOlBfKTUybd9a4O207eezKkxEZoiIj4j4uLi4WLLemvZICAqNpm7FEriXcXigIBFyI4E3/zxAg3EbuBptbOc3GITYTPmWRISg0GjqV3Zi8ZBmHB/XiXHd6rJjRFt6NaoEgI9baVa93YKzXz6d7fWK2dmwZ1BlFtqNp3vo9xxMdWe+8wjq9vuBcuXKceX8SUrWaIybczEc7W15qcntkVG3nk6cixehZjnHLO83MCTadK6W8/Ki4zqrZC4iIhHA4EwbA4CeuVIrTcsHRITAkGg61C6HlZVi5eEQUg2CtZV5OZIMBuH7DSdN75t8uYm/32jKrB3nWBt4haCxHbG1tmLK5tM8VsKeqIRkPCqUwNpKYW1lnaFj+haPuzXzpCbDzkmU++8b7JU1A8N68fsaf2qVms/KlSvx8fEBoItXBVwcjU8iY571IDzmJuuDrlI53TBWX/fS/O0fTHKqwbSa3JHgG6wJuMJTdR8z6/61+5cXQSIYqJTuvSsQkgf10LR8Z9p/Z4iMS6JOhRKUKmbLgr0XCQqJpt5d1mFO74N/DrPk4OUM23r94md6feDCDbaeCOO3tFFJ8BCzli8fMJtsXl4AACAASURBVK4Ud/UoKTW78MRCGw5t+Js2z/dn7W9fY2t7+y//73p5mV7bWlvxS9+GGASs0gU/X3dn5vpd4EhwlKkJ6+Xf9gDotSIsKC++s/uA6kqpqkopO+AFYHke1EPT8p1v1p4AjCOJfN2NHcL30+SUOUBktvtsBKuPhpreWymo9dh9BomkeFj/CfzWDuLC2V/3cxpPDMIx5hJ/r9rEprnfZwgQWVFK3fF01LhqaVMdAc6GxxKdaGwie62Z2/3VUTObpYfALgD8gJpKqWClVH8RSQGGAeuAY8DfIhJoyXpoWn4nIvx38vbAjDrlS1CuhP199UsYDMKtz909o9pxbGwnqjhnnJU8dctpQtLNR3B3KU5Ru/sYUnr2P5jWDHZNIa72C3xwqQNPvzmBd999ly2bNtKjrY/5ZWWSuV/i02XGkfXVXByoX9npbqdqD8HSo5v6iEh5EbEVEVcRmZm2fbWI1BCRaiLyhSXroGn5XUTsTVYcCeXVWXsBmP5yQ1MepCbuzuw9F3nHfAkR4eetp7kUGc+x0Gjaf/8fBy5exyAwoXs9ypWwp6idtamTuItXBVrXuD3w443W7tjZWJnf1JRwA5YNg3nPglKsq/opdUesJSwyioCAAF555ZUcWVvC1700/uevE3w9np2njcHik2f0FCpLemRmXGuadqfIuCQajt+YYVv6D25f99KmfonElFS+XnOcr3rUo+/MvYRGJfLN2hP08nHldFgsP2+9cznPAS3caeRWmvqVnVgbEGp6WvGsWIoZfRtSxdnh3pUMWg6r/wdx1wirM4D3llzGb8/PTJ8+nY4dO+bAd+G2W/0SQ/84YNrWtpZeWMiS7hkklFI1gGlAORGpq5TyBJ4VkfEWr52mFVLHQqN5atJ2yjreOffANV0q7Fv9El2m7qCsYxHCYm4y5PcDhKZrMvrbPxiAzcfDAKhR7vZ8AisrZWqq6VS3POO61eXTpQE0rOLEYyXt717JmKvG4HBsOVKuLvNsXubDt37ilVde4ehvM3FwMCPA3Kdb/RJHLxsn0H2frsNbswxzniR+BYYDvwCIyBGl1J+ADhKaZiHP/bwTMKacSO/nlxpkaLYpV+L2B/mtY0+Hx9617LulrejrWyXLYa4ZiMDB32H9x5CcyJlaQ3nj171EXl/KmjVraNCgwd3PfwjOxYtQo1xxTl6NpcXjZejewNVi19KMzOmTKCYiezNts9Aq55qmARgypShq8XgZ9oxqx9P1yt9xbO1Mq69JWvdEWccitKlp7GewtzX+qr/9RMYZzfct8hzM6wrLh5FcujZfJ79Kk2G/8tTTT7N3716LBohbbj09Zc4VpVmGOU8S15RS1QABUEr1BELvfoqmaQ8qOdUACp6oVZbNx8OY9II3Xb0rZnv8goFNaPzlJpJSDDRyc2Lf+ev0b1GVj5+uzY7T19h6IpxpLzXEo2IJyjreowkpO6kpsGcabP4CrGzYV/UtBv6wknLljrBv3z6qVq167zJyiK+7M/P8LuhV53KJOUHiTWAGUEspdRk4B7xs0VppWiEVHnOTDxcdJinFQBev8sx6rdE9zylVzI4jo58kKdXAoYs3eGXWXupVLImVlaJl9TL8ObAJTd2dH3x00ZUAWD4MQg4SW7k9n/qXZsHk3/juu+948cUXc2TU0v1oV7ssI5+qxZN19Czr3HDPICEiZ4H2SikHwEpEYixfLU0rnL7fcIItJ4wjjO6a7iITe1tr7G2tafF4GSa94M1TdY3NUkqpDOtL35fkRNg2EXb+CPalWF12GEPH/kHr1q0JCAigTJkHLPchFbGx5o3W1fLk2oWROaObUoGJwEgRudXkdEBELN/4qGmFSExiMgv23s596V7m/kcHWVmpuzZNme2CH6x4G66d5GqVrry7Oo69Bxby66+/0qFDh4cvX8s3zOm4Dkw7br1S6lZe4Nx9vtS0QiD9gj4ANtZ5kDUnMRpWfQCzOyFJCcy2H4jnR+uo7P44R48e1QGiEDKnTyJFRD5USvUCtiulXiGtE1vTtJxzLdY4hLW3TyX+17Fm7lfg5DpY+R5Eh3DKtTdv/H6C6NitrF27lvr16+d+fbRHgjl/qigAEfkb6AXMBtwtWSlNK4yOhcZQt2IJvu7paUqbnSvirsG/A+DPXiTbOPBlwos0/WgxXbp2Y/fu3TpAFHLmPEkMuPVCRAKVUi2AbparkqYVPrcW+GlfOxdTTIjA0X9gzQi4GcOesi8xcNp2Krqew9/fHzc3t9yri/bIyjZIKKWeEJHNQBWlVOYpmHef0qlp2n25Ep1IZFzSfY1oeig3Lhmblk5vIMbZm49PVeKfaUv4/vvveeGFF3J9WKv26Lrbk0RrYDPQJYt9Aiy2SI00rQASkQwfvOeuxfHJ0qN81Kk2pYrZcvKqcWS5xWcRGwyw7zfY9DmIgZXFX+LNb1fxxBPuBAQE4OzsbNnra/lOtkFCREan/dsv96qjaQXPiEVH+Mv/EkVsrLiZYmDhIF9emLEbMCbmA3ivfQ3gzhQbOSr8BCx/Cy7t4YpzM97ZbM2BgA3Mnj2bJ554wnLX1fK1e3ZcK6XeUUqVUEa/KaUOKKWezI3KaVp+JyL85W+c+3AzxZiQ6Y35++84btmhy7g5F6N4EQtk709Jgv++gektkPAT/Gbojee4/VSr7cmRI0d0gNDuypyfyNdFZJJSqiNQFuiHcYTTeovWTNPyqZNXYyhVzJZSRe2o8cmaO/ZHJSQDYGOlSElbLOjstTieySJ530ML3m98eggL5IRTO95YdJX4m4fYuHEjnp6eOX89rcAxewgs8DQwW0QOoyfTaQVcVEIyIvc/HUhEePKHbTT+YhOL9gdn2Pd8w4xprb99PuNaCDnaH5EUB2tHwcz2JMVEMv56V5qP20b353vj5+enA4RmNnOeJPYrpdYDVYGRSilHwHCPczQt3wqPuUmLrzcz8XkvnvWqkO1xBoOw8mgoLR4vQ2kHOwB2n4007R+15Kjp9dExT+Job0uqQVh88LIp6Z5DERvWBISy+MDlnAsSZ7bAinfgxgV2F3+KgXOOUsUtkgMHDlC5cuWcuYZWaJgTJPoD3sBZEYlXSjljbHLStAIlKCSag5euU87RnpspBtYGhJKQlELPhpWwtsr48BwZl0SDcRtM7/eMakdQaDT9Zu/LcJyzgx3+n7Q3jWya+LwXX3avZ1r4p0OdcihgQ+BVvFxLPdwNJFyHdZ/Aod+JLubGqLPtWbxxGz/88AO9evXSw1q1B6Ie5JH6UeXj4yP+/v55XQ0tn6o3Zh0xiXeup1WqmC3bP2zLzB3nGNjSnf9OhmdYYzn9cTfikzO8blm9DPP7N7nntQ0GwcrqAT/ERSBoGaweDvERLFcdGfbbDp58siPffPMNpUuXvncZWqGmlNovIj5Z7bPAUIoHo5RyBz4GSopIz3TbHYBtwGgRWZlX9dMKvqwCBMCN+GTeWXiIzcfDMBiEqVtOZ3scwJQ+9an5mCNP/rDN7CakBw4Q0aHGdaaPryS0WG3e2ufGkZP7mTdvPm3atHmwMjUtHYummVRKzVJKhSmlAjJt76SUOqGUOq2U+giM61aISP8sihkB/G3JemqF23y/84xbGXTXYzYfDwNg8ubTGNI9fB/+LONo8BebVKaLVwWqly3OyKdq8VLje6wX/aBEYP8c+KkJhpMbmBHfEa9vz1LLuwlHjhzRAULLMWY9SSilrIFy6Y8XkYvZn2EyB5gKzMtU1k9AByAY2KeUWi4id/yWKqXaA0HAA665qGl35z12vekJAGBom2r8vPUMw9o+ztQtp6lYqiiXbyRkee7bTzxOyWK2BI3tyIC5/uw6E2FaUlMpZbmFcSLOGDumz2/neJH6DFoTT5JcZNOmTdSrV88y19QKLXMWHXoLGA1c5faoJgHuOYZORLYppdwybW4MnE5b8Q6l1EKgK8ZgkFlbwAGoAyQopVaLiB5ZpT2UxORUbK2tiE5IzhAgAHo3qsSQNtVwsLOhf4uqFLWz5uXf9uB/4XqG46a/3JBOdY3LZxazs+GFxpXZdSbi4Tuf7yY1BfymwtYJJIktX115gin/7mD06NEMGTIEa2try11bK7TMeZJ4B6gpIhE5dM2KwKV074OBJmmjpr4A6iulRorIBBH5GEAp9RpwLasAoZQaBAwC9PA+7Z5EhFqfrsWpmC3XMwUIR3sbKpcuZhoF5JQ2rPWXvg3ZdSaCNjVdqDfGOIe0Tqb0GZ3rlaeSU1HqVrRQgr7QI8Z1pkMPs1M1YdDCC1SrnsqBAweoVKmSZa6paZgXJC4BUTl4zax66CQtCA3O6gQRmZNdYSIyA5gBxtFNOVFBreAKvm5sOkofIGqUK87vA5rgUrxIlsNEnYsXoUvafInj4zpx+NINKjsXy3CMlZWifmWnnK9wciL89zXsnESUKsnI481Ytv0IkyZNokePHnpYq2Zx5gSJs8BWpdQq4OatjSLy/QNeMxhI/6ePKxDygGVp2n0JDMn4946ve2kWDmpq9vn2ttY0cc+lTKnndxrXmY44zZLEprw1/whPP1OJgIA/cXKyQEDStCyYEyQupn3ZpX09rH1AdaVUVeAy8ALwYg6Uq2n3FBgSDUCl0kWZ/EJ9y/z1/7ASo2HjaPCfRYgqzzA/TwLPn+PPBQtp1apVXtdOK2TuGSRE5HOAtHQcIiJmLziklFoAtAHKKKWCMc51mKmUGgasA6yBWSIS+CCV17T7FRQSTY1yxVn/Xuu8rkrWTqyBle9jiA5lRmRTPl3oz5Ahvfhz2Sjs7fUgPy33mTO6qS4wHyid9v4a8Io5H+wi0ieb7auB1fdXVU17OJuOXWXT8TA6e1og2+rDig2HNR9C4GKCUt0YtM4Ng000W7f+h4eHR17XTivEzGlumgG8LyJbAJRSbYBfgWYWrJem5ZgrUYl8sfoYKw4bu76K2DxCQ0VF4PBCWDeSm/GxTLjQmJ9WHeLzzz9n8ODBWFlZdL6rpt2TOUHC4VaAABCRrWmpMjQtX5ix7awpQAC80656HtYmnesXjOtMn9nEjsQaDFxsTc06Dhw8eBBXV9d7n69pucCs0U1KqU8xNjkBvAycs1yVNC1nxd68PdzVo0KJO4av5jpDKuydAZvGcSNR+CjQm5V7zzB58mS6d++et3XTtEzMeZZ9HXABFgNL0l7rVOFavnH8Sozp9Tc983ixnbBjMKsjsmYE/4ZVwWOmQNnaBAQE6AChPZLMGd10HXg7F+qiaTkuJdXA8SsxvNK0Ch89VYtidnmU+DglCXZ8D9u+JTixKMN21+RE6A3++nsRLVq0yJs6aZoZsv2NUUr9KCLvKqVWYMzVlIGIPGvRmmlaDjgTHkdSioEGlZ3yLkBc2gfL38IQFsS0y3UZs/Qkw97qzF8ffUSRIkXypk6aZqa7/dbc6oP4NjcqommWcGuGdY6uH22um7GweTzsmU5gvDMDN7piVcyKbdu3U7t27dyvj6Y9gGyDhIjsT3vpLSKT0u9TSr0D/GfJimlaTggKiaaIjRXuZXJ5QN7pjbDiPRIjLvDlGQ+mbTzFuHHjGTRokB7WquUr5vy0vprFttdyuB6aZhGBIdHUeswRG+tc+mCOj4TFb8DvPdh2IRnvhU4EJlXg0KHDet6Dli/drU+iD8acSlWVUsvT7XIEciptuKZZjIgQFBrN0/VyYYa1CAQuhtUfcv16JB8ers6aI6FMnfoT3bp1s/z1Nc1C7tYnsQsIBcoA36XbHgMcsWSlNC0nXL6RQFRCsmm1OIuJugyrPkBOrGbRlUq8s8yW53o0IXDBl5QsaaH1JTQtl9ytT+ICcAEwP4+ypj1CgtIyvlqs09pggP2zYcNoLl2/yVC/ypy9LixavJRmzXTWGq1guGcDqVLKVym1TykVq5RKUkqlKqWic6NymvYwAkOisVJQ+zELBIlrp2FuZ1JXvMeUoFI0mG2gccfeHDx4UAcIrUAxZ+D4VIxrPvwD+ACvAI9bslKa9rDCYhKZtOkU7mUcKGqXgwn9UpNh12TY+jVHI6wYuKksdk6l2b5jEbVq1cq562jaI8Ks2UUiclopZS0iqcBspdQuC9dL0x6YwSA0/mITAGExN+9x9H0IOQjL3yIx+AjjAqvw645QvvhyBP3799ejlrQCy5wgEa+UsgMOKaW+wdiZrbPAao+si5HxptevNqvy8AUmxcPWCeD3E1tCi/HGWke8GtXi8JH1lC//CK5NoWk5yJwg0RfjCnLDgPcwrk/dw5KV0rT7ISKsC7xC+9rlsLG2IijU2GXW17cK73eo+XCFn9sGy98mMuQsww+4suF4NFN/msazz+qsNFrhcM9nZBG5ICIJIhItIp+LyPsicjo3Kqdp5th1JoLBvx9gdcAVUg3CiEXGEdofP1Mbayv1YIUm3IDlbyNzOrPwwHU8ZttSrFY7AgKDdIDQCpW7TaY7ShaJ/W4RkTzOuaxpRkcvG/Mz+Z2JoHgRa2JupgBgb/uAHdbHVsKqD7h4+QpDd5XjfKwti5ctpGlTPRpcK3zu1tzUOddqoWkPITBtPsS2k+G4OhUFoNZjjvdfUMxVWDOc1IClTD3mzLhNivc+eI3Fw4djZ2eXk1XWtHzjXpPpNO2R9v36E6alSS/fSGDiuhMUs7Nm2bDm5hciAof+gHUfc/hSDAM3l6JY2crs2r2OGjVqWKjmmpY/3LPjWikVw+1mJzvAFogTkTzIvaxpt12NTmTyZmP3WGkHOyLjkgCoUc6RIjZmNjVFnoOV75JwYgtjD5Vh5l7FhK9G8/rrr6PUA/ZnaFoBYk7HtaOIlEj7ssc4smlqTldEKeWulJqplFqUbpuDUmquUupXpdRLOX1NLX+7tVYEwNx+jU2vbczprDakwq6pMK0Zm/7bSb159px18OHI0UD69++vA4SmpbnvGUAishR4wpxjlVKzlFJhSqmATNs7KaVOKKVOK6U+Siv3rIj0z1REd2CRiAwE9JASLYNbuZm2DW9LPdeSnBz/FN0bVGT8c3XvfuLVQPitPRFLR9JvvT2vr7fnx2mz+Ovvv3nsscdyoeaaln+Y09yUfnV2K4ypObId9ZTJHIxPHfPSlWcN/AR0AIKBfUqp5SISlMX5rsDRtNepZl5TKyQCQ6Jxcy5GZediANjZWPF9L+/sT0i5CdsmItu/Z8FxWz7YaEOvPl0IWPkFjo4P0NGtaYWAOZPpuqR7nQKcB7qaU7iIbFNKuWXa3Bg4LSJnAZRSC9PKyypIBGMMFIfI5qlHKTUIGARQuXJlc6ql5SNnwmOZuPYE3/f2yrBGdUJSKltPhNO2lot5BV3cA8vf4vzpYwzZ5sTlJEeWrZpN48aN732uphVi9wwSItIvh69ZEbiU7n0w0EQp5Qx8AdRXSo0UkQnAYmCqUuoZYEU29ZsBzADw8fEx9wlHyyfWHA1lbeAVnj/jikHA3cWBkBsJ9J25F4Ba98rwejMGNo0lZfcMJh8uypfbrPjgwyH873//w9bWNhfuQNPyN3Oam9yBSYAvxmYmP+C9W08CDyCrHkERkQhgcKaNcUBOByktH7k1B2L7qWvM2XX+jv235kVk6dQGWPEuh05eZMAmB0pUrInf3t+oXr26hWqraQWPOR3XfwJ/A+WBChhThi94iGsGY8z/dIsrEPIQ5WkF1M2UVHadMa6Um1WA8HItSVfvineeGBcB/w4kfk4PRqyJpOO/trw5agKbtmzVAULT7pM5fRJKROane/+7UmrYQ1xzH1BdKVUVuIxxrYoXH6I8rQDaczaC3jN2A+BYxMaUauOWsV09eKWpW8aTRODoIlg7gg2B1xi8wZYmrVtxJGAy5cqVy6Waa1rBYs6TxBal1EdKKTelVBWl1IfAKqVUaaVU6budqJRagLF5qqZSKlgp1V9EUjBmlF0HHAP+FpHAh70RrWC5FSAAWte83Tlds5xxFFKHOpk+9KOC4c9eXPv9dV5ZEs/AzY5M+W0+fy74SwcITXsI5jxJ9E77941M21/H2Efhnt2JItInm+2rgdXmVFArfI4E3zC9rlrGge97eVOnQgmerlueyqWLcflGAuVLpvVFGAzgPxPZMJo/DsXzv82KF195gYBx4ylevHge3YGmFRzmjG6qmhsV0bRbnp26E4BRT9diUKtqAAxtc3vF3EqljfMiCD8Jy9/i3JFdDN5clKtUYMXauTRq1CjX66xpBdU9m5uUUrZKqbeVUovSvoYppfTYQc0iDIbbo5g7e1bI+qDUZNg2kZSfmvHtv3tpNA/avfwB+w4c1gFC03KYOc1N0zAm9fs57X3ftG0DLFUprfC4EZ/EtdgkHi9bnL/3XSIxxTixfkCLqlQolcXw1ssHYPlbHDh0mAHri1DarR57/GdTrVq1XK65phUO5gSJRiLile79ZqXUYUtVSCtc+vy6h2Npy42md8fQ1qR42PIFcdt+YvROK+YH2DLx+x/o27evTsanaRZkTpBIVUpVE5EzYJpcp/MoaTkiqwABUL1cuk7ns1thxTus23eawRtsaP5EBwL+mYqLi5kpOTRNe2DmBInhGIfBnsU4W7oKeha0lgOSUgxZbncqZmtcejThOqz/hPAd83hviy27wkoyfd5sOnbsmMs11bTCy5zRTZuUUtWBmhiDxHERuWnxmmmPFBFBBKzSrdUQfD2eHzeeYlzXuhS1u//1pE+FxQDQpGppSha1ZX3QVVYMa0E915IQtAxZ9T/m7wph+FZF336vcHTcFzg4OOTYPWmadm/m5G6yB4YCLTDOi9iulJouIomWrpz26Gg1cQvVXIozp19johKSafXNFqISkgFo/rgzz9V3ve8yD1y4DsCX3etRzSWteSnmCiwcypldyxm8wYYIm8qs3vg7DRs2zLF70TTNfObMuJ4HeABTMK4NUQeYf9cztAIlMTmVS5EJbD0Rzpyd5/D6fL0pQAAEXM66X+FuwqIT+XSZcaK9m7ODMaXG/rkkT2rEN7OX0WSe0LHfR+w9eFQHCE3LQ+b0SdTMNLppix7dVHiERiXQdMJm0/sxK+5c9mP32Yj7Lnfv+UjTa+vrZ2HFO/jv2sqAtVaUrebN3gPzcHfPdjK/pmm5xJwgcVAp5SsiuwGUUk2AnZatlvao+GbtiSy3P9/QlaOXo4iMSyIoNJqo+GRKFjNvjmVicirD/jyINamsaRxA7I/9+GxLAn8es+bbH6bwkh7WqmmPDHOCRBPgFaXUxbT3lYFjSqmjGNeB8LRY7bQ8F59kzL5avWxxToXFArDszeZ4VSoF3M7Wuvd85J1J97Jx4koMtdUFvradwZl/T/LkekWr9k8RsHQ6ZcqUscyNaJr2QMwJEp0sXgvtkXU2PI72tcvy26uNaPfdVs6Ex+HpWtK0/1awGDjPn7mvN2b8yiBGPVObtjXLZl1gciLWW8fzW/JvvL86Bf/IEvz6xzw6PPlkbtyOpmn3yZwhsBdyoyLaoycxOZUz4bE8VfcxAJa82ZyYxJQMTUH2treHvn646DBXo2+y4nCIKUicuxZH35l7aFDZie9947Be8Q77NgXx3iYDAwcOYvYXX1OsWLHcvTFN08xmzpOEVkjN3HEOg0CdCsYnhxL2tpSwv7Pf4e121Zm86RRXo43TZ/acjUREuBSZQNtvt1KceBrF/MR5v/W8sVY4luBCo3fGMvGL13LzdjRNewDmDIHVCqCUVAN/7btIVHxylvtFhInrjJ3WHhVK3LWs9zvU4I1Wt0ciXb6RQPD1BNYEhPKE1QFW2Qznkt8q6v+WhHO7Idj2mUoLX52tVdPyAx0kCqErUYk8/vEaRvx7lA/+OYyIMP2/M5wNjyU60Rg0gq8nmI53dcoiG2smvu7OAFQoaQ/A9kPHKL/xTQaHf03HX68w6XxVSr46jb0l26CsrKlT/u6BR9O0R4MOEoXMX/su4jthk+n9xmNX2X7qGl+tOc4T3/2H55j1HL8STWCIcYLcrNd8zBqO6uPmhJ21Fd28K9C3qB/N13dmw4attF9gwPf1sfy1+j9sSt4e/eRRoeRdStM07VGh+yQKmRH/Hr1j2yuz9mZ43+nH7abXTd3NG5LqaG/L6teq4Ob3MetOrMV3dSqqckPmbZhFt6a1iE9KoaitNQnJxgTCj5fVS4tqWn6gg0QBdzMllcmbTjGwpTulitll2Pd1j3pZBo30zErcZ0iFfb9RfMlnvLQqmu1Xi0Kn97Fzq09X35oAFLOz4di4TkTE3uRiZDx2NvohVtPyA/2bWoAdvHidmp+s5actZ5i47gQ34pNM+z55pvadC/tkMrSNGau9hR1HZnZk5hfv4vnzDdzbv86uo2eo2bA5v/RteEdTlXPxItSv7PRA96NpWu5TInLvo/KIUsoKGAeUAPxFZO7djvfx8RF/f/9cqVt+8Mzk7aa+hfTm929My+rGBXuSUw0sOxTCs14V+Mv/El29K+A5Zj0A5yY8nX1/REoS7PiBk0u+5o2V8cQVdeXX3xfh5e1tsfvRNM0ylFL7RcQnq325/iShlJqllApTSgVk2t5JKXVCKXVaKfVR2uauQEUgGQjO7brmdzbWWf/3ph9ZZGttRc+GrtjZWNHXtwol7G0Z160uwzvWzD5ABPuT9HNLvhg7mmaz4+k66BP8Dp/QAULTCqC86JOYgzHl+LxbG5RS1sBPQAeMwWCfUmo5xoWO/ETkF6XUImDTncVpWRERLkbEZbnPuXiRu57b17dK1juS4mDzeHYvmsLAVSlUrunJ/qP/UKVKNsdrmpbv5XqQEJFtSim3TJsbA6dF5CyAUmohxqeIS8CthnS9rvZ9CI1K5Hp8MmO7etDLpxJ/7bvEvvORlDIzU+sdzmwm5p9hjFp8mn/P2PHDj9Po9fJrOlurphVwj8ropooYA8ItwRizz04CpiilWgLbsjpRKTUIGARQuXJlC1cz/xi+yLjkh0eFEtjbWvNqMzdebeZ2/wXFR/6/vTsPr6K8Hjj+PQRCIBAWUUCgEiUgkVKXKCKi+U13ZgAAEfVJREFUoIBLRbRSf6i4UqmiiIJaLdqithWtVkWLiqwCBmjBsJQA/iybLEpYAklYE0CyGEAgCUtWTv+YAS+YKwkG7nY+z3OfzLx3ZnLODeTkneV9Yf5QZk+bwBPzS+ne/SZS5k2gYcOGVRuwMcYv+UuRKO/PUVXVw0C/n9pRVUcBo8C5cH0GYgsoxaVHKT16lGXbnImALm5ymk82q0JaAjlTnmHQjGzW5tVlwtTP6Nr9piqM1hjj7/ylSGQCLTzWmwPZPoolIB04XMzYZTsY8eVWOsf88ABcZM3T+BHnZ3N09hDGTElg6OIyHn34QSYMH0GtWqcensMYE1z8pUisAmJEJBrIAvoA9/o2pMCRm19Ih7/9cE1/6da9AHz1h66VO9DRo7BmApsmv8jvE/ZRWKsp/7/4c9pfdnlVhmuMCSC+uAU2HlgBtBGRTBHpp6qlwJPAfGAjME1VU892bIFqQVpuue3N6lfiL//v0yke82tee+5xrh29j96/f5HlGzKsQBgT4nxxd9M9XtrnAnPPcjgBL7+whJcTUn7UPuvJThW786isFFa8z/JPX+PRmQVc2OYS1qbOooXdBGCMwX9ON5nTtGbn/uPLr93Rji3fFTBx5U5aN6576p1zksmf8hgvfraaz7dV5713/0nvBx6121qNMcdZkQhwW3ILAHj5ttjjD8G9dke7n96p5AgsGk7CmLcZOK+Im7t1JXXBVBo0sDGVjDEnsiIRwErLjjJh+U7OrxdBv2ujK7bTjq/Invg4A+M3kZJXh0nTpnB9j9vObKDGmIBlo8AGmINFpbT641xGL82g1dBEsg4cIe9I+VOQnqAwj6MzB/HRk9351Rsbib3hHpLTc6xAGGN+kvUkAkxqVh6lR5W//Gfj8bZjI7p6tWkuG8c+Qf+puyiLbMLCpZ/T7jKbY9oYc2rWkwgwaTknDv3doHYN3rr7V+VvfHA3RZ/dz7DH7qTzyEz69B/CVynfWoEwxlSY9SQCTNpJ80Msfr4rdU5+qloVkuP56qNn6D/je2Jat2Fd6hyaX1DB6xbGGOOyIhFgUrPz6RzTiMeuv4hl2/YSFXHSqK77d5I37Qle+GQBs9KrMeLtd/jNw0/aba3GmNNip5sChKoyYPJq0nLyiT0/ik6tGvH8zRf/sMHRMlgxkhlPXMolzyeiv+hIano2dz0y0AqEMea0WU8iQCxIy2Xuhu8AaH7ycBu5aWRN6M+T475mU34t4uOn0vnW3j6I0hgTbKwnESBGL804vnxTuybOQmkRR7/8CyP7XcWlr66g/fV3sC59txUIY0yVsZ5EAEjNzmPVDmf4jfXDejjXIXZ9Q+pH/eg/cSPUPofFixYQG3eNjyM1xgQb60kEgI8X/9CLiJIiChMG86e+19HlvY307TeApWk5ViCMMWeE9ST83J6CImYlO/MvJd9XjSVDfkn/Kd8S2/oi1q2fTbML2/g4QmNMMLMi4efScvJpQD6Dy8bx3EOLmJsB77/xN+7s/wdfh2aMCQFWJPyZKoVrpjAwfRgvzN9Pr65XkrpgNvXOOc/XkRljQoQVCT/0n/U5XNHgECUJzzDyn/NI+b460yZPoFPP+30dmjEmxFiR8DN5h4pYFv86W9aP52+LDtKhYwe6PvsmnXp29nVoxpgQZEXCn+zdSuo797NwUhKHpTaj46fz9Ipwurds4uvIjDEhym6B9QdlJRR+8Tov3XUpvd5dRd1f3sCBeyYxaHk4qhB7fpSvIzTGhCgrEr6WtYZFz15O+z4vsbmoEfe/Ec+W9oMQCTu+ySVWJIwxPmJFwleKD7N/+hB+17Mj94/bxFuv/pHxi9OZnVOH61ufS8cLzwGgfu0aNK0X4eNgjTGhyq+LhIjcISKfiMhMEenh63iqiqYvYuqjbbnkoXeIaNqG1M3p3P7Eazw2aTUlZcpF59Zh3MNXEl69GpecH2WjuBpjfOasX7gWkbHAbcBuVW3n0X4z8B4QBoxW1eGqmgAkiEgD4C1gwdmOt0odOcC3kwcx4M3P2FFQnenjP6DjXQOOv710614AWjaqTUSNMIbe2pYLzqntq2iNMcYnPYnxwM2eDeKcgP8ncAsQC9wjIrEem7zkvu/3Zidns+9Q8fH1DZl53DlyGbu/nsp7fWK4/KmJXH3NtazZ9t0JBUJVqVerBufWrcl9HS4A4MFrWtKljT04Z4zxnbNeJFR1CbDvpOargG2qmqGqxcAUoJc43gASVXVNeccTkf4ikiQiSXv27DmzwZ9C1oEjDIxfyydLMzhSXMb+Q8U88sEcuq0bQs+77mVGWjHL5s3gpfELCY+sR35hCQ+P+4bEDTlk5xWSd6SEp25oRVg1O71kjPEP/vKcRDNgl8d6JtABGAh0A+qJSCtV/ejkHVV1FDAKIC4uTs9CrF6lZuUB8PHidKYn7eKa/P9w3fLRvLL2MNd1uZb2j44gPDYGgILCEtoPc86eLdy8h+6xjQG73dUY41/8pUiU96ezquoIYMTZDqaydu07THj1aqRm5wPQjFx6b3qbDxI3cUGT+sQ8MoykyLYkrcxi7Mosvhl6I9tyD55wjC/SchGBi5tYkTDG+A9/KRKZQAuP9eZAto9iqbCi0jIen7SG/27aDUAYZdxblEDmos94K6OEjjfdxqoLf8eXQ7ry9/mbSUxxph+dsSaL4YmbAKgdHsbh4jIAohtFElnTX34kxhjjP7fArgJiRCRaRMKBPsAsH8cEOKeOBkxe/aP2fYeKafPSvOMFog07eXTLU4wbPZ5D4efS5JGRfHNhf+pGhBPdKJIb2zY+vu+xAgGQ9urNvNm7PQCxTa0XYYzxL764BTYe6AI0EpFM4M+qOkZEngTm49wCO1ZVU892bOWZm/IdybsOsP9QMQ0iw4+3T1+dCUBNirn30CRWzpvFxHzhH3/9I78d8Aoz1mWzfe8hDheVIiJ0iG74o2Pff7VzF1OP2Mb8uUYYcRc0ODtJGWNMBYmqT6/1Vqm4uDhNSkqqkmPlF5YwanEGHyzcBsBHfa+gXbMonvvXegZ1i2H00gz2py0iNvltPlyym/u6teOZETNp2TLa6zEXbd5Nbn4hf5i+AYCtf72FGmFOZy43v5BzIsOpHuYvnTtjTKgQkdWqGlfee3YCvBxbcgvo8c6SE9pWZnxP0o59rMj4ng2jdvF/ez8mOXEh+2vXYsWsscR0e+iUxz32zEPjqAi25BYcLxDH2owxxt9YkfCQdeAInYb/90ftTetFkJiSQ25+EZ3LviZq+Xu8v66AB3tew5vjE5GIupX6Pl3anGcPyRljAoKd2/BQXoEAuDuuBSX5e3h41zCWj/kLO/Mgrt8wIu/9uNIFwhhjAon1JFxH3NtQPSX/qQfZBw5zeOUEvk58gXE7i+j56x7Mb/EYfTu1ZnCP1j6I1Bhjzh4rEq7NuQUA3HPVL8g+cIR7rvoFUUU5zHmtN0MmJnFju3OJ7vc67/79KbbkFtDWblc1xoQAKxKu1GxnSI0BXS6iRf0Its8czi2/eZWcgjJmj3iOxec/QPe6tQmrJlYgjDEhw4qEKy07n6iI6jQpzODtPr/h9TnbeK5nWwa/P5Ma57XiSl8HaIwxPmBFwrU563vuzRvL1XHTaRBZg5VT3qZVz6fBJvwxxoQwKxJA2YEsGk29m3fW7uXNBzvxwPB/I3Ubn3pHY4wJclYkgIM1GtI6ujnxA1/hhrsf93U4xhjjN6xIAPUiazH403LnNDLGmJBmD9MZY4zxyoqEMcYYr6xIGGOM8cqKhDHGGK+sSBhjjPHKioQxxhivrEgYY4zxyoqEMcYYr4JqjmsR2QPs/BmHaATsraJwAoXlHBpCLedQyxd+Xs4XqOq55b0RVEXi5xKRJG+TgQcryzk0hFrOoZYvnLmc7XSTMcYYr6xIGGOM8cqKxIlG+ToAH7CcQ0Oo5Rxq+cIZytmuSRhjjPHKehLGGGO8siJhjDHGKysSgIjcLCKbRWSbiLzg63iqioi0EJGFIrJRRFJFZJDb3lBEvhCRre7XBh77vOh+DptF5CbfRX/6RCRMRNaKyBx3PajzBRCR+iLybxHZ5P68OwZ73iLyjPvvOkVE4kUkIthyFpGxIrJbRFI82iqdo4hcISIb3PdGiIhUOAhVDekXEAakAxcC4UAyEOvruKoot6bA5e5yXWALEAu8Cbzgtr8AvOEux7r51wSi3c8lzNd5nEbeg4HPgDnuelDn6+YyAfiduxwO1A/mvIFmwHaglrs+DXgo2HIGrgMuB1I82iqdI/AN0BEQIBG4paIxWE8CrgK2qWqGqhYDU4BePo6pSqhqjqqucZcLgI04/7l64fxSwf16h7vcC5iiqkWquh3YhvP5BAwRaQ78Ghjt0Ry0+QKISBTOL5MxAKparKoHCPK8caZfriUi1YHaQDZBlrOqLgH2ndRcqRxFpCkQpaor1KkYn3rsc0pWJJxfmrs81jPdtqAiIi2By4CvgcaqmgNOIQHOczcLhs/iXeB54KhHWzDnC04veA8wzj3NNlpEIgnivFU1C3gL+BbIAfJUdQFBnLOHyubYzF0+ub1CrEg43a+TBdV9wSJSB5gOPK2q+T+1aTltAfNZiMhtwG5VXV3RXcppC5h8PVTHOSXxoapeBhzCOQ3hTcDn7Z6H74VzWuV8IFJE+v7ULuW0BVTOFeAtx5+VuxUJp6q28FhvjtNtDQoiUgOnQExW1Rluc67bBcX9utttD/TPohNwu4jswDlteIOITCJ48z0mE8hU1a/d9X/jFI1gzrsbsF1V96hqCTADuIbgzvmYyuaY6S6f3F4hViRgFRAjItEiEg70AWb5OKYq4d7BMAbYqKr/8HhrFvCgu/wgMNOjvY+I1BSRaCAG54JXQFDVF1W1uaq2xPk5/ldV+xKk+R6jqt8Bu0Skjdt0I5BGcOf9LXC1iNR2/53fiHPNLZhzPqZSObqnpApE5Gr3s3rAY59T8/XVe394Abfi3PmTDgz1dTxVmNe1ON3K9cA693UrcA7wJbDV/drQY5+h7uewmUrcAeFvL6ALP9zdFAr5XgokuT/rBKBBsOcNvAJsAlKAiTh39QRVzkA8zjWXEpweQb/TyRGIcz+ndOAD3NE2KvKyYTmMMcZ4ZaebjDHGeGVFwhhjjFdWJIwxxnhlRcIYY4xXViSMMcZ4ZUXCGB8TkadFpLbH+lwRqe/LmIw5xm6BNeYMcx9gElU96uX9HUCcqu49q4EZUwHWkzAhS0Redudf+MKdj+BZEblIROaJyGoRWSoiF7vbjnfH4V8uIhki0tvjOM+JyCoRWS8ir7htLd15HUYCa4AWIvKhiCS5cyAc2+4pnLGHForIQrdth4g0cpcHu/MlpIjI0ycd+xP3WAtEpNax44lImhvLlLP3aZqg5esnCu1lL1+8cJ5AXQfUwplrYyvwLM4TrDHuNh1whvYAGA/8C+cPq1ic4eUBeuBMQC/ue3Nwhu1uiTMS7dUe37Oh+zUMWAS0d9d3AI08ttsBNAKuADYAkUAdIBVnJN+WQClwqbv9NKCvu5wN1HSX6/v6c7ZX4L+qV0GdMSYQXQvMVNUjACIyG4jAGSTuXx4Td9X02CdBnVNGaSLS2G3r4b7Wuut1cMbM+RbYqaorPfa/W0T644za2hSn2Kw/RYyfq+ohN8YZQGecMXq2q+o6d7vVOIUD93iTRSQBZ3gOY34WKxImVJU3fHI14ICqXupln6Jy9hfgdVX9+ISDO/N3HPJYj8bpqVypqvtFZDxOUapsjOXFUobTIwJnwqXrgNuBl0XkElUtPcX3McYruyZhQtVXQE9x5kWug/PL9TCwXUR+C84FZxH51SmOMx94xD0GItJMRM4rZ7sonKKR5/ZCbvF4rwDnlNfJlgB3uCOdRgJ3Aku9BSIi1YAWqroQZ+Kl+jg9G2NOm/UkTEhS1VUiMgtnTuCdOCOo5gH3AR+KyEtADZx5KZJ/4jgLRKQtsMI9RXUQ6Ivz173ndskishbnukIGsMzj7VFAoojkqGpXj33WuD2OY0Naj1bVtW4vpTxhwCQRqYfTC3lHnWlMjTltdgusCVkiUkdVD7rPKCwB+qs7J7gxxmE9CRPKRolILM61gQlWIIz5MetJGGOM8couXBtjjPHKioQxxhivrEgYY4zxyoqEMcYYr6xIGGOM8ep/d23Nmmtpt2oAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure()\n", "plt.plot(bh1.pop_hist)\n", "plt.plot(pop_avg, label='average growth rate')\n", "plt.plot(pop_theory, 'k', lw=1, label='theoretical growth rate')\n", "plt.yscale('log')\n", "plt.xlabel('generations')\n", "plt.ylabel('population size')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "overall-argument", "metadata": { "slideshow": { "slide_type": "slide" }, "tags": [] }, "source": [ "## Optimal proportion of persisters" ] }, { "cell_type": "markdown", "id": "collaborative-seattle", "metadata": { "slideshow": { "slide_type": "subslide" }, "tags": [] }, "source": [ "Let us now look for the optimal proportion of persisters in the cell population, which maximizes the long-term population growth rate. To find that, we need to simulate the population growth for different proportions of persisters." ] }, { "cell_type": "code", "execution_count": 9, "id": "aggressive-smooth", "metadata": { "slideshow": { "slide_type": "fragment" }, "tags": [] }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":10: RuntimeWarning: divide by zero encountered in log\n", " lam_avg = np.log(bh1.pop_size / N0) / T\n" ] } ], "source": [ "q1_list = np.arange(0.01, 1, 0.01) # list of values for q_1\n", "lam_avg_list = [] # list to collect average growth rates\n", "lam_theory_list = [] # list to collect theoretical growth rates\n", "\n", "for q1 in q1_list:\n", " phe_dist = [1-q1, q1]\n", " bh1 = BetHedging(phe_dist, fit_mat, N0=N0, record=True)\n", " env_seq = np.random.choice(2, size=T, p=env_dist)\n", " bh1.grow(env_seq)\n", " lam_avg = np.log(bh1.pop_size / N0) / T # simulation result\n", " lam_theory = np.dot(np.log(np.dot(phe_dist, fit_mat)), env_dist) # theoretical result\n", " lam_avg_list.append(lam_avg)\n", " lam_theory_list.append(lam_theory)" ] }, { "cell_type": "markdown", "id": "piano-young", "metadata": { "slideshow": { "slide_type": "subslide" }, "tags": [] }, "source": [ "We can plot the population growth rate with respect to the proportion of persisters." ] }, { "cell_type": "code", "execution_count": 10, "id": "meaning-sentence", "metadata": { "slideshow": { "slide_type": "fragment" }, "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEICAYAAAB1f3LfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxU1dnA8d+Tyb7vIUAgAdl3iKwquKBIVcC1rlSraK1W21qlVVv7vm9btG51LyqKiLsiKCgCLSACQoAAYU8gQEiAJGRPyDJz3j9miAESSEgmN8vz/XzmMzP3nnvnOSHMk3PuueeIMQallFLqXHlYHYBSSqnWTROJUkqpRtFEopRSqlE0kSillGoUTSRKKaUaxdPqAKwQGRlp4uPjrQ5DKaValQ0bNuQYY6JO3d4uE0l8fDxJSUlWh6GUUq2KiOyvbbt2bSmllGoUTSRKKaUaRROJUkqpRmmX10hqU1lZSUZGBsePH7c6lDbH19eXzp074+XlZXUoSik30ETikpGRQVBQEPHx8YiI1eG0GcYYcnNzycjIICEhwepwlFJuoF1bLsePHyciIkKTSBMTESIiIrSlp1QbpomkBk0i7qE/V6XaNk0kSqlGKa+yM29TBknpx+p9jMNhyCooc2NUqjlpImnh7r77brZv394k54qPjycnJ+eMZf7+97+f9H706NFN8tmq9Ssoq+RvC7fzj0U7+HLTIVIOFfDqf1O54On/8tuPN3PLWz/yQ+qZf78Ayirs3Pv+Bkb94z88990u7I6T10TKKiijospx0rZKu4OtGQXkFpc3aZ1q+u+uo1z5r++ZvTodXaepYaQ9/sASExPNqXe279ixgz59+lgUUfM4cUd/ZGRknWUCAwMpLi5u8s9uDz/ftmzX4SLunZPEwbwybCJU2H/6or+wRyS3j+jC60s2k38sh5emdGdApAeUF0F5IVSUuB7FlJYUsmxLOiXFhcQFQlFxMbEB0DfKm5LSUo4VFFJZUY63VOFvM/jZ7Bh7FcZehQd2bDjw8jDYxOCBQWp8fzlfCQ4RDDYcCEY88LB54enlBR6elDts5Jcbiio98PXzIzw4EF8/f/YeqyIt347d5kNhlRcdo8IZ3ScOL98gHF4BHC6zERQSRlBwGPgEOR++IeAbDN6B0E66b0VkgzEm8dTtOmqrBSkpKeHGG28kIyMDu93Ok08+yeuvv86zzz5LYmIigYGB/PrXv2bp0qWEhYXx97//nUcffZQDBw7w4osvcs011/Duu++SlJTEK6+8AsBVV13FI488wrhx4076rMmTJ3Pw4EGOHz/OQw89xLRp05g+fTplZWUMHjyYfv36MXfu3OrEYozh0Ucf5ZtvvkFEeOKJJ7jppptYvnw5Tz31FJGRkaSkpDBs2DDef/99vS7STAqPVxLk43lOP+9jJRXsOVLE4C6h+HjaftphDJTlQUkOlGSTtH03C9du5eeeJUwZ5EuUrYTSgmwqi3IJcBThnVMInxdwubGDDVhQ92d6Y2Os8cE7IAAf30AKbDb2FxqS93tx3Hji4RVKh+gQSh020ksd5JQ68PXxISo8gKhgf3LL7KRml1Jw3I4DISLQl46hftgdhv25JZSUVzoTDMaZcIwdjANfmyHA01BRUU6gp4OYIA8OFZWQXlKBv60MT0c5QwMh0tdQVlKE/VgpsrocsOMBdDzDz9GBjUqvIMQvFPEPp8wziCIJptw7lLhOnfEKigL/CAiIotIvgtRiP2JiOhAe6OP6cRuOFpWTdrSYfh1DCPFvfcPkNZHU4q9fbWN7ZmGTnrNvx2D+cnW/M5b59ttv6dixIwsXLgSgoKCA119/vXp/SUkJ48aN4+mnn2bKlCk88cQTLFmyhO3btzN16lSuueaaescza9YswsPDKSsr4/zzz+e6665jxowZvPLKKyQnJ59W/osvviA5OZnNmzeTk5PD+eefz0UXXQTApk2b2LZtGx07dmTMmDH88MMPXHDBBfWORZ2bYyUVXPLcckYkhPParcOwedSeTDYdyGN1Wi7dIwPoE244nnOQFRs2k75vD5GOY6R7FtA3qJQuXoUE248hxUfBUVl9fCKQ6AE4gL3B4B9OoF8YRESDXy/wCwXfUPANocD48cL3h9lXZKPI+FHm4Y/dK4C8Kh/yq7wJDQrgjduGMaxrGAChwL4Debzz/V6u7B/Llf074Glz9rj3w/klWzNJdgfON4bNGQV8vzubr/blkpSeh6eHcGGPKC7pE01i1zCig30J8LZRWmFndVouy3cdJSOvjClDOjFuQAd8PG0UHa/ki42H+CE1h9tHdWVQD+dchAHAit3ZPPrZZoI8HYyN9+P8jt4cPJzNmh37sJcVEUgZQVJKuK2MECnFv6yIkOMlhOYVEyIHCaOISCnGK+3k60BeQB+g3HhyWEIp8gwn0x7MwcoQjpgwFtki6NmjF5cMH0znruc5Wz2tgCaSFmTAgAE88sgjPPbYY1x11VVceOGFJ+339vZmwoQJ1WV9fHzw8vJiwIABpKenN+izXnrpJebNmwfAwYMH2bNnDxEREXWWX7VqFTfffDM2m42YmBjGjh3L+vXrCQ4OZvjw4XTu3BmAwYMHk56eronkDD5ef4AN+/N4+rqB9W5JHK+04+tlO2nbW9/vJb+0ksXbjvD0tzv505W9oSQb8vZDvvOReyiNoh3bGG+y6Si5BIpzGHYvcLYebFBsC+FQcSjJ9hBKvHvRpesVxHTswuwtJWzM9eKiQX246/JEvIMiwdP7jHGGAI8Nt3Mwr5SoQB9C/LzwcCW4E9dCTk14Q7qE8dqtw2o9X20/HxFhcFwog+NCeZAeVFQ5EAEv2+mXfAN8PBnfN4bxfWNO2xfk68XU0fFMHR1/2r6xPaNY+8dLT/v8X9gdfL8nm7IKB31ig+gaEYAA+3JLSDlUwN6SCjqF+eMd5seh0ko+XJPKuu2phJkCYjyLGdsREqOqcBQdobLgMJ5l2fSw5THcMw2/yjznh6S5HkCx8SPThJNti6LMryOExOEX3Y3IuJ506d4Xv9CYFtGt1iISiYhMAP6F81f7LWPMjFP2i2v/RKAU+IUxZmN9jj0XZ2s5uEvPnj3ZsGEDixYt4o9//COXX375Sfu9vLyqf7E9PDzw8fGpfl1VVQWAp6cnDsdP/de13b+xfPlyli5dypo1a/D392fcuHFnvc/jTNfSTsQBYLPZqmNRpztadJy/frWd0go7Y3tG87OBsbWWszsMi7ZmsTothx/3HmNfbgnP3TCIa4d0gqIsig7tpHT1It6KLaSTIwvW7qNyQw5e9pP/ArYRRLRHFJ0TBlDiG0uaCafYpwOD+/clICIOgjoQ6OlDN7uDfTuOMHPFXpJ35MMOCPLx5J+3DGRC/9pjrIuft42eMUGnba+rxdRY3p7uGTNUWxLzsnlwSe/Tk1L3qEC6RwWetn1U9wgOFwxic0Y+I7tFEOJ3hm6rqnIoOkze4XTWbUnBVpRJaOVRQiqO0LUsi+CyVQSXFkIWsNl5SCm+ZHt2oMC3M2VBXakISYDwbvjF9qJH956E+J858TcVyxOJiNiAV4HxQAawXkQWGGNqDlW6EujheowAXgdG1PPYViMzM5Pw8HBuu+02AgMDeffddxt8jvj4eF577TUcDgeHDh1i3bp1p5UpKCggLCwMf39/du7cydq1a6v3eXl5UVlZedp0JhdddBH//ve/mTp1KseOHWPlypX885//ZOfOnQ2OsT17eVkq5VUOukb48/S3OxnfN+a0L0JjDH+dt4k1SesZ4HOYB0OPER60j8j5B3B8cxSPyhKCgKc8wFHog4TFszmgE+8X96NH7wF06Nqb8E49uOvLw6Tlw2f3jca/QxD+wGkLSbh42TyY0D+WK/p1YH16Hit2H+XGxDi6RgS4+0fS5nUI8aVDSIezF/T0gbCuhIV15Yo+Y2stYsqLOHwglaz0HRRkpiL56QSWHiSiJJ3eRWvxzvrpj7hi48tuWyeKAhNwRPTAv1NfYhIGEtG1D+LpU+v5z5XliQQYDqQaY/YCiMhHwCSgZjKYBLxnnH8WrxWRUBGJBeLrcWyrsXXrVv7whz/g4eGBl5cXr7/+Oo888kiDzjFmzBgSEhIYMGAA/fv3Z+jQoaeVmTBhAm+88QYDBw6kV69ejBw5snrftGnTGDhwIEOHDmXu3LnV26dMmcKaNWsYNGgQIsIzzzxDhw4dNJE0QHpOCR+uO8DNw+O4rE8Mv3hnPXPWpPPL/jY4ss35OLqNvPQtPFmcjpeP3XlgPtiD4/ixIoIF0p9LLhnD75cVE5PQl/+7YwJ42Oh+vJKn3l5H8tZ82AqQgc1DeOcX59Orw+mtg7qICMMTwhmeEO6OH4FqJPEJIrbHEGJ7DDltn7FXUZi9n5LMXRRn7qI0cye2vDRiCzfTsXAp7ANWwcbRrzL08tuaNi6rh/+KyPXABGPM3a73twMjjDEP1CjzNTDDGLPK9X4Z8BjORHLGY2ucYxowDaBLly7D9u8/eX0WHZ7qXq3h52uM4aVlqXy9JZP37x5BTLDvWY+xOwx/+mIr1w7txIhudV9jAnjog/Wk79jInIneBOVtZ1fyajpVpBHET91RpQGdWV0YTVVELy4fNw6P6N4Q2QO8A1i+6yi/eGc9HYJ9OVx4nG8eupA+scEnxZ9VcJyUQwVsyyxkaNcwxvasqw2i2pO8vDwyUjeTfyCF3qOvISq2yzmdpyUP/62t4/TU7FZXmfoc69xozExgJjjvI2lIgKrtM8Yw45ud/HvlXgB+/8lm3rtrePWF4rqs2H2Uj5MOsmJ3Nt/97iKCfV1dgg4H+Qe3UZD2IyHHtmI7nMwz2dvwsVXCYsDLny7hffg8cwy+nQcR0GUwm8s78N6GY/SMCeSjaaPw8D754vq4XtHcN7Y7b6xI44p+MSclEXC2JjqG+tEx1I/L+9WjK0W1G2FhYYSdPw7OH+eW87eERJIBxNV43xnIrGcZ73ocq9QZORyG//l6O++uTuf2kV3pHRvE4/NSeHvVPu65qNsZj/3gx4ME+XpSXpTDZx/O4q6u2ZCxHvvBDYRWFhGKs686xSSw22MC11/9M/y7DIOI7vh72Ej+ZDOfb8yAfeDvfYwBnUN45ZYh+J2SRE74/eU9CQ/w4qqBZ7qzQanm1RISyXqgh4gkAIeAnwO3nFJmAfCA6xrICKDAGJMlItn1OFap0xyvtPPjvmOs2pPNyt057DpSxN0XJPD4z5zdbyt3Z/PM4p2M6h5B/061jOXP20/+zpVckvoFMwL3Ekk67AdzwEZ5RB++qhzJLs9ejL14Asf8u1JSAYM7BePfOfSk0/zf5P7cMqILcWF+RAX5nHU4sJfNg2kXdW+qH4NSTcLyRGKMqRKRB3A2+G3ALGPMNhG5z7X/DWARzqG/qTiH/955pmMtqIZqJarsDj7dkMGLS3dzpLAcb5sHw7qG8fcpA7h5eFz1F/mMawcy4V8r+c1Hm1j44IX4lWTAvpWQvgr2/wAFBwkFrvLwxyt6FJUJtzF9nR8bKhMoyvPG01v4eNoo4iPPPOrJz9tWfXOeUq2V5YkEwBizCGeyqLntjRqvDfDr+h6r1KkcDsOilCye/243e3NKGNIllL9NHsDo8yLw9z79v0EYRbw3/CCbln9J+XO/wq/c1WPqHwnxY3CMepBf/McLonvz3p3OiS1v7ZbHvNdXExEofHjPyLMmEaXaihaRSJRyF7vDsHBrFi8v28Oeo8X0iA5k5u3DGN835uRuJHslZKyHPUsgbRlkbaEXhi7egXxf1pshF91H1MDLIao3iLBi51FWFq7ntat/WvVxaJcwPrhnJJ3D/Ogc5m9BbZWyhiaSFiI/P58PPviA+++/n+XLl/Pss8/y9ddfWx1Wq7Rh/zGW78pme2YhWw4VkF1UTo/oQF66eQg/GxD70x3WxdmQugR2fwtp/3XOVCs2iBsBFz8O3S+mJLgvjzz/PX33BvPhpb2rk88H6w4QGehz2tQbI88yBFiptkgTSQuRn5/Pa6+9xv333++2z6iqqsLTs23/k6ceLeaGN9YgIpwXFcgF50VyWZ8YruzfwTmUN3s37FoIOxc5WyAYCOwAfSdBj8uh29iTJsqLBB6d0Jsnvkzhy+RDDE+IYPbqdP6z8yjTLupW6/xOSrU3bftbpRWZPn06aWlpDB48GC8vLwICArj++utPm5p9w4YN/O53v6O4uJjIyEjeffddYmNjSU5O5r777qO0tJTu3bsza9YswsLCGDduHKNHj+aHH37gkksu4d1332X37t14eXlRWFjIwIED2bNnz2lTorRWLy7dja+XjeV/GEd0kK9zSvSsZPjvm7B9AeTucRaMHQzj/gg9r4DYQWec+O6W4V34dEMGj89Lody14NLEAbHce5ahwUq1F5pIavPNdDi8tWnP2WEAXFn3fJIzZswgJSWF5ORkli9fzqRJk06bmn3EiBE8+OCDzJ8/n6ioKD7++GMef/xxZs2axR133MHLL7/M2LFj+fOf/8xf//pXXnzxRcDZ2lmxYgUA6enpLFy4kMmTJ/PRRx9x3XXXtZkksj2zkK+3ZPHAuO5EF+2EtV/AtnmQf8DZZRV/AYy4F3pNhJBO9T6vh4fwjykDeOijTVzcO5qpo+PpFOrnxpoo1bpoImmhapuaPTQ0lJSUFMaPHw+A3W4nNjaWgoIC8vPzGTvWOdHb1KlTueGGG6rPddNNN1W/vvvuu3nmmWeYPHky77zzDm+++WYz1qrhSiuqWLErmy4R/pwXHXjyAkynmLtoGY/5fs20XRth7V7w8IRu4+CiR53JI+Dcr1/07RjMkt/VPpGeUu2dJpLanKHl0Fxqm5rdGEO/fv1Ys2bNSWULCgrOeK6AgJ+GoY4ZM4b09HRWrFiB3W6nf//+TRt4E3tz5T5eWLobcE5D3ic2iBduHEyPE9OUl+RAyueUrH+fv+VswSBI6IVwwUPQ5xrw18kHlXI3vVLYQgQFBVFUVHTGMr169SI7O7s6kVRWVrJt2zZCQkIICwvj+++/B2DOnDnVrZPa3HHHHdx8883ceeedTVcBNzDGsGDzIQbHhfLyzUP41djuHC4oZ+pbq8lePw8+vAWe6wXfPMrRglKek6mUPrAVpn4Fw36hSUSpZqItkhYiIiKCMWPG0L9/f/z8/IiJOX3xHG9vbz777DN+85vfUFBQQFVVFQ8//DD9+vVj9uzZ1Rfbu3XrxjvvvFPnZ91666088cQT3Hzzze6sUqPtPFxEWnYJ/ze5P1cP6sjVcXu58/h/cGyaS9TCfOz+UWT3vYu/ZQziq8NhPHlVXwIi485+YqVUk7J8GnkrJCYmmqSkpJO2tYZpzpvKZ599xvz585kzZ06zfWZdP9/i8ioWbc3ii40ZBHh78uYdidUz7j7z7U7eXrmbpOsrCdr2PqT9B8RGQdzFPJ4+hHWewzha6iAm2IfpV/Zm8uBO9V66VinVcC15GnnVjB588EG++eYbFi2yflaZt1ft45+Ld3K80kF0kA9Hi8pZsDmTyUM6YYqPErbhJdb4LSZoQQ4Ed3LeJDjkdkKCY7ltby7Jn27m/uEd+fXF5xHgo7/KSllF//e1My+//LLVIQBQUFbJs4t3MahzKI9O6M3guFAmvbqKL75ZzFX71uGR8jn3OCo4HDkKLnsJelwBtp9+XUd2i2DVY5dYWAOl1AmaSGowxmjXiBvU1n36+YYMyirtPHlVX/p3DIa0Zbzv9TyhFT9QmeLLxsirePzQaD656xfg3zbuc1GqrdJE4uLr60tubi4RERGaTJqQMYbc3Fx8fX9attbhMLy/dj/D4oLof2wJfPUiHN5KaFAsn4b+khfzx1CVH0L/HiGEaBJRqsXTROLSuXNnMjIyyM7OtjqUNsfX17f65kqAH3ZnMSxvIX8x38LnByCyJ0x6FQbcyMCccrL+tRLH8XKmXxlrYdRKqfrSROLi5eVFQkLC2Quqc2evhOS59P7mH1zodRhH0CC46u/Q62fg4bylqVcHb246vwtfbc7ksj6nD4FWSrU8mkiU+znssOUTWDED8tLJcJzHir7Pc/1Nd9U6WeL/TOrHQ5f2IMhXu7WUag00kagmU15lZ2tGAevT80g5VMAdI7swomItLPsfyNkFHQbySa/nmb4lhu8nXFrnjLteNg86hPjWuk8p1fJYmkhEJBz4GIgH0oEbjTF5p5SJA94DOgAOYKYx5l+ufU8B9wAnLmz8ybX0rnKTzPwyXlueyr0XdScu/KdVABdvO8wjn2ymqLwKgFFeqXjvmQvsgogecMNs1viM4Yl3k5jQP1pnz1WqDbF6rq3pwDJjTA9gmev9qaqA3xtj+gAjgV+LSN8a+18wxgx2PTSJuNlH6w7w/toDTHzpe77ZmoUxhlf/m8q9czbQLTqQ2dfGsGvQx3xo+zNxHOE5n/vJv3MlW0Mu5p45G+ka7s/fJg+wuhpKqSZkddfWJGCc6/VsYDnwWM0CxpgsIMv1ukhEdgCdgO3NFqWqtio1hx7Rgfj7ePKruRvp1zGYbZmFXDcgjKejl+K5+BUQDxg7nQOdb+ff725l9Zxk9uWUEOrvxZxfjiAswNvqaiilmpDVLZIYV6I4kTCiz1RYROKBIcCPNTY/ICJbRGSWiISd4dhpIpIkIkk6xPfcFB6vZHNGAVf068Cn947ingsT2Hm4kFeHHOLZo/fg+cNz0PcaeDAJLv4jQ3t05h/XDmDD/jw8BOb8coRe+1CqDXJ7i0REluK8vnGqxxt4nkDgc+BhY0yha/PrwP8CxvX8HHBXbccbY2YCM8E5aWNDPls5/bj3GHaHYcx5kXh7evD4KD+m576NbccyiO4Hd74JXUefdMx1wzoT7OfFedGBJEQG1HFmpVRr5vZEYoy5rK59InJERGKNMVkiEgscraOcF84kMtcY80WNcx+pUeZN4Oumi1yd6ofUHHy9PBgaFwirXoTlM7B5eMIV/4Dh006aC6um8X31fhCl2jKrr5EsAKYCM1zP808tIM75St4Gdhhjnj9lX+yJrjFgCpDi3nDbt1WpOVzfKQ+fWZfB4S3OGwkn/rNB658rpdoeqxPJDOATEfklcAC4AUBEOgJvGWMmAmOA24GtIpLsOu7EMN9nRGQwzq6tdODeZo6/TVmdlkN2UTmTBp+eGA4fK2Ji7ns8XPSlc+XBG99zLmWr85Ip1e5ZmkiMMbnApbVszwQmul6vAmr9tjLG3O7WANuRgtJK7p+7kfzSSrxsHkwcUGOeq6M78Zl7J7/z2k5+98mEXvuCLmOrlKpm9agt1UK88t89FJRV0jMmkEc+3cyuw0VgDKx7E2aOxbM4k0c8HiH4lnc1iSilTqKJpB16dvEu3luTXr1OyIHcUmav3s8Nwzoz55cjCPDx5JHZy6h8/yZY9Agm/kJu9HiO4+f9rHoZXKWUOsHqaySqmaUeLeKV/6a6Xhfzl6v78fS3O7F5CL+/vBcxwb7MGe8gdOHDmLRi1vT8A8f63cWOlE1MPS/S4uiVUi2RJpJ25qN1B/GyCTcmxvHemv3syCpkfXoeD13ag5ggH/jhX/Re+leKAzvxK/uTLNvSAbZsAmCMJhKlVC00kbQj5VV2Pt+Ywfi+MfxtygB6xgTx1FfbiA7y4d6RkfDxbbDza+g7icBrXuFt32COFB5n3b5jVFQ5TpqkUSmlTtBE0o4s3naEvNJKfn5+FwCmjo6nX8dgQssO4j97AuSmOm8uHPmr6mG9McG+XD2oo5VhK6VaOE0k7chH6w7QOcyPC2p0USVWbYT5d4HY4I4vIeEiCyNUSrVGOmqrnUjPKWF1Wi4/Pz/up5FX69+CuTdASBxMW65JRCl1TrRF0k58tP4gNg/hhsQ4cDhgyZOw5hXocQVcPwt8Aq0OUSnVSmkiaQeq7A4+25DBxb2iifEDPp0KOxbA+ffAhBl1TraolFL1od8g7cAPabnkFJdz08BQ+OAG2LcSLv8bjPq1zpWllGo0TSTtwILkTLr4lHLpuruds/ZOmQmDbrI6LKVUG6GJpI07Xmln47YdfOrzNzyyj8LPP4BeE6wOSynVhmgiaePWJG/lLfMUkY4CuO1ziL/A6pCUUm2MJpK2rDCTvt/dSqBHHubWLyB+jNURKaXaIL2PpK0qOozjnZ8RUJHLBz1exDNBk4hSyj00kbRy2zILmPZeknP9kBNKj8GcKTgKs7ij4jGGjLnCugCVUm2edm21YkeLjnP37CSyCo6zOi2XV24Zwrh4f/jgRshN5fmI/+WIV0+GdgmzOlSlVBtmaYtERMJFZImI7HE91/qNJyLpIrJVRJJFJKmhx7dFxyvt3DtnA/mllbw9NZEu4f7c9+5qDr5xHY6MDbwS/ifeOBjHVYNidTEqpZRbWd21NR1YZozpASxzva/LxcaYwcaYxHM8vs0wxvCneVvZdCCf528cxKV9Yvj03pHMiphDXN5aHq28h89KB3PbyK7cP/Y8q8NVSrVxVndtTQLGuV7PBpYDjzXj8a3Sgs2ZfLHxEL+9rCdXDogFIGDtc4wuXsr+Qb/lobGP6tohSqlmY3WLJMYYkwXgeo6uo5wBvhORDSIy7RyOR0SmiUiSiCRlZ2c3UfjW+HzjIeLC/fjNpa7WxuaPYfk/YNAtdJ38F00iSqlm5fYWiYgsBTrUsuvxBpxmjDEmU0SigSUistMYs7IhcRhjZgIzARITE01Djm1J8ksrWJ2awy8vTEBE4MBamP9riL8Qrv6Xzp2llGp2bk8kxpjL6tonIkdEJNYYkyUiscDROs6R6Xo+KiLzgOHASqBex7cl3207QpXD8LMBsVCYBZ/cAaFxcNMc8PS2OjylVDtkddfWAmCq6/VUYP6pBUQkQESCTrwGLgdS6nt8W7Nwaxadw/wY0MHPOR18eTHcNBf82s2ANaVUC2N1IpkBjBeRPcB413tEpKOILHKViQFWichmYB2w0Bjz7ZmOb6vySyv4ITWHnw2IRRY/Dgd/hEmvQExfq0NTSrVjlo7aMsbkApfWsj0TmOh6vRcY1JDj26ol253dWrf4rYEVb8LoB6H/tVaHpZRq56xukagGWLQ1i+Eh+XRZ8yR0GQ2XPmV1SEopZfl9JKqeCsoqWZt6mKWhryB2G1w7U5fIVUq1CPpN1Ep8t+0wDyX6iZAAABO7SURBVMindCrZDjfMdo7UUkqpFkC7tlqBKruDtf/5kl95foUZcgf0m2x1SEopVU0TSSuwYN0eHi7+F2WBXZAr2/TANKVUK6RdWy3c8Uo7lUv/SpxHNub6ReAdYHVISil1Em2RtHCLv/mSG6oWcbjX7YgulauUaoE0kbRgBUVFDNz4BMc8o+hw7T+sDkcppWqliaQF2/Hxn0kgk+LLnwefIKvDUUqpWmkiaaFM7l6GZrzHD/6XEj/iaqvDUUqpOunF9haqcP6j2IwnuaMbMtu+Uko1P22RtER7lhJyYAmv2adw0dABVkejlFJnpC2SlqaqAvPtdDIklp0JtxHqr2uMKKVaNm2RWGR7ZiGlFVWn71j3byR3D38uv40rBnZp/sCUUqqBNJFYIL+0gptmruHeORsor7L/tKMsH1b+k70hI1nJUMb3rW2FYqWUalk0kVgg1N+bJ6/qy/d7cnjwg01U2R3OHWteheMF/L38RkZ3jyA8QLu1lFItnyYSi9yYGMdfru7Ld9uP8IfPtuAoyoa1r1GQMJGl+R24sn+s1SEqpVS96MV2C905JoGS8iqe/W43o9Je4PqKUv5aNBmbh3BFvxirw1NKqXqxtEUiIuEiskRE9riew2op00tEkms8CkXkYde+p0TkUI19E5u/Fo3z64vP49krophcuYglnmP5KiuIy/pEExHoY3VoSilVL1a3SKYDy4wxM0Rkuuv9YzULGGN2AYMBRMQGHALm1SjygjHm2WaKt8mJCNeXfgzi4Ir7X2BnaDweYnVUSilVf1ZfI5kEzHa9ng2cbcWmS4E0Y8x+t0bVnIqOwMb3YMhtEJ6AzUMQ0UyilGo9rE4kMcaYLADXc/RZyv8c+PCUbQ+IyBYRmVVb19gJIjJNRJJEJCk7O7txUTeldTPBXgmjf2N1JEopdU7cnkhEZKmIpNTymNTA83gD1wCf1tj8OtAdZ9dXFvBcXccbY2YaYxKNMYlRUVHnUBM3qCiB9W9Bn6sgorvV0Sil1Dlx+zUSY8xlde0TkSMiEmuMyRKRWODoGU51JbDRGHOkxrmrX4vIm8DXTRFzs9k0F47na2tEKdWqWd21tQCY6no9FZh/hrI3c0q3liv5nDAFSGnS6NzJXgVrXoG4ERA33OpolFLqnFmdSGYA40VkDzDe9R4R6Sgii04UEhF/1/4vTjn+GRHZKiJbgIuB3zZP2E1g51eQv19bI0qpVs/S4b/GmFycI7FO3Z4JTKzxvhSIqKXc7W4N0J1Wvwzh3aHXlVZHopRSjWJ1i6R9ykyGQxtgxH3gYbM6GqWUahRNJFbYNAc8fWHgjVZHopRSjXbOiURExojIq00ZTLtQUQpbPoW+k8Ev1OpolFKq0Rp0jUREBgO3ADcCBzj5ng5VHzsWQHkBDG29l3eUUqqmsyYSEemJ847ym4Fs4DNgtOuCuGqoje9BeDfoOsbqSJRSqknUp0WyE1gIXGGMOeDmeNq2nFTY/wNc9hTofFpKqTaiPtdIrgNKge9F5E0Rudw1C69qqE3vgdhg0C1WR6KUUk3mrInEGDPPGHMT0BdYDvwGyBCRt0RkgpvjazvslZD8gfO+kSBdtEop1XbUe9SWMabEGDPXGHMV0A9YDzzqtsjamvRVUJINg262OhKllGpS5zT81xhzzBjzb2PMJU0dUJu1fT54B8J5p93Ir5RSrZrekNgcHHbY8RX0vAK8/KyORimlmpQmkuawfzWU5kDfBi3BopRSrYImkuawfT54+cN5462ORCmlmpwmEndzOJx3s/cYD97+VkejlFJNThOJux38EYqPQJ9rrI5EKaXcQhOJu22fDzYf54V2pZRqgzSRuNOJbq3zLgOfIKujUUopt9BE4k5ZyVB4CPpqt5ZSqu2yNJGIyA0isk1EHCKSeIZyE0Rkl4ikisj0GtvDRWSJiOxxPYc1T+T1tHe587m73replGq7rG6RpADXAivrKuCaIPJV4Eqc833dLCJ9XbunA8uMMT2AZa73Lcfe5RDTHwKjrY5EKaXcxtJEYozZYYzZdZZiw4FUY8xeY0wF8BFw4s6+ScBs1+vZwGT3RHoOKsvgwFroNs7qSJRSyq2sbpHURyfgYI33Ga5tADHGmCwA13Odf/qLyDQRSRKRpOzsbLcFW+3AWrCXayJRSrV5DVpq91yIyFKgQy27HjfGzK/PKWrZZhoahzFmJjATIDExscHHN9je5eDhBV1Guf2jlFLKSm5PJMaYyxp5igwgrsb7zsCJZX6PiEisMSZLRGKBo438rKazdznEDQefQKsjUUopt2oNXVvrgR4ikiAi3jjXj1/g2rcAmOp6PRWoTwvH/UqPQdZm7dZSSrULVg//nSIiGcAoYKGILHZt7ygiiwCMMVXAA8BiYAfwiTFmm+sUM4DxIrIHGO96b719KwGjiUQp1S64vWvrTIwx84B5tWzPBCbWeL8IWFRLuVyg5a0UtXc5eAdBx6FWR6KUUm7XGrq2Wp+9yyHhQrBZmqeVUqpZaCJpannpkLdPu7WUUu2GJpKmduBH53PXMdbGoZRSzUQTSVPL3OhcDTGqt9WRKKVUs9BE0tQyN0HsIL0+opRqNzSRNCV7FWRtgY5DrI5EKaWajSaSppS9E6rKNJEopdoVTSRNKXOT81nvH1FKtSOaSJpS5ibwCYbwblZHopRSzUYTSVPK3AgdB4OH/liVUu2HfuM1lapyOJyi10eUUu2OJpKmcmQbOCo1kSil2h1NJE1FL7QrpdopTSRNJXMT+IVDaBerI1FKqWaliaSpZG6CTkNBalsZWCml2i5NJE2hohSO7tDrI0qpdkkTSVM4kgLGrtdHlFLtkiaSpnB4i/M5dqC1cSillAWsXrP9BhHZJiIOEUmso0yciPxXRHa4yj5UY99TInJIRJJdj4m1ncPtcveCpx8Ed7Lk45VSykpWz3WeAlwL/PsMZaqA3xtjNopIELBBRJYYY7a79r9gjHnW3YGe0bE057QoeqFdKdUOWZpIjDE7AOQMX8DGmCwgy/W6SER2AJ2A7XUe1Nxy0yBaF7JSSrVPreoaiYjEA0OAH2tsfkBEtojILBEJO8Ox00QkSUSSsrOzmy4oe5Vznfbw7k13TqWUakXcnkhEZKmIpNTymNTA8wQCnwMPG2MKXZtfB7oDg3G2Wp6r63hjzExjTKIxJjEqKuoca1OLgoPOqVEiNJEopdont3dtGWMua+w5RMQLZxKZa4z5osa5j9Qo8ybwdWM/q8GOpTmftUWilGqnWnzXljgvoLwN7DDGPH/Kvtgab6fgvHjfvHL3Op+1RaKUaqesHv47RUQygFHAQhFZ7NreUUQWuYqNAW4HLqllmO8zIrJVRLYAFwO/be46cCwNvAMhMKbZP1oppVoCq0dtzQPm1bI9E5joer0KqHVYlzHmdrcGWB+5aRCeoEN/lVLtVovv2mrxju3V6yNKqXZNE0lj2Ksgf79eH1FKtWuaSBojfz84qpx3tSulVDuliaQxjrlGbGnXllKqHdNE0hi5rntItGtLKdWOaSJpjGNp4B0EAU14p7xSSrUymkgaIzcNInTWX6VU+6aJpDGOpen1EaVUu6eJ5FxVVUD+Ab0+opRq9zSRnKv8/WAc2iJRSrV7mkjOlY7YUkopQBPJuau+h0RvRlRKtW+aSM5VURbYfMA/wupIlFLKUppIzlVpLgRE6tBfpVS7p4nkXJVkOxOJUkq1c5pIzlVJDvhrIlFKKU0k56okR1skSimFJpJzV5qjc2wppRTWr9l+g4hsExGHiCSeoVy6a232ZBFJqrE9XESWiMge13NYswReUQKVpTpiSymlsL5FkgJcC6ysR9mLjTGDjTE1E850YJkxpgewzPXe/UpynM/aIlFKKWsTiTFmhzFmVyNOMQmY7Xo9G5jc+KjqofREItFrJEopZXWLpL4M8J2IbBCRaTW2xxhjsgBcz9F1nUBEpolIkogkZWdnNy4abZEopVQ1T3d/gIgsBTrUsutxY8z8ep5mjDEmU0SigSUistMYU5/usGrGmJnATIDExETTkGNPcyKR6DUSpZRyfyIxxlzWBOfIdD0fFZF5wHCc11WOiEisMSZLRGKBo439rHopcbVotEWilFItv2tLRAJEJOjEa+BynBfpARYAU12vpwL1beE0TmkOePqCd0CzfJxSSrVkVg//nSIiGcAoYKGILHZt7ygii1zFYoBVIrIZWAcsNMZ869o3AxgvInuA8a737leS62yN6DxbSinl/q6tMzHGzAPm1bI9E5joer0XGFTH8bnApe6MsVYl2Xp9RCmlXFp811aLVKrToyil1AmaSM5FiU6PopRSJ2giaShjXDP/ateWUkqBJpKGqyiBqjJtkSillIsmkobS6VGUUuokmkgaqiTX+awtEqWUAjSRNNyJu9p1dUSllAI0kTScdm0ppdRJNJE0VPU8W5pIlFIKNJE0XEkOePrpPFtKKeWiiaSh9GZEpZQ6iSaShirNgQC9GVEppU7QRNJQ2iJRSqmTaCJpqJIcHfqrlFI1aCJpCGN05l+llDqFJpKGqCiGquOaSJRSqgZNJA1RcuJmRL1GopRSJ2giaYhS1zxbeo1EKaWqWb1m+w0isk1EHCKSWEeZXiKSXONRKCIPu/Y9JSKHauyb6NaA9a52pZQ6jaVrtgMpwLXAv+sqYIzZBQwGEBEbcIiT13l/wRjzrDuDrFai82wppdSpLE0kxpgdACJS30MuBdKMMfvdFtSZ6My/Sil1mtZ2jeTnwIenbHtARLaIyCwRCavrQBGZJiJJIpKUnZ19bp9emgteAeDtf27HK6VUG+T2RCIiS0UkpZbHpAaexxu4Bvi0xubXge44u76ygOfqOt4YM9MYk2iMSYyKOsdRV5E9of+153asUkq1UW7v2jLGXNZEp7oS2GiMOVLj3NWvReRN4Osm+qzaDZvqfCillKrWmrq2buaUbi0Ria3xdgrOi/dKKaWakdXDf6eISAYwClgoIotd2zuKyKIa5fyB8cAXp5ziGRHZKiJbgIuB3zZT6EoppVysHrU1j5OH8p7YnglMrPG+FDht7nZjzO1uDVAppdRZtaauLaWUUi2QJhKllFKNoolEKaVUo2giUUop1SiaSJRSSjWKGGOsjqHZiUg20JD5uiKBHDeF05JpvduX9lpvaL91b2i9uxpjTpsapF0mkoYSkSRjTK3T3LdlWu/2pb3WG9pv3Zuq3tq1pZRSqlE0kSillGoUTST1M9PqACyi9W5f2mu9of3WvUnqrddIlFJKNYq2SJRSSjWKJhKllFKNoomkBhGZICK7RCRVRKbXsl9E5CXX/i0iMtSKOJtaPep9q6u+W0RktYgMsiLOpna2etcod76I2EXk+uaMz13qU28RGSciySKyTURWNHeM7lCP3/MQEflKRDa76n2nFXE2Ndcy5EdFpNb1mprke80Yow/ndSIbkAZ0A7yBzUDfU8pMBL4BBBgJ/Gh13M1U79FAmOv1le2l3jXK/QdYBFxvddzN9O8dCmwHurjeR1sddzPV+0/A067XUcAxwNvq2Jug7hcBQ4GUOvY3+ntNWyQ/GQ6kGmP2GmMqgI+AU9eVnwS8Z5zWAqGnrNLYGp213saY1caYPNfbtUDnZo7RHerz7w3wIPA5cLQ5g3Oj+tT7FuALY8wBAGNMW6h7feptgCARESAQZyKpat4wm54xZiXOutSl0d9rmkh+0gk4WON9hmtbQ8u0Ng2t0y9x/vXS2p213iLSCecSzm80Y1zuVp9/755AmIgsF5ENInJHs0XnPvWp9ytAHyAT2Ao8ZIxxNE94lmr095qlKyS2MFLLtlPHRtenTGtT7zqJyMU4E8kFbo2oedSn3i8Cjxlj7M4/UtuE+tTbExgGXAr4AWtEZK0xZre7g3Oj+tT7CiAZuAToDiwRke+NMYXuDs5ijf5e00Tykwwgrsb7zjj/MmlomdamXnUSkYHAW8CVxpjcZorNnepT70TgI1cSiQQmikiVMebL5gnRLer7e55jjCkBSkRkJTAIaM2JpD71vhOYYZwXDlJFZB/QG1jXPCFaptHfa9q19ZP1QA8RSRARb+DnwIJTyiwA7nCNchgJFBhjspo70CZ21nqLSBfgC+D2Vv5XaU1nrbcxJsEYE2+MiQc+A+5v5UkE6vd7Ph+4UEQ8RcQfGAHsaOY4m1p96n0AZysMEYkBegF7mzVKazT6e01bJC7GmCoReQBYjHOExyxjzDYRuc+1/w2cI3cmAqlAKc6/YFq1etb7z0AE8Jrrr/Mq08pnSq1nvduc+tTbGLNDRL4FtgAO4C1jTK1DR1uLev57/y/wrohsxdnd85gxptVPLS8iHwLjgEgRyQD+AnhB032v6RQpSimlGkW7tpRSSjWKJhKllFKNoolEKaVUo2giUUop1SiaSJRSSjWKJhKllFKNoolEKaVUo2giUcpCItLDNTlikog8IyKpVsekVENpIlHKIiJiA94DfueaKcAP2GZtVEo1nCYSpawzGdhujNnoer8D2CIi3UTkbRH5zMLYlKo3TSRKWWcIzmnLTxgEbHYtvvRLi2JSqsE0kShlnVyc05QjIiOAO3BOlKhUq6KJRCnrzAESXbPNXoszsejFdtXqaCJRyiLGmBxjzAhjzACcy7weMsY4RCRCRN4AhojIHy0OU6mz0vVIlGoZBuHq1nKtQHmfteEoVX+6HolSSqlG0a4tpZRSjaKJRCmlVKNoIlFKKdUomkiUUko1iiYSpZRSjaKJRCmlVKNoIlFKKdUo/w/CgrzDQQwOGgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure()\n", "plt.plot(q1_list, lam_avg_list, label='simulation')\n", "plt.plot(q1_list, lam_theory_list, label='theory')\n", "plt.xlabel(r'$q_1$')\n", "plt.ylabel(r'$\\Lambda$')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "automotive-glory", "metadata": { "slideshow": { "slide_type": "subslide" }, "tags": [] }, "source": [ "It can be seen that the simulation results follow the theoretical values pretty closely. We will use the theoretical values to find the optimal proportion of persisters. This can be simply \"read off\" from the plot above." ] }, { "cell_type": "code", "execution_count": 11, "id": "banned-reproduction", "metadata": { "slideshow": { "slide_type": "fragment" }, "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "optimal proportion of persisters = 0.75\n", "maximum growth rate = 0.058892\n" ] } ], "source": [ "imax = np.argmax(lam_theory_list)\n", "print(f'optimal proportion of persisters = {q1_list[imax]}')\n", "print(f'maximum growth rate = {lam_theory_list[imax]:.6f}')" ] }, { "cell_type": "markdown", "id": "f14cb87f-ca67-468b-b3cf-f54e95d559d3", "metadata": { "slideshow": { "slide_type": "fragment" }, "tags": [] }, "source": [ "The fact that the optimal proportion is reached at an intermediate value between 0 and 1 means that the cell population diversifies into a mixture of normal and persister cells. That is, a bet-hedging strategy is optimal!" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.10" } }, "nbformat": 4, "nbformat_minor": 5 }