{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Birth-Death Process"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Previously we have simulated the growth of bacterial populations by simply assuming that cells double synchronously, so their number increases geometrically and deterministically as $N = 2^G$, where $G$ is the number of generations. In this note, we will make a more realistic (yet still very idealistic) model of population growth by incorporating the stochastic aspect of the process.\n",
    "\n",
    "In particular, we recognize that the dynamics of the population size consists of discrete events of births and deaths of individuals, which can happen asynchronously. Therefore, we will study the process in continuous time, instead of discrete generations, and let each individual have a certain probability of giving birth to a new individual (i.e., asexual reproduction) or dying at any given time. Note that we specify only the probabilities, not the time of actual events. The process will be *stochastic* in the sense that, if we repeat the simulation, the events can happen at different times. If we plot the population size as a function of time, then the trajectories from those simulations will fluctuate in time and differ from each other.\n",
    "\n",
    "Let us assume that every individual has a probability $\\beta$ per unit time of giving birth to a new individual, and a probability $\\delta$ per unit time of dying. That means, in a short period of time $dt$, the probability of having a birth or death event per individual is $\\beta\\,dt$ and $\\delta\\,dt$, respectively. Here $\\beta$ and $\\delta$ are called the birth rate and the death rate, respectively. The population dynamics that this model describes is called a \"birth-death process\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Deterministic approximation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can represent such a process symbolically by the following \"reactions\":\n",
    "\\begin{alignat}{2}\n",
    "\\mathbb{I} &\\xrightarrow{\\beta} 2 \\, \\mathbb{I} \\,, \\quad &&\\textsf{i.e., one individual becomes two at rate } \\beta \\,; \\\\\n",
    "\\mathbb{I} &\\xrightarrow{\\delta} \\mathbb{0} \\,, \\quad &&\\textsf{i.e., one individual becomes none at rate } \\delta \\,.\n",
    "\\end{alignat}\n",
    "This representation is similar to that for chemical reactions. We can use this analogy to derive the equation that describes the dynamics of the population size, similar to the \"mass action\" equations for chemical reactions. The resulting equation will determine the *average* behavior of the population dynamics.\n",
    "\n",
    "To derive the equation, suppose that at time $t$ the population size is $N(t)$. Then in a short period from $t$ to $t+dt$, the *average* number of births will be $N \\beta \\, dt$, and the *average* number of deaths will be $N \\delta \\, dt$. Therefore the population size will change by $dN = N \\beta \\, dt - N \\delta \\, dt$. This can be written as a differential equation:\n",
    "\\begin{equation}\n",
    "\\frac{dN}{dt} = (\\beta - \\delta) N\n",
    "\\end{equation}\n",
    "The solution is simply:\n",
    "\\begin{equation}\n",
    "N(t) = N_0 \\, \\mathrm{e}^{(\\beta - \\delta) t}\n",
    "\\end{equation}\n",
    "That is, the population is expected to grow exponentially with time, with a net growth rate equal to the difference between the birth and the death rates. This is rather intuitive, and gives qualitatively the same behavior as geometric growth for synchronous generations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an exercise, let us numerically solve the differential equation above and compare to the analytic solution. In choosing values for $\\beta$ and $\\delta$, note that they both have units of 1/\\[time\\], so we can always rescale time to set one of them to 1. We will set the death rate $\\delta = 1$, so that the average lifespan of an individual is $1/\\delta = 1$, i.e., the generation timescale is $1$. In order for the population to grow, we need $\\beta > \\delta$, such as $\\beta = 2$. We will start from $N_0 = 1$, i.e., from one individual (think about a seed that is dispersed to a new patch of soil). To solve the equation, we can discretize time into very small time steps $\\Delta t$ and calculate the population growth $\\Delta N$ within each time step. (This method is simple but not as accurate as some more sophisticated methods, which we will use later in the course.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta = 2.    # birth rate\n",
    "delta = 1.    # death rate\n",
    "\n",
    "dt = 0.01    # time step size\n",
    "T = 10.    # total time to solve the equation for\n",
    "t_array = np.arange(0, T, dt)    # time points to evaluate the variable N\n",
    "\n",
    "N0 = 1    # initial population size\n",
    "N_array = np.zeros(len(t_array))    # initialize values of N at every time point\n",
    "N_array[0] = N0\n",
    "\n",
    "for i in range(1, len(t_array)):\n",
    "    dt = t_array[i] - t_array[i-1]    # this is true even if t_array is not uniformly spaced\n",
    "    dN = N_array[i-1] * (beta - delta) * dt    # increment to first order in dt\n",
    "    N_array[i] = N_array[i-1] + dN    # calculate value at next time point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAuG0lEQVR4nO3de5yM9fvH8ddlWWvXWseIjd0iOZMVFaWkKIRypiNCDqUSpeigVOqXSoqSdCCJEEURJfRFikI5s47rtBZ7mpnr98eu/fqKNdjZe2fmej4ePezcM3Pf79nYa+/7c3+uj6gqxhhjDEA+pwMYY4zJO6woGGOMyWJFwRhjTBYrCsYYY7JYUTDGGJMlv9MBLkbJkiU1JibG6RjGGONXVq1adUBVS53pOb8uCjExMaxcudLpGMYY41dEZPvZnrPLR8YYY7JYUTDGGJPFL4uCiLQUkXGJiYlORzHGmIDil2MKqjobmB0XF9fj9OfS09OJj48nJSXFgWSBKSwsjOjoaAoUKOB0FGOMj/llUchOfHw8kZGRxMTEICJOx/F7qsrBgweJj48nNjbW6TjGGB/zy8tH2UlJSaFEiRJWEHKIiFCiRAk78zImSARcUQCsIOQw+34aEzwCsigYY0zASj0G856Go3t8snsrCgFo1qxZjBw58oLeGxMTw4EDB3I4kTEmR2xagGdMA3TZGI6sneuTQwTcQHOwc7lctGrVilatWjkdxRiTU04cgvlD4ffP2CHleCLtWdoXvJV2PjiUX54p5PV5Ctu2baNKlSr06NGDatWqceutt5KcnEzjxo2z2nIcOHCAk32bJk6cSOvWrWnatCkxMTG88847vPHGG9SpU4cGDRpw6NAhADZv3kyzZs2oW7cujRo1YsOGDQDcd9999OrVi/r16zNo0CAmTpxI3759Adi3bx9t2rShVq1a1KpVi6VLlwLQunVr6tatS7Vq1Rg3blwuf4eMMV5bNxPPO9fg/n0yb7ta0zfyLZ7u/QDt4i7zyeH88kwhu3kKp3pu9l+s2300R49dtWwRhrWsds7Xbdy4kcmTJzN+/Hjat2/PV199le3r//zzT1avXk1KSgoVK1bklVdeYfXq1Tz66KNMmjSJRx55hJ49e/Lee+9RqVIlfv31V/r06cPChQuBjFtxly5dSkhICBMnTszab//+/bnxxhuZMWMGbrebY8eOATBhwgSKFy9OcnIy9erV46677qJEiRIX/o0xxuSspL3o3MeR9bP5m1ieTB9Ik5uaMr3xFYTm993v835ZFPxBbGwstWvXBqBu3bps27Yt29ffdNNNREZGEhkZSVRUFC1btgSgRo0arFmzhmPHjrF06VLatfvvCWNqamrW1+3atSMkJORf+124cCGTJk0CICQkhKioKADeeustZsyYAcDOnTvZuHGjFQVj8gJV+P1zPN8NwZ16gtfTO7KibBdG3V2HK0tH+vzwAV0UvPmN3lcKFiyY9XVISAjJycnkz58fj8cD8K/7/k99fb58+bIe58uXD5fLhcfjoWjRovz+++9nPF5ERITX2RYtWsQPP/zAsmXLCA8Pp3HjxjYPwZi84PB2dPYAZMuPrNarGOp5iPa338TUa2MIyZc7t4b75ZiCv4qJiWHVqlUATJs27bzeW6RIEWJjY/nyyy+BjJnGf/zxxznf16RJE8aOHQuA2+0mMTGRxMREihUrRnh4OBs2bGD58uXn+UmMMTnK44Ff38czpgEpW5czNP1+Rl/2JuMe7cj918fmWkEAKwq56vHHH2fs2LHUqVPngm77/Oyzz/jwww+pVasW1apVY+bMmed8z+jRo/nxxx+pUaMGdevWZd26dTRr1gyXy0WVKlUYPHgwDRo0uJCPY4zJCQl/45lwG3w7iCVplWitr1O77eN8/GADLisenutxRFVz/aA5JS4uTk9fZGf9+vVUqVLFoUSBy76vxuQwdzr8MhrPolc4pgUZltqVtGrtGN6qOqUiC577/RdBRFapatyZngvoMQVjjMmTdv+O5+uHybf/T+a6GzAmrAeP3t2QW6uVcTqZFQVjjMk16cmwaCS69G0OU4QhaY9SIq4tU5pXIapQ3mhNb0XBGGNyw/aluGf2JeTQZr5wNeaTIt15uuu1XHdFSaeT/Y88VRREJAJYDAxX1W+czmOMMRctNQl+GA4rPmAvlzA4bQhVG97JtFuupFDov+cWOc2nRUFEJgAtgP2qWv2U7c2A0UAI8IGqnuze9iQw1ZeZjDEm12z8HvesAUjSbj5yNWdWiQd4od011Iwu6nSys/L1mcJE4B1g0skNIhICjAGaAvHAChGZBZQD1gFhPs5kjDG+deIQ+t1gZM0XbCOawa7nadzkdqbdcDkFQvL2TACfplPVn4BDp22+BtikqltUNQ2YAtwJNAYaAJ2BHiKSt79zueDUxnbZvWb37t1Zj7t37866det8Hc0Ycyaq8Od03G/H4V4zjdGuNjx9yRhe7v8AD99UMc8XBHBmTKEcsPOUx/FAfVXtCyAi9wEHVNVzpjeLSE+gJ0D58uV9m9QPTJw4kerVq1O2bFkAPvjgA4cTGROkju5B5zyG/D2H9Xo5z+iTtL3jNj6vX4F8uTgj+WLlubKlqhOzG2RW1XGqGqeqcaVKlcrNaOflTK2pCxcuzNNPP02tWrVo0KAB+/btA2D27NnUr1+fOnXqcMstt2RtPykpKYnY2FjS09MBOHr0aFbLi5UrV9KlSxdq1679r/bc3333HVdffTW1atWiSZMmufjpjQkiqvDbJNzvXEPa39/zUnon3qjwLu8MvIdu18b4VUEAZ84UdgGnNgKPztzmNRFpCbSsWLFi9i/8djDsXXu++bJXpgY0P/eqZmdqTX38+HEaNGjAiBEjGDRoEOPHj2fo0KE0bNiQ5cuXIyJ88MEHvPrqq7z++utZ+4qMjKRx48bMmTOH1q1bM2XKFNq2bUu7du0YM2YMo0aNIi7ufycnJiQk0KNHD3766SdiY2Oz1mQwxuSgQ1vxzBpAvm2LWempwkshvbn/7lsYUrus365t7kRRWAFUEpFYMopBRzLGEbzm7XoKTjpTa+rQ0FBatGgBZLTT/v7774GMtRA6dOjAnj17SEtLIzY29l/76969O6+++iqtW7fmo48+Yvz48dkef/ny5dxwww1Z+ypevHhOfjxjgpvHDb++j3vB86S4YET6gxyv1oUJrapTorBvW1T4mq9vSZ1MxgBySRGJB4ap6oci0heYR8YtqRNU9a/z3K93Zwpe/EbvC2drTV2gQIGs3x5CQkJwuVwA9OvXj4EDB9KqVSsWLVrE8OHD/7XP66+/nm3btrFo0SLcbjfVq1f/12uMMblg/3rcXz9MyO5VLHLX4c2w3jzSoTFNqpR2OlmO8GlRUNVOZ9k+F7jgVafz+pnC+bamTkxMpFy5cgB8/PHHZ33dPffcQ+fOnXnmmWeytkVGRpKUlPSv1zZo0IA+ffqwdevWrMtHdrZgzEVwpcEvb+JZ/CpJnjCeTXuYIvU68nnzKkSG5Y0WFTkhzw00B4LzbU09fPhw2rVrR926dSlZ8uxT3rt06cLhw4fp1Om/tfbk+swnB5pPKlWqFOPGjaNt27bUqlWLDh06XPwHMyZY7VqF+/0b4ccRfJMex/0RY+jS/TFebFMzoAoC+Gnr7FMuH/XYuHHj/zwXyC2ep02bxsyZM/nkk09y/diB/H015qzSTsCil9Gl75BAUYam30/FRu3p36QSYQXyXosKbwVc6+y8fvnIF/r168e3337L3LkXfNXNGHM+ti3B9XVf8h/ZymTXzcwo+RDD2l1L9XJRTifzKb8sCsHo7bffdjqCMcEhJRH9fhiy6iN2U5qh7me49pY2TG4US34/mJF8sfyyKJzr7iNV9dt7hPMif7zEaMwF+WcerpkDyHd8H+Ndd/BTdA+ev6sel5cq7HSyXOOXZU9VZ6tqz6iof5/GhYWFcfDgQftBlkNUlYMHDxIWZn0KTQA7fgDPtO7weXu2HMtPZ32RiJYvM6ln46AqCOCnZwrZiY6OJj4+noSEBKejBIywsDCio6OdjmFMzlOFP7/CNecJNOUoo9PvYkOl7vxfmzpcGlXI6XSO8MuikN3lowIFCpxxRrAxxvyPo7txz36UkI3f8afnCl7KP5Ru7ZvzSM1Lg/rys18WhWC8+8gYk0NU4bePcX83lPT0VF5K78qRGg/wfssaFIsIdTqd4/yyKBhjzAU5uBn3rP6EbF/Cr56qvBnWlz6dmtK48iVOJ8szrCgYYwKfxw3L38W94EWS3fl4Mb0HBevdy4TmVShc0H4Mnsq+G8aYwLZvHa4Zfci/dzUL3VfzQZG+DGp/E3UrWC+wM/HLouB1l1RjTPBypaE/j0J/eoOjWojhrv5UaNSFj2/27xYVvuaXRcEGmo0x2YpfRfqM3hQ4+Ddfu69nWqmHGdquEVUuLeJ0sjzPL4uCMcacUdpxdOEIdPlYDmoxhnkGEde0Ex9fHxMULSpyghUFY0xg2LKY9K/7UeDodj513cKP0X0Y1q4BFUpEOJ3Mr1hRMMb4t+QjeOY/Q77Vk9ilZXhOhtPszrv4MO6yoJ6EdqH8sijYQLMxBoANc0mf9Qj5TiTwnqslayv2ZmTbupQuYr26LpRfXmTLriGeMSYIHEvAPfU+mNKJzcdCuS/kZcp3eI137r3OCsJF8sszBWNMkFKFtV+SPmcQmprE/6W3I6FmL95uWZOi4daiIidYUTDG+IfEeFyzBpB/8w+s9VTkjULP8VCX22lUqZTTyQKKFQVjTN7m8cCqCbjmPUu6y8UIVzfyXfMQ45pVITzUfoTlNPuOGmPyroObSZ/xMAXil7HMXZ3xRfvzSLtbubp8MaeTBSwrCsaYvMftQpe9g2fhS6R4QnjW1ZPSN3bng5sqEZrfL++P8RtWFIwxecvetaRN70Po/jX84I5jSqkBDOlwE1eWjnQ6WVCwomCMyRtcqejiV9Elb3LUE8ELnkepees9fHB9LCH5bBJabvHLomCT14wJMDv/k3F2cHgj092N+D56AE+3u57Lioc7nSzo+OXFOZu8ZkyASD2Ge+6T6Ie3knDoML15Cs+dYxnb8xYrCA7xyzMFY0wA2LyQtBn9CT22k49dTfmj8iM81zqOSyJtRrKTrCgYY3JX8mFc3z1N/j8+I14vZWT+F7irXXvurVbG6WQGKwrGmNy0fjZpsx4lJPkgY1yt2FOrP6+1qE1UoQJOJzOZrCgYY3zv2H7Sv3mcAhtmstFTgf8LH8wDd7fm4YolnU5mTmNFwRjjO6rwxxTS5z6JJ+0Er7k64G7Qj7dvrUqhUFsnOS+yomCM8Y0jO0j7uj+h237kD8+VvBf1CP3a30Gty4o6ncxkw4qCMSZneTzoig9wzR+Gy+1mpPs+it3Yh3cbW4sKf2BFwRiTcw5sJHV6Hwru/g/L3DX4vPRjPN7+FipeYi0q/EWeKQoiUgUYAJQEFqjqWIcjGWO85U7H88vb6KKXSXEX4DntQ+VmPXn32hjyWYsKv+LToiAiE4AWwH5VrX7K9mbAaCAE+EBVR6rqeqCXiOQDJgFWFIzxB3v+IOWrPoQd+JM57muYV/4xBt19A9HFbEayP/L1Bb6JQLNTN4hICDAGaA5UBTqJSNXM51oBc4C5Ps5ljLlY6Sm45w/H8/5NJCXsZKA8Tlrbjxjd/TYrCH7Mp2cKqvqTiMSctvkaYJOqbgEQkSnAncA6VZ0FzBKROcDnZ9qniPQEegKUL1/eV9GNMdnZsTzj7CBxM1NdN7Lyqsd5qnV9ShYu6HQyc5GcGFMoB+w85XE8UF9EGgNtgYJkc6agquOAcQBxcXHqs5TGmH9LTSJ9/nDyr/qQA1qCVws8S6v2XXm1ammnk5kckmcGmlV1EbDIm9da62xjHLDpB1Jm9CP0+B4mum5lR+3HeLHF1RQJsxYVgeScYwoicqWILBCRPzMf1xSRoRdxzF3AZac8js7c5jVrnW1MLjpxiLRpD8GndxGfBP0LvUyVB8Yy7O76VhACkDdnCuOBJ4D3AVR1jYh8Drx4gcdcAVQSkVgyikFHoPMF7ssY40vrZpI6cyAhqYcZ427NiQYDGXVrdcIKWIuKQOVNUQhX1f+I/M+9xi5vdi4ik4HGQEkRiQeGqeqHItIXmEfGLakTVPWv8wltl4+M8bGkvaTMGkjYxjn844nh/ahn6dWxNdXL2dl5oPOmKBwQkSsABRCRu4E93uxcVTudZftcLuK2U1WdDcyOi4vrcaH7MMacgSq6+lPSv30KSU/mNXcnCt/0CP9345UUCLEWFcHAm6LwMBl3+1wlIruArUBXn6Y6BztTMMYHDm8neXpfCu38idWeq/j8ksfp36E5V5Qq7HQyk4tE1bu7OkUkAsinqkm+jeS9uLg4XblypdMxjPFvHjeeX8fh/uE50lzKG3Qhtlk/Ote3FhWBSkRWqWrcmZ4755mCiLiB14AhmllBROQ3Vb06Z2MaY3Jdwt8kT+tDoX0rWeyuxdwKg3j07iaULVrI6WTGId5cPvqLjFtX54tIB1U9BDj664NdPjLmIrnTcf38f7D4VVI8BRmRrz/X3NWLV2uV5bSbSkyQ8aYouFR1kIh0AH4WkXvIHHR2ig00G3MRdq/mxJe9CD+8gW/cDfj1qsEMbH09xSNCnU5m8gBvioIAqOoXIvIXGT2JrOmQMf4mPZm0BS+Rf/kYkjSS5ws8yW2duvNC5UucTmbyEG+KQveTX6jqnyLSiIwGdsYYf7HtF0581YfwpG1McTVmS53BDG1Rj8IF80ynG5NHnPVvhIjcrKoLgQoiUuG0p4/5Nlb2bEzBGC+lHCV13rMUXP0RBzylGB3xPB07dKNjTHGnk5k8KrtfE24EFgItz/CcAtN9ksgLNqZgzLnpP/NJmdGfgsl7meC+naRrn2RE0xrWosJk66xFQVWHZf55f+7FMcZctOMHSf5mEIXWT2Onpxzjio7igY7tqVq2iNPJjB/wZp7CAOAjIImM5nhXA4NVdb6Psxljzocq+tcMUmc9Rv60RN7x3EXYzU8wslFl8luLCuMlb0aZHlDV0SJyG1AC6AZ8AjhWFGxMwZjTHN3DiRkDCN86j789l/NZ6Rfp06E1MSUjnE5m/IzXt6QCtwOTVPUvcXh2i40pGJNJFfeqSbi+e4p86amMoivRtz/GK/VjbRKauSDeFIVVIjIfiAWGiEgk4PFtLGPMOR3ayvFpDxOx+xdWeKowq8IQ+t99G2WiwpxOZvyYN0XhQaA2sEVVT4hICcAGn41xisdN+tJ30YUvom5hRMhD1G4zgBE1rUWFuXjnLAqq6gF+O+XxQeCgL0MZY85i/3qOTe1F4QO/s8Bdh1+uepp+rW+kmLWoMDnEpjMa4w9caaQufp2QJaNI9RRiVOhAbu7Um2etRYXJYX5ZFOzuIxNUdq3i2NTeFE78m5nu6/inztM8cUcDIqxFhfEBrxbZEZEQoDSnFBFV3eHDXF6xRXZMQEs7QfL3L1JwxVj2aVHeDe9D647dqVuhmNPJjJ+72EV2+gHDgH38964jBWrmWEJjzP/QrT9xfNrDFD6+g8nuJhy6bihDm9aiYH5rUWF8y5vzzwFA5cwBZmOML6UkcnzOUCLWTuKApzQji71Mt07dqFwm0ulkJkh4UxR2Aom+DmJMsPNs+JbkGf0JS0ngQ21JvpuH8NwNVQmxdZJNLvKmKGwBFonIHCD15EZVfcNnqYwJJscPcOzrxyi88Wt2ei7j8zKj6d6hHeVLhDudzAQhb4rCjsz/QjP/M8bkBFVca74k/ZsnCE1LYoy0p3SLITxXz1pUGOd4M3ntOQARKZz52NEFdowJCIm7OPpVP4rsWMCfniuYXeF1HmrfgksirUWFcZY3dx9VJ6MravHMxweAe1T1Lx9nyy6TzVMw/snjIX3FR7jnP0MBVzpvhNxH1bue4Jka0U4nMwbw7vLROGCgqv4IICKNyVhX4TrfxcqedUk1fungZo5O7UORfcv5j7saS6o8S6/WTYgKL+B0MmOyeFMUIk4WBABVXSQi1qTdGG+5XaQseZuQRS+BJ4RXCj5Mw3aP8mSlUk4nM+ZfvLr7SESeIeMSEkBXMu5IMsacy76/SPziIaIOreV7d13W1XmWfi0aEh5qLSpM3uTVymvAc8D0zMc/Z24zxpyNK5UTC16h4LI3SdcIXooYxO0d+zCgvLWoMHmbN3cfHQb650IWYwKC7lzB0am9iEraxExPQ/ZdN5zHb6lLaH5bJ9nkfWctCiLypqo+IiKzyeh19D9UtZVPkxnjb9KOk/TtcCJWj+e4FuPtos/RofOD3FnaWlQY/5HdmcLJMYRRuRHEGH/m2bSI41/1ITJ5F5O1KZ6bhzGkUXVrUWH8zlmLgqquyvyytqqOPvU5ERkALPZlMGP8QvIREmcNJmr9ZBI8ZXirzOvc26kz0cWsRYXxT94MNN8LjD5t231n2GZMUHGt+4aUrwdQOPUgH8mdFGv5LE/FXWEtKoxfy25MoRPQGYgVkVmnPBUJHPJ1MGPyrGMJHJn+KEW3zGanpzyzY17k/nZtKRVZ0Olkxly07M4UlgJ7gJLA66dsTwLW+CKMiLQG7gCKAB+q6nxfHMeYC6JK2uopuOYMopDrBGPzd6LS3UMZVN1aVJjAkd2YwnZgO3DtxRxARCYALYD9qlr9lO3NyLgEFQJ8oKojVfVr4GsRKUbGALcVBZM3JMZz+IuHKbZ7EWs9lfipyjAebNOMImHWosIElnPeOC0iDURkhYgcE5E0EXGLyNHzOMZEoNlp+wwBxgDNgapAJxGpespLhmY+b4yzPB6Sl75P6ug4Cu5axluh3Um/Zy6PdmppBcEEJG8Gmt8BOgJfAnHAPcCV3h5AVX8SkZjTNl8DbFLVLQAiMgW4U0TWAyOBb1X1tzPtT0R6Aj0Bypcv720MY87fgU0cmvIQxQ+sZImnOmvqPE/PFo0JK2DrJJvA5dUUS1XdBISoqltVP+K03/wvQDkylvk8KT5zWz/gFuBuEel1lizjVDVOVeNKlbKGYsYH3C6OLRhF2phrCUlYxxvhA4jqMYc+bZpYQTABz5szhRMiEgr8LiKvkjH47JP5+qr6FvDWuV5n6ykYX9E9azgy5SGKJa5jvqceu657gX5N61MgxFpUmODgzd/0bmQMBvcFjgOXAXdd5HF3Ze7npOjMbV5R1dmq2jMqKuoiYxiTyZVK4pxhuN9vjOtIPK8XfYor+s7g/mbXWkEwQcWbhnjbM79MJqNbak5YAVQSkVgyikFHMuZEeMXOFExOcm9fztGpvSl2fAtf642k3vwCjzaqST5rUWGCUHaT19ZyhkZ4J6lqTW8OICKTgcZASRGJB4ap6oci0heYR8ZZyITzWd7TVl4zOSL1GIe/eYaotR9xXEswoczLdOp0P2WLFnI6mTGOye5MoUVOHEBVO51l+1xgbk4cw5jzlf7PAk581ZdiqbuZIs2JvON5BsZVshYVJuida/JanmSXj8wFSz7MwelPUGLjlxzwXMrEmLfo2r4jJQpbiwpjAET1rFeIMl4gksR/LyOFAgWA46paxMfZzikuLk5XrlzpdAzjJ1LXfk3azIEUSj/MZ/lbU6HtczSuZnNdTPARkVWqGnem57wZaM5aIUQyzq3vBBrkXDxjfCxpHwe+HEDJHd+yyVOBxVeNolvblkTajGRj/uW87rXTDF8Dt/kmjndEpKWIjEtMTHQyhsnrVDnxn0848WZdIrf/wPjQbpy493v6dG5rBcGYszjnmYKItD3lYT4yWl2k+CyRF+zuI3NOR3aQMLk3pfYtYaWnMn9c/Tzd7mhqM5KNOQdvZjS3POVrF7CNjEtIxuQ9Hg9Hfx5L6KLnCfco70U8RKMuQ3iwXDGnkxnjF7wZU7g/N4IYc7E04R8OfN6TUodX87OnFjuuH0H3W64nv81INsZr3rTOvlxEZotIgojsF5GZInJ5boTLJpONKZj/cqdzeN5I0sdcR4FD/zAm6nGi+82ly22NrCAYc568uSV1ORlrG0zO3NQR6Keq9X2c7ZzsllTjiv+dI1N6UvLY38zXBhxr8jKtG9axFhXGZOOibkkFwlX1k1MefyoiT+RMNGMuUHoKB+Y8T7Hfx+LRIrxbejhtu/SmTFSY08mM8WveFIVvRWQwMIWMSWwdgLkiUhxAVQ/5MJ8x/5K25ReOTe1FyZQdzJSbCGvxMr3rXmUtKozJAd4UhfaZfz502vaOZBSJXB9fsDYXQSo1if0zhnDJhk847inF1JhRdOhwL8UiQp1OZkzAOOeYQl5mYwrBI3n9fFKn96VI2n6m5b+dS9u+RKNqMU7HMsYvXdSYgogUAHoDN2RuWgS8r6rpOZbQmLM5cYi9UwdSZtsMdnnK8tVV79HxrruJKOjNSa4x5nx58y9rLBlN8N7NfNwtc1t3X4UyBlWOrf4KnfM4JVxH+bRge6p2fIEHLy/jdDJjApo3RaGeqtY65fFCEfnDV4GM0aN72Du5L5fu+YG1nlhW13mbDi1vp2B+a1FhjK95UxTcInKFqm6GjMlsgNu3sUxQUuXIso8o8MMzFHOnMjHifup3eZZ7yhV3OpkxQcObovAE8KOIbAEEqAA42vrC7j4KPJ6DW9n3eS8uPbicFXoVW68bSbemjQmxSWjG5Cqv7j4SkYJA5cyHf6tqqk9TecnuPgoAHjcHf3yHiCUjcHmEqUW706TbYCqUjDz3e40xF+Ri7z4KA/oADcmYl/CziLynqo62zzb+z7V3PQc+70mZo2v4mTokNnmV+xvVs0loxjjIm8tHk4Ak4O3Mx52BT4B2vgplApw7nX3fjqT4yjcJ1TAmlH6KFl36c0lUIaeTGRP0vCkK1VW16imPfxSRdb4KZAJb6o5VJE7pSekTm5gv1xNyx6s8UK+607GMMZm8KQq/iUgDVV0OICL1AbuQb85PejK7Zw6j9J/jUS3CxAojaNPxIaLCbVlMY/ISb4pCXWCpiOzIfFwe+FtE1pKxbHNNn6UzAeHEP4s58VUfyqbGMzvkFkq2fY37qjm6JIcx5iy8KQrNfJ7CBKaUo8RPe5LoTZ+ToJfwbeW3uOvuLoSHWosKY/Iqb5bj3J4bQc6HzVPI+46umYNr1iNcmp7AVwVbU7HDS3S7opzTsYwx5+CXaxWq6mxV7RkVFeV0FHMaPX6AnR90pcj0zhxML8CXtT+i5aCPqGUFwRi/YOfxJmeocug/U8g/bzBl3ElMiehM3a4v0LFsSaeTGWPOgxUFc9E8ibvZ9WlvLktYxFq9nM3XjqfdrU2tRYUxfsiKgrlwquxfPJ6IxcMp5Unj86I9adTtGVqXLOJ0MmPMBbKiYC5I+oEt7Pu0J9FHVrCCahxo8hqdGl1nLSqM8XNWFMz58bjZPe9Niv/6CkU0H5+XHkjTroOoV8RaVBgTCKwoGK+lxK/l0OSelD2+jp+lLp473qBzvdpOxzLG5CArCubcXGnsnPUiZdaMoaAWYkr5YTTv1Jeo8FCnkxljcpgVBZOtY1t+5djUXlyWsoXvQ24gqu0oOlarfO43GmP8Up4pCpnLfD4NRKnq3U7nCXppJ9g27Wku+2ciSVqUL698jRbtHqRQqK2TbEwg8+mMZhGZICL7ReTP07Y3E5G/RWSTiAwGUNUtqvqgL/MY7xxZt5CE1+oS888Evgu9lYP3LqZdl55WEIwJAr5uczGR0xrqiUgIMAZoDlQFOolI1X+/1eQ2TT7ClgndKTq1DSfS3Myo+R63PjmZ6peXdzqaMSaX+PTykar+JCIxp22+BtikqlsARGQKcCfg1cI9ItIT6AlQvrz9sMopCSu/JuTbgVRwHWJmxF1U7/IKbcqVcjqWMSaXOdEQrxyw85TH8UA5ESkhIu8BdURkyNnerKrjVDVOVeNKlbIfWhfLnZTA5vc6UOqbe0lwRfDdtZ/S8vEPucIKgjFBKc8MNKvqQaCXN6+11tk5QJW9Sz4hfOHTXOY5zvSi91D/nhe5o4R1njUmmDlRFHYBl53yODpzm9dUdTYwOy4urkdOBgsWaYd2Ev9Jby4//DNrqETCzaNoc8ON1qLCGONIUVgBVBKRWDKKQUegswM5go/Hw84FYyn+ywuUUQ/TL3mYG7sNpWaRcKeTGWPyCJ8WBRGZDDQGSopIPDBMVT8Ukb7APCAEmKCqf53nfu3y0XlK3vsP+z57iJik31ghNUm74w3a1qvndCxjTB4jqup0hgsWFxenK1eudDpG3uZ2sXXOa5T97Q1SNT8Lyvfnls6PEVnIWlQYE6xEZJWqxp3puTwz0Hw+7EzBO0nbfufIFw8Rm7yBJSH1CW/zJm2q25QQY8zZ2RrNgciVyqYvhlBo4s2En9jNzEojiHtyLldbQTDGnINfnimYszu0YQkpX/WhYvp2FoTeRNkOb3LnFTFOxzLG+Am/LAp2+ejfNPUY/0wZQqWtn7BXizOn5tvc1ror+UP88mTQGOMQG2gOAPt+n4fM7s8l7r3MC2/BlV1GEVvuUqdjGWPyqIAbaDYZ3CcOs+nTR6i8+2u26aV832ACTW9rS758NgnNGHNhrCj4qfhlXxL+/SCucB9hblQH6tw7kqYlijsdyxjj5/yyKATzmELqkT1s+6QvlQ/+wN/EsPamcTS/8RZrUWGMyRF+OQoZlLekqrJ14QRSR9cj5sAi5pbqTqmBS7mxcVMrCMaYHOOXZwrB5sT+bez6tBeVji5jrVQm+Y43uf2a65yOZYwJQFYU8jKPh41zR1Nu5UjKqjL3ske4octTFC5U0OlkxpgA5ZdFIRjGFBJ3rufA5J5UOrGGlSG1KdjmbW6vXtPpWMaYAGdjCnmMutNZ/+VzhH3YiJLHN/NdxWep/uQCalhBMMbkAr88UwhUBzat4vjUh6iStpFloddSsv3bNKtYyelYxpggYkUhD/CkJbPui2eovGkCSmG+r/EaN7Xpbi0qjDG5zoqCw/asXYRnZl+qu3ayOLwpl3d9k6blop2OZYwJUlYUHOJKPsr6T5+gWvwX7JUSLL7mPW5o3tHmHBhjHOWXRcHf7z7a/p/ZhH03kGruBBZF3Un1e17nxpIlnY5ljDF291FuSjl6kLVjOlNhbleSPfn5tfGn3PToRC6xgmCMySP88kzBH21a/DnFFj1FFU8iC0p1pe49LxNTpIjTsYwx5n9YUfCxYwfi2f7Jw1RLXMTfcjnbmk+kSYPGTscyxpgzsqLgK6qs/+59yv36AhU1lQXRvWnQZRgR4YWcTmaMMWdlRcEHjuzexN7PelHl+ArWhlQlX+u3aVLjjIscGWNMnmJFIQepx83aGaOouPYNohUWXvEEDTsNJrSAfZuNMf7BflrlkP1b1pD4RS9qpv7FbwXqUrTDGG6uWMXpWMYYc178sijkpXkKnvQ0/pj6PNX+GUsoYSyu9gIN7+pLiLWoMMb4IVFVpzNcsLi4OF25cqVjx9+5bhmu6Q8T69rMr4VuILrLO5SLruBYHmOM8YaIrFLVMw50+uWZgtPSU46z5rOnqLVjEoelCL/UfZPrWtxnLSqMMX7PisJ52rLqe0LnDKCuZxdLo5pTqdubXF+qjNOxjDEmR1hR8FLKsSP89clj1N03jd1cwopGH3Fdk7ZOxzLGmBxlRcEL63+eTrGFg6jjOcCSku2occ9r1Isq5nQsY4zJcVYUspF0eB8bJ/Xn6sPfsU2i+bPZVBpee6vTsYwxxmesKJyJKmvmf0z0smepocdYUu4B6nYdQUx4uNPJjDHGp6wonObg3h3Ef9KbWseX8E9IRRJaTqFh7eucjmWMMbnCikIm9Xj4bdY7XPn7y1TWdJZc3p9rOj1DaGio09GMMSbXWFEA9mzbwMEpvamb8ht/FahO4Xbv0vDKWk7HMsaYXJdnioKIRADvAmnAIlX9zNfH9LhcrPhyJDU2vEUk+VhW9Wnq3/0Y+UJCfH1oY4zJk3zaoEdEJojIfhH587TtzUTkbxHZJCKDMze3Baapag+glS9zAWzfsIp/XmlI/b9f459CtUh68Geu7TDICoIxJqj5umvbRKDZqRtEJAQYAzQHqgKdRKQqEA3szHyZ25ehfv1yFJdOvpXS6fGsuPoVag2ax6XlK/nykMYY4xd8WhRU9Sfg0GmbrwE2qeoWVU0DpgB3AvFkFIZsc4lITxFZKSIrExISLihXxKVXsrZIIzy9l1OvVS8kn3U0NcYYcGZMoRz/PSOAjGJQH3gLeEdE7gBmn+3NqjoOGAcZXVIvJED1hq2goc+vUBljjN/JMwPNqnocuN+b1+al9RSMMSaQOHHdZBdw2SmPozO3eU1VZ6tqz6ioqBwNZowxwc6JorACqCQisSISCnQEZp3PDkSkpYiMS0xM9ElAY4wJVr6+JXUysAyoLCLxIvKgqrqAvsA8YD0wVVX/Op/92pmCMcb4hk/HFFS101m2zwXm+vLYxhhjzp9f3otpl4+MMcY3/LIo2OUjY4zxDb8sCsYYY3xDVC9o/leeICIJwPYLfHtJ4EAOxvEH9pmDg33m4HAxn7mCqpY60xN+XRQuhoisVNU4p3PkJvvMwcE+c3Dw1We2y0fGGGOyWFEwxhiTJZiLwjinAzjAPnNwsM8cHHzymYN2TMEYY8y/BfOZgjHGmNNYUTDGGJMlKIvCWdaIDlgicpmI/Cgi60TkLxEZ4HSm3CAiISKyWkS+cTpLbhGRoiIyTUQ2iMh6EbnW6Uy+JCKPZv6d/lNEJotImNOZfOFM692LSHER+V5ENmb+WSwnjhV0RSGbNaIDmQt4TFWrAg2Ah4PgMwMMIKMTbzAZDXynqlcBtQjgzy8i5YD+QJyqVgdCyGjFH4gmctp698BgYIGqVgIWZD6+aEFXFDj7GtEBS1X3qOpvmV8nkfGDopyzqXxLRKKBO4APnM6SW0QkCrgB+BBAVdNU9YijoXwvP1BIRPID4cBuh/P4xFnWu78T+Djz64+B1jlxrGAsCmdaIzqgf0CeSkRigDrArw5H8bU3gUGAx+EcuSkWSAA+yrxs9oGIRDgdyldUdRcwCtgB7AESVXW+s6lyVWlV3ZP59V6gdE7sNBiLQtASkcLAV8AjqnrU6Ty+IiItgP2qusrpLLksP3A1MFZV6wDHyaFLCnlR5jX0O8kohmWBCBHp6mwqZ2jG3IIcmV8QjEXhoteI9kciUoCMgvCZqk53Oo+PXQ+0EpFtZFwevFlEPnU2Uq6IB+JV9eRZ4DQyikSgugXYqqoJqpoOTAeuczhTbtonIpcCZP65Pyd2GoxF4aLXiPY3IiJkXGder6pvOJ3H11R1iKpGq2oMGf9/F6pqwP8Gqap7gZ0iUjlzUxNgnYORfG0H0EBEwjP/jjchgAfWz2AWcG/m1/cCM3Nipz5djjMvUlWXiJxcIzoEmHC+a0T7oeuBbsBaEfk9c9tTmcuimsDSD/gs8xeeLcD9DufxGVX9VUSmAb+RcYfdagK03UXmeveNgZIiEg8MA0YCU0XkQTKWEGifI8eyNhfGGGNOCsbLR8YYY87CioIxxpgsVhSMMcZksaJgjDEmixUFY4wxWawoGOOlzA6kfTK/Lpt5O6QxAcVuSTXGS5l9o77J7MhpTEAKuslrxlyEkcAVmRMANwJVVLW6iNxHRofKCKASGU3aQsmYMJgK3K6qh0TkCjLatpcCTgA9VHVDbn8IY7Jjl4+M8d5gYLOq1gaeOO256kBboB4wAjiR2ZRuGXBP5mvGAf1UtS7wOPBuboQ25nzYmYIxOePHzLUqkkQkEZiduX0tUDOzQ+11wJcZbXoAKJj7MY3JnhUFY3JG6ilfe0557CHj31k+4EjmWYYxeZZdPjLGe0lA5IW8MXP9iq0i0g4yOteKSK2cDGdMTrCiYIyXVPUg8Evm4umvXcAuugAPisgfwF8E+DKwxj/ZLanGGGOy2JmCMcaYLFYUjDHGZLGiYIwxJosVBWOMMVmsKBhjjMliRcEYY0wWKwrGGGOy/D/Kfu5BG1ZkcwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N_analytic = N0 * np.exp((beta - delta) * t_array)    # analytic solution\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(t_array, N_array, label='numerical')\n",
    "plt.plot(t_array, N_analytic, label='analytic')\n",
    "plt.yscale('log')    # plot in log scale\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('population size')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that, in this solution, the population size is treated as a continuous variable, not an integer. This is because the solution is meant to describe the *average* trajectory over different realizations of the stochastic process. When the size $N$ is large, the deviation of an actual realization of the process from the average trajectory will be small. In that case the deterministic solution can be considered as a good approximation of the true stochastic process. However, when $N$ is small, fluctuations will be significant, and we will have to simulate the stochastic process more seriously."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stochastic simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why should we care about the fluctuations in the stochastic process? Well, the fluctuations can have serious consequences. In the case of population growth, when $N$ is small, there is a real chance that all individuals happen to die before any of them gets to reproduce. Then the whole population goes extinct and will never recover. Such a scenario cannot be captured by the deterministic solution above (because $N$ was treated as a continuous variable and so would never reach zero exactly, even if it becomes $\\ll 1$.) Therefore, to calculate things like the probability of the population going extinct, we need to treat the process stochastically.\n",
    "\n",
    "One method for simulating stochastic processes like those described by the reactions above is the **Gillespie algorithm**. We will explain and use it in the following."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider first the birth process alone. Suppose at the moment the population size is $N$, then the *total* birth rate of the population is $\\beta N$, meaning that the probability of having *one* birth event among *all* individuals in the next short time period $dt$ is $\\beta N dt$. This is similar to the Poisson process that we learned before, where there is a rate $k$ that certain events happen. Recall that the waiting time $\\tau$ before the next event happens is exponentially distributed, $P(\\tau) = k \\, \\mathrm{e}^{-k \\tau}$. Therefore, to simulate the occurrence of the next birth event, we just need to draw a random number $\\tau$ from the exponential distribution with $k = \\beta N$.\n",
    "\n",
    "Now let us also include the death process. Like above, the total death rate among the population is $\\delta N$. Therefore the time before the next death event is exponentially distributed with $k = \\delta N$. We can draw a random number from this distribution to be the time at which the next death would happen. However, we need to account for both birth and death events. There are several ways to do that. One intuitive way (the so-called \"first-reaction\" method of the Gillespie algorithm) is to first generate two random numbers, corresponding to the *putative* times for the next birth and death event respectively. We then compare these numbers to see which event would happen first, and execute only this first event. That is because, after this event, the population size will change (plus or minus one depending on whether the event is a birth or death), and hence the total birth and death rates among the population will change accordingly. So we will have to start over by calculting the new total birth and death rates, drawing two new random numbers as the putative times for the next birth or death event, etc. In this way, the simulation proceeds one event at a time. Note that, in such a simulation, the population size $N$ is always an integer and time is a continuous variable, as compared to the deterministic solution above where $N$ is continuous and time is discretized."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us write a Python class to simulate the birth-death process using the Gillespie algorithm as described above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BirthDeath:\n",
    "    \"\"\"\n",
    "    simulate the birth-death process using Gillespie algorithm.\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, birth_rate, death_rate=1., N0=1):\n",
    "        \"\"\"\n",
    "        initialize the population.\n",
    "        inputs:\n",
    "        birth_rate: float, birth rate per individual.\n",
    "        death_rate: float, death rate per individual.\n",
    "        N0: int, initial population size.\n",
    "        \"\"\"\n",
    "        self.birth_rate = birth_rate\n",
    "        self.death_rate = death_rate\n",
    "        self.pop_size = N0    # current population size\n",
    "        self.time = 0.    # time since beginning of simulation\n",
    "        self.pop_hist = [N0]    # list to record history of population size\n",
    "        self.time_hist = [0.]    # list to record time of all events\n",
    "    \n",
    "    def next_event(self):\n",
    "        \"\"\"\n",
    "        generate the waiting time and identity of the next event.\n",
    "        outputs:\n",
    "        tau: float, waiting time before next event.\n",
    "        event: int, 0 means birth and 1 means death.\n",
    "        \"\"\"\n",
    "        k_b = self.pop_size * self.birth_rate    # total birth rate\n",
    "        k_d = self.pop_size * self.death_rate    # total death rate\n",
    "        t_b = np.random.exponential(1/k_b)    # draw a random number from exponential dist as putative birth time\n",
    "        t_d = np.random.exponential(1/k_d)    # draw a random number from exponential dist as putative death time\n",
    "        if t_b < t_d:    # birth happens first\n",
    "            event = 0    # use 0 to label birth\n",
    "            return t_b, event\n",
    "        else:    # death happens first\n",
    "            event = 1    # use 1 to label death\n",
    "            return t_d, event\n",
    "    \n",
    "    def run(self, T):\n",
    "        \"\"\"\n",
    "        run simulation until time T since the beginning.\n",
    "        inputs:\n",
    "        T: float, time since the beginning of the simulation.\n",
    "        \"\"\"\n",
    "        while self.time < T:\n",
    "            if self.pop_size == 0:    # population is extinct\n",
    "                break    # exit while loop to end simulation\n",
    "            tau, event = self.next_event()    # draw next event\n",
    "            self.time += tau    # update time\n",
    "            if event == 0:    # birth happens\n",
    "                self.pop_size += 1    # increase population size by 1\n",
    "            elif event == 1:    # death happens\n",
    "                self.pop_size -= 1    # decrease population size by 1\n",
    "            self.time_hist.append(self.time)    # record time of event\n",
    "            self.pop_hist.append(self.pop_size)    # record population size after event"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let us test this class by running the simulation several times and plotting these realizations of the stochastic process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "current time = 10.000008613756664, current population size = 58530\n",
      "current time = 10.00000353379543, current population size = 28281\n",
      "current time = 10.000004011786409, current population size = 21343\n",
      "current time = 10.000001743706804, current population size = 70183\n",
      "current time = 0.3460562000994244, current population size = 0\n",
      "current time = 2.8762221984809395, current population size = 0\n",
      "current time = 10.000002400358154, current population size = 58945\n",
      "current time = 0.6621550854798959, current population size = 0\n",
      "current time = 1.6860558927602003, current population size = 0\n",
      "current time = 0.29505707715923196, current population size = 0\n"
     ]
    }
   ],
   "source": [
    "beta = 2.    # birth rate\n",
    "delta = 1.    # death rate\n",
    "\n",
    "T = 10.    # total time for running each simulation\n",
    "trials = 10    # number of simulations to repeat\n",
    "bd_list = []    # list to save all simulations\n",
    "\n",
    "for i in range(trials):\n",
    "    bd1 = BirthDeath(beta, delta)    # create a simulation\n",
    "    bd1.run(T)    # run simulation until time T\n",
    "    print(f'current time = {bd1.time}, current population size = {bd1.pop_size}')\n",
    "    bd_list.append(bd1)    # save the simulation in a list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEJCAYAAAB7UTvrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABf8ElEQVR4nO3dd3gURR/A8e/cXXolIQ0ChN57AOldUZoNBWwowouFUATpUkSlKAKiAoqIKKggRRAUBREE6b33EkhIJb1cmfePDZeEeiEJKczneXy8md3bmQ3J/W6nCikliqIoigKgK+gKKIqiKIWHCgqKoiiKlQoKiqIoipUKCoqiKIqVCgqKoiiKlQoKiqIoipUKCoqiKIqVCgqKoiiKlaGgK3CDEKIN8D5wFPhRSrn5Xu8pWbKkDAoKytd6KYqiFDd79+6NklL63O5YvgYFIcQ3QBcgQkpZK0t+J2AWoAe+llJOASSQCDgCobZcPygoiD179uR5vRVFUYozIcTFOx3L7+ajb4FON1VGD3wOPA7UAHoJIWoAW6WUjwMjgIn5XC9FURTlNvI1KEgptwAxN2U3Bs5IKc9JKdOBH4HuUkpLxvFYwCE/66UoiqLcXkH0KZQGLmdJhwJNhBBPA48BnsCcO71ZCNEf6A9QtmzZ/KuloijKQ6jQdDRLKVcAK2w4bz4wHyA4OPiWJV6NRiOhoaGkpqbmfSWV++Lo6EhgYCB2dnYFXRVFUe6hIILCFaBMlnRgRp7NhBBdga6VKlW65VhoaChubm4EBQUhhMhVRZXck1ISHR1NaGgo5cuXL+jqKIpyDwUxT2E3UFkIUV4IYQ/0BH7NyQWklGuklP09PDxuOZaamoq3t7cKCIWEEAJvb2/15KYoRUS+BgUhxFLgP6CqECJUCNFXSmkC3gb+AI4DP0spj+bwul2FEPPj4uLudDyXNVfykvr3UJSiI79HH/WSUgZIKe2klIFSygUZ+euklFWklBWllB/cx3Xv+KRQGE2YMIGPP/74jsdXrVrFsWPH8q38J554guvXr9/1nPfee4+//vrrjsdvruO9zlcUJX+Y4+K4NmUqxoiIfLl+oelozom79SkURatWraJLly7UqFHD5veYTCYMBtv++datW3fPcyZNmnTX4zfX8V7nK4qSt6TJxLUPPyJ22TIwmXCoXBnPZ57O83KK5NpHReFJ4YMPPqBKlSq0aNGCkydPAnD27Fk6depEw4YNadmyJSdOnGD79u38+uuvDB8+nHr16nH27NnbngfQp08fBgwYQJMmTXj33Xfp06cPb7zxBo888ggVKlRg8+bNvPbaa1SvXp0+ffpY6xIUFERUVBQXLlygevXq9OvXj5o1a/Loo4+SkpJivfby5csBGDlyJDVq1KBOnToMGzbstnXMev7u3btp1qwZdevWpXHjxiQkJDzAn7SiFH8pBw5wolZtYpcswbFKFcqv+CVfAgKgjQ4pqv81bNhQ3uzYsWO35D1oe/bskbVq1ZJJSUkyLi5OVqxYUU6fPl22a9dOnjp1Skop5Y4dO2Tbtm2llFK+8sorctmyZdb33+28zp07S5PJZE0///zz0mKxyFWrVkk3Nzd56NAhaTabZYMGDeT+/fullFKWK1dORkZGyvPnz0u9Xm/N79Gjh1y8eHG2OkRFRckqVapIi8UipZQyNjb2tnW8kU5LS5Ply5eXu3btklJKGRcXJ41G4y0/k8Lw76IoRU16+DV5vHYdeaxqNXmsajUZOmSoTExIk9vWnpShZ2Pu+7rAHnmHz9Vi3XyUXx2c2s/0zrZu3cpTTz2Fs7MzAN26dSM1NZXt27fTo0cP63lpaWm3vDcxMfGu5/Xo0QO9Xm9Nd+3aFSEEtWvXxs/Pj9q1awNQs2ZNLly4QL169bJdv3z58ta8hg0bcuHChWzHPTw8cHR0pG/fvnTp0oUuXbrc9V5PnjxJQEAAjRo1AsDd3f2u5yuKcm/SaCRs7DjiVq+25lXZtZOTv4UiJ++kLHDBBKUrlMjzsotkUJBSrgHWBAcH9yvoutjKYrHg6enJgQMHcnWei4tLtrSDg7YiiE6ns76+kTaZTLe8P+s5er3e2nx0g8FgYNeuXWzcuJHly5czZ84cNm3adNc6K4qSdxK3/svlfpkfbaWmTSWtXmvCPjzIja9cF8q78MijFfKl/CLZp2CrOz0e5fa/e2nVqhWrVq0iJSWFhIQE1qxZg7OzM+XLl2fZsmXWuh08eBAANzc3azu8u7v7Hc97EBITE4mLi+OJJ57g008/vW0ds6patSphYWHs3r0bgISEhNsGI0VR7i715CmOV6tuDQhO9epRavd+Fm93If2LQwBE2gvsRjakxf8aYHDKn+/0RTIo3GueQkFr0KABzz//PHXr1uXxxx+3Nq388MMPLFiwgLp161KzZk1WZzwa9uzZk+nTp1O/fn3Onj17x/MehISEBLp06UKdOnVo0aIFM2bMuG0db7C3t+enn35i4MCB1K1bl44dO6qJaoqSA9JkImbRIs537w5AiZdeosrBAxx9ZAgJH+zmSewBMPasTP1JLfDzdM7X+ghbvvkWVsHBwfLm/RSOHz9O9erVC6hGyp2ofxdFuVXEJzOI/uorAOwrVqT0jE+4YvDB4avM+bxb6njQu3edPC1XCLFXShl8u2NFsk9BURSlKEu/dIlrU6aSmNFfV3rWLJIbNSdq2j4cyJyU5jkimN4lnB5o3VRQUBRFeUAs6emcCm6ETE9HODri1LAh/rNns3/5RZz+3kPJjBZ9+7fq4lumYEbyqT4FRVGUByB+wwZO1qmLTE8HoOIff/B1s7eJmn6UMucTKYmO2MfLEjilZYEFBCiiTwr3GpIqpVSLsBUiRbnfSlFyK+3sWS7364/x6lUwGAiYOJG0Zh05P+sAAzL2mwx11dP43cYE2hf8R3LB1yCPOTo6Eh0drZbPLiRkxn4Kjo6OBV0VRXmgzNevc+qRpta0z9ChuDzfm+iP9sHuA9wYQ+Q6PJhHvB9sv8HdFLugEBgYSGhoKJGRkQVdFSXDjZ3XFOVhII1GYhZ/T8S0ada8Cps2sntPAkEf7bPmXWkdQJPHC9+insUuKNjZ2akdvhRFKRDXV64ibNQoa9pnyBC2lmqF2+yzBGXknW9fipYdK1JYvyYVyaBQ3JbOVhSlaEs5cpQLzz5rTXsP+B/XOvbi4s+naXRea7WI9rKn5pBgAu30d7qMbRLC4fdR0HY0lKycu2vdRrGbvKYoivKgSJOJE7VqZ8tz+eMfdv9wmibxWi/yZXcDVV+rhZe/W+4Ks1jgxBr4+WUt3XMpVHvivi6lJq8piqLkscRt27jc93Vr2n/iRHY516XmlydpAiToQHYqR9NWZXNf2MX/YGGnzHStZ+47INyLCgqKoig5kLBxI6FvvW1Ne73yMqefeAW3pWepmTEb+Xold2r2rZP7EZCpcTDlpqAyLhr0+ffRrYKCoiiKDYxhYZxp286aLhkyEJ7uxYalJ2m9NHORSJ8xjQl0c7jdJWxnNsHvI2D315l5Iy6CkydSSvbu2UNw8G1bf3KtSM5oVhRFeVAs6elEzJyZLSBU2LqVL+NrkzbzCK2vGQl10SFfqErglJY45CYgmE1waBl81iAzIPT9EybEgZMnR48epX2H9jRq1IiNWzbm8s5uTz0pKIqi3EHqqVOc79bdmi41dQoHghrj9slxXsvIO1bGiY5vNsx9U9HBH2Hl/7TXJYLgucVQvSsIQWxsLOPHj+fzLz7HYragd9Gzdv9a2rdqn7syb6NIBgU1JFVRlPxkiooi4pMZxK1cCUCJF1/Eacgwwt7fScX/zljP8x/bhEBX+9wVlhQNX7eH2PNa2s4FBu4HnQ6z2cxXX33F6DGjiY2JBQFe7bwY894YhrYemrty70ANSVUURckgpeRKSAgJf/4FgGPNmpT+8ku++uEM3a9l7ijo1L823hU8c1dYXCh8WjMzXSIIXtsAbn4AbNmyhYEhAzl0UNt1zbmqM4+88Qi/vvkrLnYut7mg7dSQVEVRlHtI3r2bS/3/h8zYt7zMV/M561yBI58dp7tRaxpKaupPlW6VctdUZEyB75+Bi9sy8/qsg6DmAFy6dInhw4fz888/A2DnbUfllyrzx6Q/CHTL/3nQKigoivJQS967l8hPZ5Kc0ergN2YM6Y88TvSCY3hzDmfgWDln2r9aB72j3f0XJCVs/QQ2vZ+Z12MR1HwSgJSUFKZNm8aUqVNITUlF2Al8Ovsw+t3RhDwS8sAW+FRBQVGUh5IlPZ2T9eprM4UBp/r1Cfh0DgfWXKL0gmPogRM6CzUHNeRRP9fcFbbrK1g3LDPdcRI0CwEhkFKyfPlyhr4zlNDLoQC4N3an5f9asvKVldjpcxGI7oMKCoqiPHTifvuNq+9kfkgHzvmMyKMliJx5lNLAFp2Jml0q0aFZmdwVlJ4EH5bKTLcbBy2Ggk6bDXDo0CFCBoXwz+Z/AHAs48irY15l2qvTcLXPZSC6TyooKIry0DDHxXGqySPWtN/oUUQEtCTur1DsSQZgnHM6C8a2Q6fLRXONlLDqDTi4NDPvjf/ArwYA0dHRjBs3jnnz5mGxWNC76vF72o9PRnxCzxo977/cPKCCgqIoxZ40Gomc/Rkx335rzXNespron8NxPaY12VzwcaDugHosdMnlENNrx2D1m3B1v5buMAFaDAHAZDIxd+5cxo4bS9z1ONCBVwcvprw/hdebvF4oNgZTQUFRlGLt0v/+R9I/WwBwqFEd34mT+HdZDLV+voYTgr8xUuPlWrSo4Zu7gpJjYFqWvVw6TIRmA0GnLZW9adMm3g55m+NHjwPgUsOFgN4BbBq4iTLuuWymykOFKigIIVyAf4AJUsq1BV0fRVGKrtSTJznf/Ulr2n/8e5wr1ZT4Hy5SK2OFnzOPl+Gl1kG5K8higUM/wdohmXnvngdnLwDOnz/PsGHDWLFiBQB2Je0I6BXA7+N/p7p39dyVnQ/yde0jIcQ3QogIIcSRm/I7CSFOCiHOCCFGZjk0Avg5P+ukKErxJs1mro4cZQ0Irm3b4r95J38dL4Xf+suURMf2Cs6U+qA5bXITEKSEfYthUglYNUCbfDbgX22dImcvkpKSGDduHNWqV2PFihUIe4Hv077M3TCXi7MvFsqAAPn/pPAtMAf47kaGEEIPfA50BEKB3UKIX4HSwDFA7fCuKMp9Sd6zh7Bx75F+XlsyIuDT2fx5zI06Mw/QFB179WY6vNuU5zxyuYrpxe2w8PHMdL0XodtnoNMhpeTHH3/k7SFvE3MtBgCPRzzwf86fQyGHsNfnss8in+VrUJBSbhFCBN2U3Rg4I6U8ByCE+BHoDrgCLkANIEUIsU5KacnP+imKUjwk7djBpT6vAiAcHSk1YwbrIwNI/ieOphjZjpnaPavTvZ5/7gqKOg1zsqwOUbIqDNgKBi3I7N+/n75v9GX/Tq2T2bGcI6+OfZVpLxfcENOcKog+hdLA5SzpUKCJlPJtACFEHyDqTgFBCNEf6A9Qtmwe7GikKEqRJY1GLvTqTeqRzBZq99m/kLDuCi1I4CKSxaUNjHy7Re5G9pjSYMlzcG5zZt7Q4+CuzUGIjIxk9OjRfL3ga5Cgd9PT/PXmfDXmK6p4V7n/cgtAoepoBpBSfnuP4/OB+aAtiPcg6qQoSuGTtGMnl/r0saa9351F+iknLOuuAHDIx442bz5Cc6dczgjePBU2f5iZzrJOkdFo5PPPP+e98e+REJ8AevDu4M3qOatpXql57sotIAURFK4AWcdfBWbk2Uwtna0oD6+08+e59uFHJG3dCoBb1+fY59oBt1Pa8bUldHR4pQ5P+LvlrqDwIzA3ywd7o37Q+WNrcsOGDQwcNJBTJ7SCXWu58tKol5jdczYGXaH7vm2zgqj5bqCyEKI8WjDoCfTOyQWklGuANcHBwf3yoX6KohRC5sRETgU3ArR+A+//vc15pya4XUihZopkPUZSmgcQ0rVG7gqyWOC7bnBha2beO6esS1qfPXuWgYMHsn7tegDsfe15cuiTzA2ZSwmnErkruxDI16AghFgKtAFKCiFCgfFSygVCiLeBPwA98I2U8mgOr6ueFBTlISFNJiKmTydmkXUQI+4zfyb993BKk8JeTOyu4c6EF+pj0OdilL2U8O8M2DgpM+/JL6Ge9p01MTGRyZMn88mMTzAZTegcdfh09eHnaT/TKqjV/ZdbyKhNdhRFKbTSTp8mdMgQ0s+cxaFaNXxHj+fiNhMe17Q9D9Y39KRn12q45WZJa4DYizCrTma6xpPw7ELQ6bBYLPzwww+8M/wdIq9FAlCpQyXmfzqftrXa5q7cAlLsNtlRTwqKUrxJi4XoefOI+nIuMj0d/0nTOH/WB/sVCXgAf2Ok7oB69AvyymVBEjZPgX+mZOaNjbAOMd29ezcD3hrAvt37AHAq70THQR1ZGbISncjXub8FRj0pKIpSaEizmZhvFxExfToALi1akRj8Ou4XUq3nLK3hypAX6uauqQggIRxm1gFzmpZ+bQOUbQJAeHg4o0ePZuHChQAY3A349fBj16e7KOVW6k5XLDKK3ZOCoijFT+LWrVzu19+admozEJ1nbWtA+LOCM71ercdwO33uCoq/qj0d7F8M0gKVOkDPpWCwJz09ndmzZzNh4gSSEpMQeoH3o96s+GwFLSu2zF25RUSRDAqq+UhRio+EjRsJfetta9rxkS4cDniCBlL7eEpH4vXeI7zqnMvlIcxG2D47syO5Ugfo+L51j4N169YxePBgTp8+DYBbXTea9G/C8r7L8XDwyF3ZRYhqPlIUpUCknTtHxLTpJG7erGXYOePWeSYAyUh+II1ebzeiWqBn7gqyWGDjBNg2KzOv7VhoPRyAU6dOMWTIENatWweAvb89Ab0DWDZ8GY38G+Wu7EJKNR8pilJomKKjOf/kU5hiYxE6HS6tWhJZ9zX8Lhut5yT2rcHEyiVzX9ivIbBvUWa6WYi2P7IQxMfH8/777zNr1iyMRiM6Rx2+3X354cMfaFehXe7LLqKKZFBQzUeKUvRIKYn55hsipmuzgu0rVSRp6Ee4brhuDQgRFdyo368ugbndgSziOHzbBZKjtHSljtD7J9DpsVgsLPr2W0aOGknEtQgQ4NnSk0ffeJTvn/8eO10uh7cWcar5SFGUfBf/559cGRhiTetGvY/uhB9OGR8/kfaC8sMb4emWyyWtpYS5LeBaxgJ5ZZrAiyvAQVuhdMeOHYSEhLB7924AnCo50fbttizqv4iSTnnwZFJEqOYjRVEKhDEsjDNtM5ti9CV90DefjMPxzCeBtNdrUL+Sd+4LO7IClr+ame73N5RuAMDVq1cZOXIkixcvBsDgacD/OX8Wj1lMm7Jtcl92MaKCgqIoec547RqhA0NIPXQoI0dw5bUvqRaTObcg+tEy1G0XlPvC9nwDf4wFY5KWDmoJL60CvYG0tDQ+/fRTJk+eTFJSEsIg8O7kTb3n6/Frz19xMjjlvvxipkgGBdWnoCiFkzSbiVm8mIgpU615J4fNxHDGnnoZASHRXlB1YvPc9xtc3A6/9IP4UC3tXQn6/AZu/kgpWfPrrwwdOpSzZ88C4FbfDf9e/nzd62taBRaftYrymupTUBQl16SUJPz5J1dCBgGgL1ECj8lTSftT2wg+EckSZzMTxrRGl9uZyElRML1iZrp8a3j6K+sqpsePH2fw4MFs2LABAIdSDgT0DmDsy2PpW7tv7souJlSfgqIo+UJKSdK27Vx+/XVrnveAAWyt8QRN/gznxse/y+u1mFQpD9Yp2vEF/DE6M+9/WyFAW8ju+vXrTJw4kTlz5mAymdA56fB9yhfvdt7899J/uNnncn+Fh4QKCoqi3JfEbdu43DczGLi2bUtUyHi2Lz9Nk03hAKRXdKdCv7q5L+zCNvhnKpz/BwyO0P1zqP0sAGazmW+++YYxY8YQGRkJAkq0KYHf037M7jabjuU65r78h4gKCoqi5IgpOprTzVtY065t26IfNIId356j8aJT+ACXve0Jfr0udiUcc1fYvsXwa+YSGLQaDm3HQEZ/xL///ktISAj79+8HwLmKMwEvBNAkuAmLH19cbFcyzU9FMiiojmZFefCkycSl1/qSvGsXAE7161N67nxWLznBI99eojEGIrEQ8GINmtbyyV1hKddharnMdMX28MzX4Kw1QYWGhvLuu++ydOlSAAxe2hDTFp1b8F7T96juXT135T/EVEezoih3JaXk6vB3iV+71prn99441sdUonVomjUvopk/DbpVzn2BB5bAqjcy00NPgHsAACkpKXzyySd89NFHJCcno7fX49XJC5/OPjQPas7cjnPV04ENVEezoij3Je3cOc490dma9n9/EgnNO/H3d4dofV0LCMfre9G+R3UCdbn8MA7dA1+31147ecFTc6HKY4AWmFauXMk777zDhQsXAHAPdiewVyCvtXqNoQ2HInI7xFUBVFBQFOU2TLGxnG7azJq2K1OGUst/4YtFR+k+ez9NEOxzkLQZ0piOnrnsNzClw29DYP/3WrpsU3j5VzBoS2UfPnyYwYMHs2nTJgAcAh0IeCEA1+qu/P7M75R2LZ278pVsVFBQFCWbqK++IvKTGdZ0ue8Xc86vElM+2kF/HAHB9acr0K1xLj+MpYQTv8FPL2jpwEbQ60dw0dYgiomJYfz48XzxxRdYLBZc3F1w7+aOV1svJrWcxNOVn85d+cptqaCgKApSSiKmTOX6smVYkpMBcG3dmsC5XzJ2zn+8ceUw/XHEIiDg3WACS+RyeYjos1q/weWdWvqpeVC3J6ANMZ0/fz7jxo0jOjoanU6HVzsvfJ/2pWaZmnzc+mOCPIJyV75yRyooKMpDLm7NWq4OH25NO9WvT9lvF3L0ShJ/vb+VN5K1tnpzQx/K9aiWu8JC98LOuXB8jbY3csM+0Gkq2GlNUP/88w8hISEcylgzqWLDiognBY5lHHmm8jNMaDYhd+Ur93TPoCCEqAJ8CfhJKWsJIeoA3aSUk/O9dneukxqSqii5lPjvNq6++y7mmBgA3Lt0odTUKUQmGVnwxX4ev2akI4LLPg406l8Pg1sutsOMOg1zMga7CL02C7nTVCjbBICLFy8yfPhwli1bBoBvaV+8enihr6tHCMGUllPoXKHzna6u5KF7DkkVQvwDDAfmSSnrZ+QdkVLWegD1uys1JFVRcs6SksLJRo3BZAKgxMsv4fPWW+g9PFix7QIN1lzCgCAVSYk+NfGulstlrQ/+BCv7Z6aHnQFXbR5DcnIy06ZNY+rUqaSmpuLk5ET9nvW53vQ6Do4O9KzWk2HBw9Qw0zyW2yGpzlLKXTcN9zLlSc0URXlgzIlJxHyzgKgvvrTmVfpnM3Z+fhw4FUnJj7bSGABBpI8jdQc1RGfIxYdx5EnYMh0OLwO9PfRZB2W0PY+llCxbtoxhw4Zx+fJlAGp3qE3a42nEe8fTwKcBn7b99KHa+KawsCUoRAkhKgISQAjxLBCWr7VSFCXPSCmJ/upror/6CktCAgABH0zG85lnAPhx4X5anEy0nu/2VEUCm5S6/wLjr8JX7SHhqhYM6vSErrOs/QYHDx4kJCSELVu2AFCzdk2SuyQjq0o8DB68VOMlBtYfeP/lK7liS1B4C5gPVBNCXAHOAy/ma60URckTluRkTjZoCIDO3Z2SA9/Gu29fdI6OpEckc2jxEVpEapPQYhr7UufpqvdfWHoSrB0Kh37MzPvfFvDVlpyIiopi3LhxzJ8/H4vFgoeXB81fa87FGhdx0bkAsL3Xdgw6Nf6lIN3zpy+lPAd0EEK4ADopZUL+V0tRlNyQFguXXn2N5J3akE9DQACVNm1ECIE5Lo1Ln+9GF5mKP/ATafQb25JA1/vcHzniOPw5Hk7/kZnXZaY2skgITCYTX375Je+99x7Xr19Hp9fh3dEb3yd9ueRyiYoeFfm07adU8KiQ29tW8oAto4/MwHRglMzolRZC7JNSNsjvyimKknOxP/1M+Pjx1nSZr+bj2rIlADG7w0j+5Yx1n4OV1VwZ/HIL9Lr7WCJCSphVF65fzMx7ZgFU6wx22jyGjRs3MmjQII4ePQpAqfqlcH7WGcfSjlT3qs60VtPUnINCxpbntKNomydtEEI8L6WMAdQiI4pSiFiSk4mcNZuYRYuseTo3N6r8tx1hMJAemkDowiPYJ2ljRBaRxhNvNmRg2RL3V2BaAiztlRkQnlsMNbpZD58/f5533nmHlStXAmDnY0dAzwDcGrjh4+zD2qfW4mLncn9lK/nKlqBgklK+K4R4HtgqhHiZjE5nRVEKXurJk5zv/qQ1rXN3p/K/W9HZa/MK0sOTiJhzAHvgOhbmB9oz++2W91/gznmw/l3tdXBf6PyJdX+DxMREpkyZwscff0xaWhrOLs64dnKl5GMl0dnrmNV2Fu3Ktrv/spV8Z0tQEABSyp+EEEeBJUDZvK6IEKI6MAgoCWyUUn55j7coykNNSknUF18Q9dkcQJtv4DdqlHW10GuXrhP2wwlKxhkBmGqfztSxrZltf58duTHn4O+P4PDPWrrnUqj2hLUuS5YsYcSIEVy5cgWAtt3bcrXNVexK2PFslWcZ33T8na6sFCK2/HZY99uTUh4RQrQEuttycSHEN0AXICLrZDchRCdgFqAHvpZSTpFSHgcGCCF0wHdos6gVRblJ6okTRM2dR9LWrViSktC5u1P+l+XYlykDQFpkMpGf7AW0b1hmJMMc0lj6XnsM+hzOO0hP0uYZnFwPp37X8uo8D52mWDe82bdvHyEhIWzbtg2A+g3qU6tvLfa77McOOyY1m8RTlZ/Kk3tX8t8dg4IQop2UchNQTghR7qbDibd7z218C8xB+5C/cV098DnQEQgFdgshfpVSHhNCdAPeABbbfguKUvxJKTFdu8aZNm2z5TtUrkz5lSsQBu1P+cKJKAzfHgfAiOR6p7LUbxPEspwWGBcKawbBmb+y5/f5DYK0rTgjIiIYM2YMCxYsQEqJr68vA0YOYE/ZPexP1LbH3NZrG+727jm+X6Xg3O1JoTWwCeh6m2MSWHGvi0sptwghgm7KbgycyRjqihDiR7Qnj2NSyl+BX4UQv6E1UynKQ02aTER+NofoefOy5Vf65x/s/HyznGfh4OIjlDwZB8DfAfa8GNI45xvPxJyDbx6HxPDMvKCW8OSX4Kk9iRiNRubMmcPEiROJi4vDYDDQd0BfttTYwgrnFZAIM9vOpF2ZdmrjmyLojkFBSjk+4/+v5nGZpYHLWdKhQBMhRBvgacABWHenNwsh+gP9AcqWzfOuDUUpNJL37SP0jTcxx8VZ88ou/AaXpk2znZdwIY64uYe4sSBEYo+KvNTwPmYkH1gKqwZkpt/YDn41s53yxx9/MHjwYE6cOAFA07ZNKdW7FNvstqFHD8CW57dQwvE+RzUpBc6WeQqDgIVAAvAV0AAYKaXckJcVkVJuBjbbcN58IUQY0NXe3r5hXtZBUQoDKSUXnnue1MOHAfAZOpSS/fvdct7qP07T8O/Mb/R/l3ag55vBBOa030BKbRvMK1o/RNYO5BvOnDnD0KFDWbNmDQDeZbxxecaF+LrxJAhtPuvixxdTz7dezspWCh1bOppfk1LOEkI8BngDL6G1+d9vULgClMmSDszIs5mUcg2wJjg4+Na/FEUpwuI3bOBKyCBruvK/WzGUzL4onJSS+V/spvPlNGve7zXdeP2lejkv0JgC/0zVAoKDB7xzAuydrYcTEhKYPHkyn376KUajEQdnBzy7eOLd0RudnRZ8preezmPlHlNNRcWEzUNSgSeA76SUR0Xu/vV3A5WFEOXRgkFPoHcurqcoRV7qyVNEfvopiZs3W/OqHTls7UC+ISw6maOf7qGzKeNPsGUpAjtXzBwiaCuzSQsGW6Zp6bJNtU5kndYEZLFYWLx4MSNHjiQ8XHsa8Wvth+dTnth52vFclecY2Xgkdnq7+7hbpTCzJSjsFUJsAMoDo4QQboDFlosLIZYCbYCSQohQYLyUcoEQ4m3gD7Qhqd9IKY/mpNJqkx2luEg9fpzwCRNJOXjQmld+5Qocq1fPdp4pMZ2/Z+0hIMFINXSsN5h4dUxLDE45nHNgTIUdn8PGSZl5rd6FdmOsyV27dhESEsLOjHWTytQsg93TdjhXdOb12q8TUj9EPRUUY7ZssqMD6gHnpJTXhRDeQGkp5aEHUL+7UpvsKEVVyuHDhE+eTOrBzD+joGXLcKqdfe+qfZdisXxxCH8y+wlOlXWm3Zv30Z12aSf80hfiMsZ5NO6vzTfIeDoIDw9n1KhRfPvttwC4ebvh/pQ7ns08ETrB9FbT6VS+U87LVQqdXG2yI6W0APuypKOB6LyrXs6pJwWlqJIWCxd7v0DKgQPWvLLfLcKlcWNr2pyYTuzy06SeiEEbdKoFhKNeBpq/HUw75xw22UgJ/82BDWO1dMNXofMM0GnXTU9PZ9asWbz//vskJCQg9ALvx7zx6eqD3klP7ZK1+e7x79SS1g+Jez4pFGbqSUEpCqTFQsr+/VxfuZK41b+CUVt2osK6dThUKG89zxSXRvTCIxjDk2+5RsC4R9C73Ef7/a6vYN2wzHSfdRDU3Jr87bffGDJkCKdPnwbAra4b/r38cfB3oJpXNZZ0XoKdTvUbFDe53Y5TUZT7ZAwL40zb7AvAuTRrRpkFX2drlzfFpBI+bbc1/SNpzCGNoxMfw8XhPv5M48Pg6w4QH6qlg1rCy79anw5OnjzJkCFDWL9+PQD2/vYE9A7ArY4bS55YQm2f2jkvUykWbPpty1iawi/r+VLKS/lVKRvqo5qPlEIv5rvFXPvwQwAca9Sg9IxPsA8KuuW81FOxRH1zBIDlpPOlPp0NQ1sxzPs+lpaOvQAr+sPlnZl5w8+CizasNS4ujkmTJjF79mxMJhM6Jx2+T/ri3d6bXS/vwtnO+fbXVR4atkxeGwiMB66ROepIAnXysV53peYpKIWV8epVrk2dRsIfmbuQefboQcD7k245VxotRC89QeoxrYvuHZJwr1GSky/f9qn+HgWnwJ/vwa75WjqwEZRpAo9OBiGwWCwsXLiQ0aNHExERAQJKtCqB37N+jGk/hheqv6BGFCmAbU8Kg4CqGR3MiqLcwfXlywkbOy4zw2Cg8tYtGEpkLvkgpcSSkM7FPy9iv/uaNb87CbzdtTqvNi9PjkSegrktwJw5kY3Hp0OT/tbk9u3bCQkJYe9ebcaycyVnAl4MoEGDBnz/xPfY6+1zVqZSrNkSFC4Dcfc86wFSzUdKYSKl5Ez79piuhgFQ+rPZuLZsic7R0XqOJcVEwpZQEv7WhoPe+Bj+hjS+IY1do9vj6+5486XvzJQGc1tC1MnMvG5zoMFL1uTVq1cZMWIE33//PQAGTwP+z/lT97G6zGo3i4qeFe/vhpVizZZ5CguAqsBvgPXriJRyRv5W7d7U6COloMX9+itX3x1hTZdfvQrHqlWt6bSL8UR+efCW900mBZ9aJZnYuz66nOyPHH4E9nwDexZoaZ0dNOwDj30IBi3UxCbG0uatNhz56QiWNAvCICjZqSQlu5RkUfdFNPJvdF/3qhQfuR19dCnjP3syv+AoykPv8ptvkbhpEwAe3bsR8MEH1mUpjJHJXMvY6OaGKCyMJoU2bYKY26lazgpLT4IPb1r5tEJbeHmVNZlmSmPGohmMHzUeY6Q27NWtgRsBPQOY+exMHi//eM7KVB5KtkxemwgghHDNSNu6wY6iFDvJ+/aTuGkj0d8sBIs27qLiX39iHxhoPceSbMwWED71EfwSqbXAvt22EsMeq0qOnP0bFj+Zme48Axr1tSZjU2P5YPUHzH1/LklHkwDwLOfJsPeH8UjrR2gV2ErNNVBsZsvoo1poq6J6ZaSjgJdzul5RXlJ9CkpBuPRaX5K2b7emPZ56Cp/Bg62b3UizhbCPdmFJ1L6lxzbxpevOMxCpnb9leFvKeudgyKfZBH+N12YjAwT3hS6ZrbZ7r+3lpV9eImJVBNEbo8ECBhcD741/j1FDRmEwqGlISs7Z0qewHRgjpfw7I90G+FBK2Szfa3cPqk9BeRDSQ0M526GjNV1+xS/YBwWhc878gDeGJ3FtpnU1GL4glSWkA9C4vBc//y/7xjh3JSUs7Zm5JzJAyAHwyhyZtOPKDp4Z8wzXVlzDnGAGAf379+eDyR9Q8qalthXlZrntU3C5ERBA2wxHCHEfs2oUpeiJ+vJLImfNtqarHT6EsMtsirGkmoj48iCma9rSFMcw058k6/Fpz9ThuUZZtw+5h+iz8NOLEHFMSzd9Gzq+b52JPOTvIaz+czVhP4SReikVgJYtWzJ79mzq1at3n3epKJlsCQrnhBDj0JqQAF4EzuVflRSlcLjQs5d14Tr/CeMp0bNntuOpZ68T9dVha3oMyRx11TGsWRV6NymHl0sOxmVYzLC0F5zOmPRWroXWiZyxX4HZYqb1vNYcWXSEuJ1a/4R3gDeff/o5zz33nJp4puQZm3ZeAyYCKzLSWzPyFKXYuvbRFGtAuHn3M3OSkbD3d2Q7vwPxpAJ7B3fA29XBtkIsFjj8M6z8X/b8F5ZD5czmqpSUFKq/Up1Lqy8h0yWOjo68++67jBgxAmdntSyFkrdsGX0UC4Q8gLrYTHU0K/kl9dQpznfrbk2XX73aGhCklIRtvozlj4vW46+TyAkszHupIY/V9Le9oONrtGairFx8YNBBsHexlrfkpyW8HvI6qZFaU1GPHj2YPn065cqVu887VJS7u2NQEELMlFIOFkKsQVvrKBspZbd8rdldqLWPlLx2czAAqHroIDp7rQno0LbLeK25YD2WguQxEhjQpiJrO1bBoNdhE2MKfJAleOgdoO8GCKgLWZqADh06xIv/e5HDO7TmKacyTqxauIpH2z96fzeoKDa625PCjT6Ejx9ERRSlIMQsWULCHxtI3pm5qmjpmTNx7/QYAGaLZPKkzfRNzfxTWRHsyVtP1+KcrTORY87B7Pq35g85Ch6B2bKio6MZO24sc+fOBQl6Fz09QnqweNJiNcRUeSDu+Fsmpbwx+6aelHJW1mNCiEHAP/lZMUXJb6dbt8F0LXNRulIff4xHl84ApIcmEDHnAAB9M/5MYp6qQJ0mpXPWlnp5NyzokD2vbi94am62LJPJxKhpo5g3fR4J1xNAgFd7L76Y9gXPN3g+p7emKPfNlq8erwCzbsrrc5s8RSkS0i9e5OxjmXsNB/24FIeqVdE5OSGlJPlKArFzsq9X5PlqTQKretlWwM0jiW4YdAhcfcHOKVv2hr820O3VbqSFakuLuVR3IaB3AKfHns75zSlKLt2tT6EX0BsoL4T4NcshNyAmvyumKHnNeO0aSVu3Wpe3titblgpr11j7DeI2XiLhz8xO5CQkF7sF8WizsrYXcn4LLOqaPa/Ht1DzqVtOPXjyIE/3e5pzW7UR3nbedtTrW48hfYbQq3qvnN2couSRuz0pbAfCgJLAJ1nyE4BD+VkpRclr0d9+S8SUqda0a4f2lJmjLR9hjEgmduUZ0s9r4///xcg+zIx9rzVVnW2YayAlhB2A+W2y578Xa510ltWRq0cYOXEk6xauQxolwl4Q0DWA09+dVkNMlQJ3tz6Fi8BFIAfz8x8MNSRVsVV66BXOdshs03d+5BFKPP8czs3akfDPZeLWX8h2/jMksHnSY/S019/74qZ0SIuHlQPgzJ+Z+Y9+AM3evuX0KwlXaDu6LeeXnMcUYwKgYpuKrFuwjioVqtzX/SlKXrNlQbxHgM+A6mhLZ+uBJCmlez7X7Y7UkFTlXiypqYSNHUf82rXWvPKrV+NYtQrGqBTCP9qV7fx5pLKUdA5OfAwnWwLC7UYUNR8MbcdY9zW4QUrJrDWzGDN8DMmntOUwAqsGMnP2TJ559Jn7uj9FyS+2dDTPAXoCy4Bg4GVAfa1RCi1TTAynmzW3pktNnYJHd20OQtZ9DhaTxkaM+Pi7UrVCKba3rYiLg43DPrMGhFbvQqthYLh1JvOuM7vo1K8Tsf/EggQXTxc+nfYpr732Gnq9DcFHUR4wm/4CpJRnhBB6KaUZWCiE2A+Myt+qKUrOpR47xvmnM799V966BYOPD+mXE4j4/IA1PxoL80jL+XLWoC1pfcOE2+9UazQamT5rOuPGj8OSbAEdtO7ZmhVzVuDlZeMoJkUpALYEhWQhhD1wQAgxDa3z2cbpm4ry4Fz76CNiFn0HgNcrr+A3aiQACVtCiVt33nreScz0JYnDEx7FzfE+Np/5XbsuT9x+Xueff/7J4MGDOXZMW+m0bHBZ1i9aT40aNXJelqI8YLYEhZfQ+hHeBoYAZQDVEKoUGukXLnC2U/atJn2GDif5YAQxy06DSdshbSNGxpNCGS8nzg57An1O9ka+YcNY2P2V9rpW9j+Dc+fOMXToUFavXg2AnY8dAb0COP/peXS3GYWkKIWRLQvi3Ri4nYK2WqqiFBqW5ORsAcG+YkUqrFnDldH/ZjvvB9L4kjRm96pPt7qlbr6MbY79Cts/014//RU4a81AiYmJfPjhh3zyySekp6ejc9Dh09UH78e82f3KbhUQlCLlbpPXDnObhfBukFLWyZcaKYqNrq9cRdiozK6tsgu/waVpU8I/zr4b33skswnT/fUf3PDL63B4mfbavTTUeQ4pJT/88AMjRozg6tWrAHg288Svhx92Jew48NIB9DrVmawULXd7UujywGqRQQjxJNAZcAcWSCk3POg6KEVD5OzPiPriC2u6yoED6B0duL7uHKaoFAAMfWvyyIL/AFg7sMX9B4T4sMyA0LAPdJ3Fnj17CAkJ4b//tOs7lXci4IUAnCs5s+7pdZRxy8Fua4pSiNxr8lquCSG+QQswEVLKWlnyO6Gtn6QHvpZSTpFSrgJWCSFKoK3OqoKCYpV2IY74vy9y/YdP0LmVwu3J+dZjYROyzzuY5WRiWUZA6FwngFqlPe6v0IRwmFFNe924P9caDmN0374sXLgQKSV6dz3+PfzxbO7J2w3eZkDdAfdXjqIUErZMXksgsxnJHrAjZ5PXvkWb6/Bdlmvqgc+BjkAosFsI8auUMmNjWsZmHFcUUg4eJGLuYfQe2sYyTsGv3/X8M5hZlpJsTX/eu8H9FRxzHmbXAyDdLPnsmB+TelUhPj4e9ODd0Rvfbr7onfU8X/V5FRCUYsGWjma3G6+FthFsd+ARWwuQUm4RQgTdlN0YOCOlPJdx3R+B7kKI48AUYL2Ucp+tZSjFU+rx40TNm4/UPWYNCNJiRmS00ztU9sSzcwWkRbIrMYU/vzlADfRMR9ulbNeY9vi6OdpeoMWs/f/qAfi6nTV7/WkjQ/aW5eTJ4QC41nEloFcADgEOLHh0AY0DGuf+ZhWlkMjRrh1SSonWvDMeGJmLcksDl7OkQ4EmwECgA+AhhKgkpZx78xuFEP2B/gBly+Zg9UqlSEg9cx3TtSQi5y3HrlQDcHyaGwNHUw9+h/crHSjRM/v+AkN/PsCKfVes6enP1qFHcA7a9C0WOL8ZFmdfyfR0tJmQDen8fsoInMTez56A3gG41XXjrXpvqScDpViypfno6SxJHdpSF6n5URkp5Wxg9j3OmQ/MBwgODr7j6Cil8JJGM5Z0C3qX7BPHwj7ahTlO21PArlT2Jh/nRn4ETvnKmrZYJJdikmnz8eZs520b2Y7Sntn3K7ij/T/A6jdvyY5Pk0zeksaMXSbMJgs6Rx2+3X3x6uiFzqBjVfdVVPSsaFsZilLE2PKkkHVxeBNwAa0JKTeuoE2CuyEwI88mapXUointUjzXV57BGJYEgN7bEb+36pFyPIbYFafBnD3G6xzj8HmzHXa+maOG4pKN1J106/iDZQOa0ijIxuUjpIQfnoUzf2XLttR/he9iGzBq1CjCwxMA8Gzpid8zfrza9FVGNB6Rk9tVlCJJaC1C+VyI1qew9sboIyGEATgFtEcLBruB3lLKozm5bnBwsNyzZ8+9T1QKVPrVRCJm77fp3JRd8zBd1Rasq3b8GEIIjlyJY9iyg5zI+KC+2aTuNXm5aZDtFdr/Pax+S3vtVAJeWcvOM1GEvDOCXbu0UUxOFbUhpq8/8TqjmqhlvpTiRQixV0oZfLtjtjQfVUAbOvoI2iik/4AhNzqJbXj/UqANUFIIEQqMl1IuEEK8DfyBNiT1m5wEBPWkUPhJk4X4jZdI+PvyLcd8Q+qDWWZboA4g/dzfmK7upfyvq3GsUgWzRaIX0OWzf2+5xtCOVejdpCwlXW9dmfSuTv+VGRD6byZM+DNyxEi++04bHGfwNODXww/Ppp74ufqpgKA8dO75pCCE2IE2PHRpRlZPYKCUskk+1+2e1JNC4RU6custeQFjm6Bz1BH5+efEfL0AiT0l3x5F1KdjredUPbAfi509lcasB6CEsx2xyUYA3mhTkSEdqmBvuM9lI9IS4KNA7aVJMtNtDJMnTyYxMRFhEHg/6o1PVx/0TtropsOvHL6/chSlkMvVkwLgLKVcnCX9vRBieN5UTSmOLBkf4gDSmEjyttlU3rgc4WzHiepZVwo1WgOCz+BBlBygjeYZuzLzw/hGQPjyhQY8Xjvg/iokJUz0zHgpWXvKxNA9pTlzRhtA51bfDf+e/jj4OfBi9RdV34HyULMlKKwXQowEfkRrPnoeWCeE8AKQUsbkY/1uSzUfFU7SLEnaHc71VWcAcKxqIHLqUABONrztlxIrxxdf4dtt55mw5pg1L6RdJeJTTQzpWAUPp/tY4vqG/dp3mhNRZkJ+T+XPs2bgDA4BDvj39setthtv1H2DN+vdOhJJUR42tjQfnb/LYSmlrJC3VbKdaj4qXK5O3oElMfMpIWFV/9ueV2bB17g2b47xWgRCJ4i0d6XpR5uyndOppj9zX2qY+0pdO0rcp02Z+E8as3ebMJst6Jx0+D7pi3d7b4RB0LVCVz5s+WHuy1KUIiJXzUdSyvJ5X6XcUU8KhYcl1YQ5Pp1rM/Za85wb+JL0d+acguonjt/2vXZ+vgC8Nit7/8OxSY/hbJ+jeZW3ZY4NZWG/YEZvTCMyWYKAEq1L4PeMHwZ3A9W8qrGs67Jcl6MoxYkto4/sgDeAVhlZm4F5UkrjHd+Uz6SUa4A1wcHB/QqqDsWJKS4NmWrCzs/lrueZk4yY49NJ2n6V9PAkdPY60s5m347Sp38dHCp4EP2ltkpJ1b33fpI7HhYPwIUpne/zDm61bfFHhIwYy74wbYMd58rOBLwQgFOQEwdfPkhEcgT+Lv55Vp6iFBe2fB37Em0RvBvrFL+UkXf3VcmUIiH6xxOkHIjMlufSyB/P7hWJ/+sSwiCI/+vSPa9jH+SOT/86CJ0g/dIlTGFhCAcHdC53DzS7L2hdUjUCbF1f8e5CQ0MZMWIES5YsAaC0pz09R/fnd5+/EUIwsvFIdEKnAoKi3IEtQaGRlLJulvQmIcTB/KqQ8mBIo4Xwj3djjku/5VjS7nCSdoff9f0lnqmMsNcjHPUYPB1IObCNU02b4jt4EFFfaktW+QwefNdrHAq9To+52vLW47rkbv/i1NRUPvnkEz788EOSk5Nx0MPwZva8uy6KZr80QyBY99Q6yrirfQ4U5W5sCQpmIURFKeVZsE5mM+dvte5O9SnkXti03VgStIDg3acmjlVKAJC0M4zrq8/ecr7nU5VwaewPEsRNexsnbd/OlUGDAAifOMma7/HknVdD+fd0FC8u2AlAj4aBNK3ofV/3IaVk1apVvPPOO5w/r42JeLq6gY87OlJ+wFJarc1sklIBQVHuzZagMBz4WwhxDhBAOeDVfK3VPag+hdyRUloDQqkJTdE5Zv4auDYthXNdH6QEvYsdlnQzwk6Htmo6cJu97i+91tf6Wjg4UPLttyjRqxcmRycMgNFswU6fOeEsKc1kDQj9W1Vg9BPV7+s+jh49yqBBg9i4cSMANatVZlajK7SvoN3PkRIBxKbFArDrhV13vI6iKJlsGX20UQhRGaiakXVSSpmWv9VS8oM0mjHFpHLtU60T2L6cA0JnueU8nXPmnACd/d33GI5bs8b6+sYoo5ikdCq+/2e2856qX5ppz9bBTq9j0I8HAGhb1ee+AkJsbCzjx4/niy++wGw2U8JJz6Q2dgwIDseg036lYwdsYfg2bXLa9FbTcTLYuHKqojzkbBl95Ai8CbRAm7y2VQgxV0qZL8tnK3lLSknqsWgcKpXg6vjt2Y5Ff/4/omelEvjlF7i1bZvja6edPs3V4e8C4PXKy4C2immDmwICwMr9V3jxkXJEJ6bx1/FrAHz5Ys7mIZjNZr766ivGjh1LdHQ0Op2ON1qXYlLjeOI87PnNwYG2yclsfnwSY/540fq+TuU75fjeFOVhZUvz0XdAAvBZRro3sBjokV+VuhfVp3BvUkoi5x4i/WL8LcfMsRdI/idzslboG9pM3vKrVuJYrZpN17/U93WStm2zpv1GaQvHbTp5zZo3sVtNjGYLvu6OhCzdzzNfZgalZxsG4mh396eQrLZs2UJISAgHD2pjHFq3aMasRyKo6xLBJYM9TwaWyjjTG45lzpEIdA20uQxFUWwLCrWklFmHhvwthDh2x7MfANWncHvm+HRMMSk4BHmQsPHSHQLCeZL//QTPXj3xHzeO082aY75+HYDzTz5F+RW/4Fjj7iOB0kOvZAsIVQ9oy2LHJKUz/feTAPw5pBWV/bSdXOOSb53SMqFbTZvu6dKlSwwfPpyff/4Z0Hbb+3joizwb85m1n+On9kPg1E/Z3tevdj+C/YOp71vfpnIURdHYEhT2CSEekVLuABBCNAHU2hKFTHpYEhGztL4Cv3caZptbYE4IRyZeI+3EaizxV3Hv0pmA8eMBqLx9G/Fr13L1XW0RuPNPP4Pv8OF4930t2/WT9+8n4uNPcKxaldiMOQCePZ7Ff9Ik1h0OZ8afJzkbmWQ9v6KPq/W1h7Md/41qR0KqiTIlnIlPNeLqcPdfvZSUFKZNm8bUqVNJSUnB0dGRkSNHMnz4cJynBRBhMGCQktg3tvLdul4A7Oy9k0sJlyjvUR4HfQ6X1FYUBbBt7aPjaJ3MNz5lygIn0XZhk1LKOvlaw7tQax9p7rSJjX2QO7EL/4clLnPWcfnVq3CoUiVzNFGG1JOnON89cwip/4QJ1r2QzfHxXHzpZdJOnsz2nio7/uMaDjSbkn3dor1jO+Cd030OMkgpWb58OcOGDePSJe1X7rnnnmP6lI8oG1gKzGkwpSz1g8pgynIPtbxrsbTL0jtdVlGULHK7dLbqpSukpJQkbA4l/o8LADhWLUHqyVjr8aQts7DExeFYqxalP/kY4eCAnf/tZ/I6Vq1C9RPHiV+/nitDhhI+YQLOTRqT9N9/XJv0vvW8Ml9/TdL1OP5McOCvvZEkpJmsx5b2e4Rq/m6UcLG/r/s5dOgQgwYNYvPmzQDUrVuXWbNm0bqyB8zPaAZqNZxEIbIFBIDvHv/uvspUFCU7W4akXnwQFckJ1dGsSdp+1RoQ9J4OlHy1FumXE0g5EUX0nPdIO6k9RfmPH499uXI2XdP98ccxXr1KxPSPOff4E9mO2ZUtywbncgxdexBIRdtRVfPLG01pWM7GPZJvEh0dzbhx45g3bx4WiwVvb28mT55Mv3790Ov1MMEj8+Qt0/nNLbNpqk1gGz5s+SF2+lwsra0oilXul6IsAKqjWZOw5Yr1tf+IRgDYBbpytmNja77vyBE41a6Vo+t6vfYa5vgEjGFXAXCsXgNTZCSHWnVj6M/a6J+2VX3o16oCvb/SJqFVyehUzgmTycTcuXN57733iI2NRa/XM3DgQCZMmICXV0aAOfUHANKjDK85peEgJductTkHO3vvxNnOOcflKopyZ0UyKCja7mb2Zd1IOZyG//BgwkaOIm71agK//MJ6TtlvFuDSrFmOry2EwHfI4Fvy52UMKQ3wcGThq1rg2fROa1KNFtwcc/ZNfdOmTQwaNIgjR44A0L59e2bOnEmtWlkCWHoyLHmOAw72vOQlAEcA7HR21POtpwKCouQDFRSKmPTLCaSevU787xeseQZvJ+JWrwYy5xwEfvH5fQWEm20+GYGrg4GoxHT2XNT6K7aPbGc9XiHLKCNbnD9/nmHDhrFixQoAgoKCmDFjBk8++WT2zu/4MJhRjQsGAy+VyuwH+aT1Jzwa9Ggu7khRlLtRQaEISDkRQ+qxaBCQtDP76qV6Tz1xa3+75T2urVvnutzYpHT6LNydLa9Hw8BbRi7ZIikpiSlTpjB9+nTS0tJwdnZm9OjRvPPOOzg6OmY/2ZgCM7RJdC+ULmXNVqucKkr+U0GhkJMmC9HfHr0lP/XA9xhKmklYtZ3r3DqsWOhtny18J28v3XdL3vtP5qx/QkrJjz/+yLvvvktoaCgAvVtVZur3GwksUwZTYgTjvmvPeUsKTs6+jH/8awLP/8cRe3vmeXsTn7GO3mftPlMBQVEeABUUCjFpkYRn2ebSqaY3EVN7g0lbj9B4Ifv5JV54AVNEBAEfTM512WkmM9vORAMwpEMVwuNT+fCpWjl6Sti/fz8hISH8+++/ADQoX4LZbdJoXvYarOjCclMkE0tmLJktgJQrrF3Rm8rXTjG4dGaTURP/JrQp0ybX96Qoyr0VyaDwsAxJtSQZMcekIuz1+A8PRudiIOKDzAVqK2/fhs7ZmZP1tDH8PoMHoXfL+Sig24mIzyxnUIfKOXpvZGQkY8aM4euvv0ZKiY+z4MP2jrxaz4Q+YxXTjWnXmOjnY31PkzQzOx30fE4sZMkH6Fqxay7uRFGUnCiSQaG4D0mVZkna2etEfaONzPHsVhG9mz2W5GTrORX/3IAhY9hmma++Qqan5UlAsFgkey/FWndEm9Wzns3vNRqNfP7550yYMIG4uDgMOhjYxJ73Wjvg+ai2jAZ+NTFHn2bCmYUAfFzheR7Te0CDl2m2pBkJWfZd+LDFh6SYUuhWsVuu70tRFNsUyaBQnBnDk7g2M7Mt3611IE61tCaWxIyZvn6jRmJfJrN93bVlizwpe93hMIYtO0hyurax3stNy9Gumq9N792wYQODBw/m+HFtT4XHWjbi09rHqF65AnSdCRW1EUspphQ6/fIp1/V6PPROPNZyrPUa33f9kb8Pf4djQD2ervy02gNBUQqACgqFjCk2yzYVOvB4vLw1Gf6+1lfg2qZNnpXXbc6/tKxckmYVS/LmD5nByNfNgYndat6zD+Hs2bN07tOZk/9q6yKVr1CeWTNn0eVgf4RJD11mMOD8z+zfMYp+dfrx75V/iUmNoV2ZdrxZ781s16rgU5sK7abn2b0pipJzKigUMsZwrYnIoZInrk29ufhKH8wJ8VgSEjHHxiLs7GxesuJuElKN9P9uL4dC4zgUGscvezNnR/dvVYHnG5W5a0BITEzkgw8+YMaMGaSnp6Nz1OHT1YfqPaoTVvoiYm8yfzs7EbJtiPU9s/bNws1Oa+Ia13QcJZ1K5vo+FEXJWyooFDKpJ2IA8H6hOmlnjpO8cyeOtWtjvHwZgBIvvZQn5ew8F8N/57TRRY52OmqVdqdjDT8mdb/704HFYuGHH35gxIgRhIWFAeDZ3JPW/VpTKagSu8N38/2JJQS4OvObi8st708wJtCzak8VEBSlkFJBoRCRFmndGEc46Ek7dRqAkm+9yfVly0k9cgTPHs/mSVkXY7QnkhVvNqNB2RI2vWf37t2EhISwY8cOABo1akSXoV1YlraMOpXqMKnxaOZvm8i8i+sY7VOSDiVqQewRhgcPp45PHd755x1MFhN1fApstXVFUe5BBYVCJPW49s3dqa4PMi2VsDFjANC7uVHm8zl5WtaNfTQq+d57mYrw8HBGjx7NwoXaiCE/Pz+mTJnCyy+/TKMftIX4+tTqA7u/pv/mLzjj4816VxcOJF6iconKvFxT2795Y4+NeXoPiqLkvUITFIQQFYAxgIeUMm++DhcBllQTUQuOYEkxYYpKAcCtTRkS/vwTAI8nn8SpQYM8L3fr6SgA7HS6O56Tnp7O7NmzmTRpEgkJCdjZ2TF48GDGjh2Lu7s7n+3/jHRLOsF+wVTwqADn/gbgyYQk1ru6EGWMx0eUzvO6K4qSf/I1KAghvgG6ABFSylpZ8jsBswA98LWUcoqU8hzQVwixPD/rVFhY0s2knb2OKSaV9MsJ2Y7Z+TkTe0IbzVNywP/ua62huzGZLRwKvQ6Ak/3tl8NYt24dQ4YM4dQpbc+Ezp07M2PGDKpUqQJAkjGJ+YfmAzCodn84uR7Ste0467UeR4/QDST6VqNlYMs8rbuiKPkrv58UvgXmANZtsYQQeuBzoCMQCuwWQvwqpTyWz3UpVJL3XOP6r2etaZ2zAUuytouZ0AmQEuHkhH1QUJ6Xvf1sNLHJRiqUvLUj+NSpUwwZMoR169YBUKVKFWbOnMnjjz+e7bw2P7UBoIFvA+qd3wF/TdAOlG6Ic9O3eY+387zeiqLkv3wNClLKLUKIoJuyGwNnMp4MEEL8CHQHinVQSNoVjs7VDodKnkR+dRhjxtOBzxt10TkbwCy1SWt6QdTceSRu3QJ5+IQQEZ/KzI2nSTdZuBKrNVN9+HRt6/H4+Hjef/99Zs2ahdFoxM3NjfHjxzNw4EDs7e05HXuaxccWU8atDH1q9SHVrM2neDnFAnunaBf531bwLJtndVYU5cEriD6F0sDlLOlQoIkQwhv4AKgvhBglpfzodm8WQvQH+gOULVt0PoBiV2gjiXwG1LEGBIcKHtiXcUPoBJY0M3alXXF9xI9LL/VD5+aGa4u8makMWh/Ckp2X8HVzwE6vo3qAO9X93bFYLCxatIhRo0Zx7do1hBC89tprfPjhh/j5+Vnfv/78elaeWQmQrUko+MgasGSs0hqgRhUpSlFXaDqapZTRwAAbzpsPzAcIDg6+dc3oQsgYmblm0fVVmU1GPv3rcD05nSE/HSApzUyP42up9bm2+Yz3a69S8o03clTO5ZhkRq88TJc6ATzfSAuYy/ZcZs2hMEJjtTosH9CMst7ajmU7duwgJCSE3bu1PROaNm3KY4MfI9EvEW8fbxYdXcSmS5s4FXuKRGOitZz3tr0HwMd2QXhYLmmZ9V7MUV0VRSmcCiIoXAGyLowfmJFns6K2Smr6xcyOZJ2z9iN3eSQAgFPXEvn7ZCSuDgZqbVphPc+udM5H7RwKjWPr6SjORiRag8LHG05yLcuKp/4ejly9epWRI0eyePFiAEqVKsXUqVN54YUXqPNdHbgAwxsN57dzv3E85rj1vXV96uJkcMIiLbQs3ZLaCabMwt0yl7pWFKXoKoigsBuoLIQojxYMegK9c3KBorRKqik6hdiVWtOR/7uNMHhpu4z9vOcym3/YS3RiOgB9Qrdle59jjRo2Xf9afCpT1p8gzWTm6nWtnf9qXCp9Fu7C2V6fLSC0rujJjI+nMXnyZJKSkrC3t+edd94huHcwF9MuZrvuxP8mEpoQmi3v+ye+z0yc/hN+eyczXaGNTfVVFKVwy+8hqUuBNkBJIUQoMF5KuUAI8TbwB9qQ1G+klLduLXb36xaZJ4XUM9fBLDH4OqN3s7fmL9p+gQtRSZTydKJuGU8e+/P3bO+zCwiw6fq7zsewcv8Vynk7Y59l2enNJyOp7OuKr5sD1+JTKRF5iN+/n893V7QP/+7du/PJJ59QsWJFai/SOpwH1M1svbuScAVfZ18e9X2Uo9FHqeldM3vB+7+HBG2ZC0o3hJI523NBUZTCKb9HH/W6Q/46YF0urltonxRSjWb+OxeNySyp6eqAU5g2dt/n9VoIO+1DOy7ZyNGr8XSo7svXr2gzgo9/loBnr57Mq/8sP++5zMJrqTQu73zXOQpmi2TbGW0S2oJXgqnk68b7a4+x4N/zjHmiOv1aVeD48eMMHjyYDRs2AFC9enVmvT+CjgEJpKcdZeuxPdbrbTm6BEedHc+W78KIFpOyF5aeDKf+AIu2rDYJYVCiPLy9K09+boqiFA6FpqM5Jwrzk8LqA1cY8cthAL5z8qBCigSdQGSZJDZ7k9ac5O5oB4Dx6lUALHFxuDsaSE438/z8Hfw1tBWVfO+8cc7Oc9H8uFsbyOXqoF2rUVAJFvx7noqegiFDhjBnzhxMJhMeHh5MnDiRN998E7sPSsJh+M3Vhfd8vK3XG3ZgplavC9vg5oFPe7+FP0ZlzwtsnKOfjaIohV+RDAqF+UkhJWODmiBvZwwJEvty7nj1qorOMfNHnWLUzvngKa3ZxpKq9QW4tm7NwPaV8XS2Z9LaY6SkW+5eVsZ1Pu/dAH8Pra+iY3VfBpW+SM8OrxIZGYkQgv79+zN58mR8fLJvc5nqWRZIYl7HeTgbnLH/phM6i4mKpSvcWpgxYwTV6xshY0tNSuR+CW9FUQqXIhkUCqM/joaz/UwUx8O0kUYtTXo8jJIryen4emof2ElpJub8fYY9F2Io6epgXWLCkpgx3NNgwE6vo6yXNmQ0LC6Fq3EpbD8TRcvKPnSokTlvYNOJa/ywUxsOWsZL26Hs33//JSQkhP379wPQokULZs+eTf362h7OHPsVLmxljYsz80p4cNFOa9qqtnYkXoFNIDUVkHB+M4QdhIC62vsO/QynMvo8/OuAIbNvRFGU4qVIBoXC2Hw066/TnLqWgKujgeoB7vSO1eEgzSyPvE51oxlHOz0HLl/ny81ncXUw0KJS5n4CpihtddQb/QfertqH7uXYFJbvDeV4WDz7L1/PFhQ+23SGw6FxWkBIiqZ37zdZunQpAIGBgUyfPp3nn38+e5/E5o8g6jSjy2XvxHa7ehAiToGzFxgcIf4K7FmobaMJ8NdESIqA0sGZTwmKohRLRfIvvDA2H0mgbTVfvno5GICwabu54CVYeDWBETfOyZhqt/DVRjQK8rK+V2SMGrILDAQgyFtbk0gvMpe4ljdN05MSGpd1pUbUFpo26EpycjKOjo4MHz6cESNG4HKbDW6QEqo+DqkHrFkfX4vEDqDdGGg2MCOzSsYdZbm7Os9B989z9kNRFKXIKZJBoTAxmi1sPhlJfIoRSmjNOOlXEpFpZsjoR9h7MZbopHROhsff9hppZ87eNv/wlXjtusD1lHR+Pah1SEspOb3zLy6un8fSKG1Y6DPPPMPHH39M0J0W0Is8BZHHkd42Pl3FnINLO8DBTXtyUBTloVAkg0Jhaj7adiaKft9pwzqbV9JG8kQtPIIlyUhSCe3H+8LXO7O9x9PJLls65ttvAdB7eADgaKfH0U7HL/syJ49djkkhZOl+0iMvELtxPqkXDwFQq1YtZs2aRbt27e5e0RWvA3DGwQEy57PhacnozHbKfHLByQvOb4GL28Gt1K3HFUUptopkUChMzUdpJu1Dde6LDWlXzRcAabLg3MAXXS0P+O4aAH1blKdX47K4OhisI4VuEI6OODdujH05bTSPk72ebSPaEZusPSWU8nTk2PmrzJ72IUsWfY3ZbMbDswQTJ07krTffwGCw4Z/RmAq+NUlvHgLrX2R48HDalGlDWbMEsxG8K2ae23cDbJkG2z+DtDht6Gn78bn9USmKUgQUyaBQGJXxcsLekDmjWOdkQGTZ1czb1f7OW1/qBAZf32xZ3q4OeLs6YDabmT9/PuPGjSM6OhqdTsebb77JpEmT8Pb2vv31bkcIKFnJ2lFczr0cZd3vsMqsozu4ZAxfFXpttrJe/aooysNA/aVnseNcNGW8nCntqfUNmOPTSA9LQudkwBSdinMtb/ZeieNMROaKoUevxuOJwO5oDImXkpFpJqRJEmn4DUvUFsY1iSLV7EBpsz1r99fAw6cLLQNbYoqJIfHvv5EWC5aExFvqciXxCgtWLeCbyd9w8aS2NEXl0j689GJzarR3ZPO+D+94HzohaOtWCU+Ddh/EhULkCfCpysHIgzn7oZjS7n2OoijFRpEMCvnVp9Bz/g6q+bvx++BWAMT8eJK0c3GZ5eqrMeDXA0RlLGJ3wwAccN4YyvUsebEOmyH9IEEekGx0wg4zF8NOMOrQOrb32k7K4sVEfznXer6hZOYQ1YsXL9Lp1U6c+PsEAI4l7elYvwqtPMtyLS6BCVc33PNeQmJ+ol/cTR3brv58tv8zALwc79FH4Jox/NWYBK6+dz9XUZRio0gGhfzsUzgRnrnMddaAACDNknSThR4NAxn6aBVrvu7vUEz7IpAZM5C9X6zOtXQHbkQJH69aRKaEQ5o2isdkMYHRiLC3p+KGPwAw+PqSnJzMtGnTmDp1Kqmpqejsdfh09qHcE35U3euDiBHU8arGO499fNd76PjHSxibvgXVX9KeEL5/WjvQaQqGn1vTsVxHavvUvus1qNsTKrQFiwncS93rx6YoSjFRJINCQXNxMBDg4WRNXzfoMWWZJKZzsYP0zLRdltVLsxECO39/pJQsW7aMYcOGcfmytpZR1bZVCewZSLRjNFnfbW/vir9/vXtX0tEdPEpDSkxmnk6HTugo4VDCltsEN797n6MoSrFyh08rJUdunlmWAwcPHqRNmzY8//zzXL58mbp16/LPP//QdXxXnHyc7n0BRVGUPKSeFIAT4fFcvZ5iTf9+9CLn4g7RDc9s5+3bcRCTwzGum4xA5v4CkYbfiKz7h3Wv4rAwV5JNp6zHE9ITSDWnUsnRQlMXEydiTuB2OZTJly6xrEEDLBYLnq4uPNu8IW8OHk65xrWY9+l46l5xoZpZG2HkFa8tfRF54Tw/TRx51/vpFO5H3IHNzP/5sDYUNSZjDaPRb9D0uiP6ehHQ9H5/WoqiFGdFMijkZUezxSJ58vNtpBozVyQd+NuXPOYYTjf6Zju3/EVHqgX9xaaEb0k2tsPZTlu4LsbwN2mul3BMLoc0WRA6gatrDczmJNIsgsWXTwLwvBfUczDyzLvPkPRzFPFGI3q9npAXnqCOJZrrZjdObd3AUuctlLqkxy0B0jJWzo51T6eWRw0cnJzveU96oScmNYaY1IymI/uMBexiT+MVb4/9qdhc/tQURSmuimRQyOuO5lSjhV6Ny/JEbX+c7Q2svXiG+P3aB+qW2udZnvATtdKr0v/qU7Qu+Rink77CaDFmu4ZDYlkaP7YSS7IRO5/MD+4DEQf478RLjGw8kmOb5jH94z2EXdCeSpp5eTF/yxZqpu/np9lzuZ4MSAupplQcgMCK1fm+grY66R9OdSn13Gyb7qdN0jUiUyIzM+JCwakE2Luwa+Z8nA2qWUpRlNsrkkEhP/i7O9KysjZha891R24M5vSq5cXpA5fw1mudsxW8AiDp9tfQu9ihd7G7JT89Mp35Q+ez5fetANj52DGnXms6uLpToWZNyFjqOisngxMGXebGPKU8bN+7wM/FDz+XLJ3EJWtZXx61u/eThqIoDy8VFPJRYmIin3/0OafnnOaU6RROTgZ6vFiGPY2ceHR3EDJaNeMoilK4qKCQF6S2ioQ1KSVLlixhxIgRXLmizU149KlHGdzfCYtLHPvORRVQRRVFUe5OBYU8tm/fPkJCQti2bRsA1etWJ61LGuPfGI/DtTmEJ8Xd4wqKoigF56ENCuYkI9JoISXdSA3iKBkfjumyHVJKksLDcDNnbFKTqn2Iy4wlpu0TUilpLEHEmUMkZ7TPm03pxF5P5YVnnmDpyt+RUlLSy5ORA1+jXtvKvH9oDuLiaaQlAV1aOt7xkrSEaIQxlfCw/RB3iXS0R410YzKpqSZt7oNad0hRlAdMyFxMvCpowcHBcs+ePTl+X3poAhFzDth07okS/Rjib0eN5Ap8cnHYLccj/Tcy/9RnfL8onsRkEwbghRIleNO7JG56fbZzrwy3IAPMpEdqQ0Rjz7hyMMod0OYh2Jt0pBssxLinUzLejiD7BCa30YLS4cAe0P69HN/rzX6aOJLIC+fxCSp/39eo3rwNdTp0ynVdFEUpGEKIvVLK4NsdK5JPCrmdp2BO0oaTurUtw65/v6eJPIS51nPYG3Rsjj/LpoQzvOwdTBVPT1r4D+ejuNOkm81In3O4WJw5nhqOWUp2HT7OZ58s4lKoNsS0hYsLr73agrI+7iQAN1ZRchAGqjv7cTziAAbHKHR2YO+ViodDKkE7MkYV+eqxGM04GsAdF5z97KleqzF/VvMl2WKC2r1z90PLUL15m1y9P/LCeQAVFBSlmCqSQSGv5ik4VvcifOdJHOU/2Pf8GYDLB+ey4cB2pvX+An3GkNAuN73P78wZhg4dypo1awAoU8aVD9o/S/C+/VT/7Pc7lmecOBLjGXh+/BT27usNvtC155Lc3EKO1enQKVcf6PeaTa0oStGm1j7KgYSEBEaMGEGNGjVYs2YNrq6uDAypy08/P8GjNWogsg5BUhRFKYJUULCBxWJh0aJFVKlShWnTpmE0GnnllVc4deoUr7xSA3t7/b0voiiKUgQUyeajB2nXrl2EhISwc+dOABo3bszs2bNp0qQJAFfDCrJ2iqIoeUs9KdxBeHg4r776Kk2aNGHnzp34+/vz7bff8t9//1kDgqIoSnGjnhRuYjKaiFwXSfWB1UlISMDOzo4hQ4YwduxY3NzcCrp6iqIo+UoFhSx+++033n/zfSIuRQDQpUsXZsyYQeXKlQu4ZoqiKA9GoQkKQggX4AsgHdgspfzhQZV9MsrEiCeeYP369QDY+9uzYsEKOj/R+UFVQVEUpVDI1z4FIcQ3QogIIcSRm/I7CSFOCiHOCCFuDHx/GlgupewHdMvPet0QFx/Hkt/30uCLGNavX4+7uzvPvvMslSdXplMnNTlLUZSHT353NH8LZPt0FULogc+Bx4EaQC8hRA0gELiccZo5PytlsVj48eBaardtyO//Hcdsgb59+3L69Gnav9QeYVDzDRRFeTjla/ORlHKLECLopuzGwBkp5TkAIcSPQHcgFC0wHCAfg5XJZKLcU02IPnkOALcaVaj4Zh9OVq7AMwf+BkrgGTCWNn/9QqO0XbRI+xchJXD7NaKc3eNIjvdg+3E9rmVKcvAuM35vXnMoMfGYNrM5D/j7daV06V55cq17sWWfaLU+kqIUTQXRp1CazCcC0IJBE2A2MEcI0RlYc6c3CyH6A/0BypYtm+PCDQYD3lXKkRAVS9DrL+DTtnnGTOTsTweh+kBwgBZp/yKF4E7PDskJnkREViDJx58SdnePZT5B5a1rD/n7dSU8x7W/vcTEY4TDAwkKtqydpNZHUpSiq9B0NEspk4BXbThvPjAftFVS76esPT+sRgiBq6vrHc95av9pwAdx+Ast/U6DO567cOFC/KrBi69Ot7kOpUv3yrMP8bx62rCFLWsnqfWRFKXoKoigcAUokyUdmJFns9yukqrmGyiKotxeQcxo3g1UFkKUF0LYAz2BX3NyASnlGillfw8Pj3ypoKIoysMqv4ekLgX+A6oKIUKFEH2llCbgbeAP4Djws5TyaA6v21UIMT8uTm1tqSiKkpfye/TRbRvNpZTrgHW5uG6e7KegKIqiZFckF8RTTwqKoij5o0gGBdWnoCiKkj+KZFBQFEVR8keRDAqq+UhRFCV/CCnva/5XoSCEiAQu3ufbSwJReVidokDd88NB3fPDITf3XE5K6XO7A0U6KOSGEGKPlDK4oOvxIKl7fjioe3445Nc9F8nmI0VRFCV/qKCgKIqiWD3MQWF+QVegAKh7fjioe3445Ms9P7R9CoqiKMqtHuYnBUVRFOUmD2VQuMMe0cWWEKKMEOJvIcQxIcRRIcSggq7TgyCE0Ash9gsh1hZ0XR4EIYSnEGK5EOKEEOK4EKJpQdcpvwkhhmT8Th8RQiwVQjgWdJ3y2u32uhdCeAkh/hRCnM74f4m8Ku+hCwp32SO6ODMB70gpawCPAG89BPcMMAhtJd6HxSzgdyllNaAuxfzehRClgRAgWEpZC9CjLcVf3HzLTXvdAyOBjVLKysDGjHSeeOiCAln2iJZSpgM39ogutqSUYVLKfRmvE9A+LEoXbK3ylxAiEOgMfF3QdXkQhBAeQCtgAYCUMl1Keb1AK/VgGAAnIYQBcAauFnB98pyUcgsQc1N2d2BRxutFwJN5Vd7DGBRut0d0sf6AzEoIEQTUB3YWcFXy20zgXcBSwPV4UMoDkcDCjCazr4UQLgVdqfwkpbwCfAxcAsKAOCnlhoKt1QPjJ6UMy3gdDvjl1YUfxqDw0BJCuAK/AIOllPEFXZ/8IoToAkRIKfcWdF0eIAPQAPhSSlkfSCIPmxQKo4x29O5oAbEU4CKEeLFga/XgSW0IaZ4NI30Yg0Ku94guioQQdmgB4Qcp5YqCrk8+aw50E0JcQGsebCeE+L5gq5TvQoFQKeWNJ8DlaEGiOOsAnJdSRkopjcAKoFkB1+lBuSaECADI+H9EXl34YQwKud4juqgRQgi0tubjUsoZBV2f/CalHCWlDJRSBqH9+26SUhbrb5BSynDgshCiakZWe+BYAVbpQbgEPCKEcM74HW9PMe9cz+JX4JWM168Aq/Pqwvm6HWdhJKU0CSFu7BGtB77J6R7RRVBz4CXgsBDiQEbe6IxtUZXiYyDwQ8aXnXPAqwVcn3wlpdwphFgO7EMbYbefYjizOWOv+zZASSFEKDAemAL8LIToi7ZS9HN5Vp6a0awoiqLc8DA2HymKoih3oIKCoiiKYqWCgqIoimKlgoKiKIpipYKCoiiKYqWCgqLYKGMV0jczXpfKGA6pKMWKGpKqKDbKWDdqbcaKnIpSLD10k9cUJRemABUzJgCeBqpLKWsJIfqgrVLpAlRGW6TNHm3CYBrwhJQyRghREW3Zdh8gGegnpTzxoG9CUe5GNR8piu1GAmellPWA4TcdqwU8DTQCPgCSMxam+w94OeOc+cBAKWVDYBjwxYOotKLkhHpSUJS88XfGXhUJQog4YE1G/mGgTsYKtc2AZdoyPQA4PPhqKsrdqaCgKHkjLctrS5a0Be3vTAdcz3jKUJRCSzUfKYrtEgC3+3ljxv4V54UQPUBbuVYIUTcvK6coeUEFBUWxkZQyGtiWsYH69Pu4xAtAXyHEQeAoxXwbWKVoUkNSFUVRFCv1pKAoiqJYqaCgKIqiWKmgoCiKolipoKAoiqJYqaCgKIqiWKmgoCiKolipoKAoiqJYqaCgKIqiWP0fVyCfxWDtJBoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "for bd1 in bd_list:\n",
    "    plt.plot(bd1.time_hist, bd1.pop_hist, drawstyle='steps-post')    # stochastic realizations\n",
    "plt.plot(t_array, N_analytic, 'k', linewidth=2, label='deterministic')    # deterministic solution from above\n",
    "plt.yscale('log')\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('population size')\n",
    "plt.legend(loc='upper left')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hopefully you see a few trajectories that did not go extinct. For those trajectories, you should see that asymptotically they grow exponentially at the same rate as the deterministic solution we found above. However, you may also see a few trajectories that went extinct, usually pretty early on. These are due to stochastic fluctuations that are significant when $N$ is small. Once the population size has reached a certain level, say $N > 100$, extinction becomes extremely unlikely."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Extinction probability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The stochastic simulation allows us to estimate the probability that extinction happens --- something that could not be calculated from the deterministic solution above. To do that, we need to run the simulation many times to collect the statistics. To save some running time, we will run each simulation only until $T=5$, by which time most extinctions would have happened."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Extinction happened 518 out of 1000 times.\n",
      "estimated extinction probability = 0.518\n"
     ]
    }
   ],
   "source": [
    "trials = 1000    # number of simulations to run\n",
    "T = 5.    # total time for running each simulation\n",
    "extinct = 0    # number of extinctions among all trials\n",
    "\n",
    "for i in range(trials):\n",
    "    print(f'running simulation #{i+1}/{trials}', end='\\r')    # print progress\n",
    "    bd1 = BirthDeath(beta, delta)\n",
    "    bd1.run(T)\n",
    "    if bd1.pop_size == 0:\n",
    "        extinct += 1\n",
    "\n",
    "PE = extinct / trials\n",
    "print(f'Extinction happened {extinct} out of {trials} times.')\n",
    "print(f'estimated extinction probability = {PE}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mathematically, one can show that the probability of extinction in such a birth-death process is given by $q = \\delta / \\beta$. In our example this should be $0.5$. You can see that the estimate from our simulations is pretty close!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let us examine how the extinction probability depends on the initial population size. In the example above, we assumed $N_0 = 1$, i.e., the population is founded by one individual alone. With more founders, we may expect the extinction probability to be smaller, because it takes a larger and rarer fluctuation for all of them to die."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us run simulations with different initial sizes of the population."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N0 = 1: estimated extinction probability = 0.488\n",
      "N0 = 2: estimated extinction probability = 0.23\n",
      "N0 = 3: estimated extinction probability = 0.135\n",
      "N0 = 4: estimated extinction probability = 0.056\n",
      "N0 = 5: estimated extinction probability = 0.034\n",
      "N0 = 6: estimated extinction probability = 0.013\n",
      "N0 = 7: estimated extinction probability = 0.008\n",
      "N0 = 8: estimated extinction probability = 0.004\n"
     ]
    }
   ],
   "source": [
    "N0_list = np.arange(1, 9)    # list of initial sizes to use\n",
    "PE_list = []    # list of extinction probability for each initial size\n",
    "\n",
    "trials = 1000    # number of simulations to run\n",
    "T = 5.    # total time for running each simulation\n",
    "\n",
    "for N0 in N0_list:\n",
    "    extinct = 0    # number of extinctions among all trials\n",
    "    for i in range(trials):\n",
    "        print(f'running simulation #{i+1}/{trials}', end='\\r')    # print progress\n",
    "        bd1 = BirthDeath(beta, delta, N0=N0)    # specify initial population size N0\n",
    "        bd1.run(T)\n",
    "        if bd1.pop_size == 0:\n",
    "            extinct += 1\n",
    "    PE = extinct / trials    # estimate extinction probability\n",
    "    PE_list.append(PE)    # save extinction probability for given initial size\n",
    "    print(f'N0 = {N0}: estimated extinction probability = {PE}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAYrklEQVR4nO3df3RlZX3v8feHMMgBYQIy1+UExwGxQSo/IhFroVZFbrBlYEpR/FmrLLloS5HexkW8rKpdVWYZ9SpV1AEVETouwHEUscS2I/5oaWF+AAGH3Gv5IZNBx18ZoetgQ/j2j73PeJJJdnaSs7PPOfm81srKOc/e55xPhmG+2c/z7OdRRGBmZjaT/coOYGZmzc2FwszMMrlQmJlZJhcKMzPL5EJhZmaZ9i87QBGOOOKIWL16ddkxzMxaxtatW38WESumO9aWhWL16tVs2bKl7BhmZi1D0iMzHWurridJaySt37NnT9lRzMzaRlsVioi4JSIuXL58edlRzMzaRlsVCjMzazwXCjMzy+RCYWZmmdpy1tN8bNo+yuDQCLvGqqzsrNDf183anq6yY5mZlc6FgqRIDGwcpjo+AcDoWJWBjcMALhZmtuS1VdfTfKfHDg6N7C0SNdXxCQaHRhoZz8ysJbVVoZjv9NhdY9U5tZuZLSVtVSjma2VnZU7tZmZLiQsF0N/XTWVZx6S2yrIO+vu6S0pkZtY8PJjNbwasPevJzGxfLhSptT1dLgxmZtNw15OZmWVyoTAzs0wuFGZmlqmtCoX3ozAza7y2KhTej8LMrPHaqlCYmVnjuVCYmVkmFwozM8vkQmFmZplcKMzMLJMLhZmZZXKhMDOzTC4UZmaWyYXCzMwyuVCYmVmmtioUXuvJzKzx2qpQeK0nM7PGa6tCYWZmjedCYWZmmVwozMwskwuFmZllcqEwM7NMLhRmZpZp/7ID2Nxt2j7K4NAIu8aqrOys0N/XzdqerrJjmVmbcqFoMZu2jzKwcZjq+AQAo2NVBjYOA7hYmFkh3PXUYgaHRvYWiZrq+ASDQyMlJTKzdudC0WJ2jVXn1G5mtlAuFC1mZWdlTu1mZgvlQtFi+vu6qSzrmNRWWdZBf193SYnMrN15MLvF1AasPevJzBaLC0ULWtvT5cJgZovGXU9mZpaprQqFNy4yM2u8tioU3rjIzKzx2qpQmJlZ47lQmJlZJhcKMzPL5EJhZmaZXCjMzCyTC4WZmWVyoTAzs0wuFGZmlsmFwszMMs1aKCRdLOmwxQhjZmbNJ88VxbOBuyTdKOlMSSo6lJmZNY9ZC0VEXA68APgc8KfA/5f0IUnPLzibtYlN20c5dd1mjrrsVk5dt5lN20fLjmRmc5BrjCIiAvhx+vUUcBhws6QPF5jN2sCm7aMMbBxmdKxKAKNjVQY2DrtYmLWQPGMUl0jaCnwY+Bfg+Ih4J3Ay8McF57MWNzg0QnV8YlJbdXyCwaGRkhKZ2Vzl2eHucODciHikvjEinpZ0VjGxrF3sGqvOqd3Mmk+erqejpxYJSV8CiIgdhaSytrGyszKndjNrPnkKxW/XP5HUQdLtZDar/r5uKss6JrVVlnXQ39ddUiIzm6sZu54kDQDvBSqSflVrBv4LWL8I2awNrO3pApKxil1jVVZ2Vujv697bbmbNT8mEpowTpCsiYmCR8jREb29vbNmypewYZmYtQ9LWiOid7ljWFcWxEfEAcJOkF089HhHbGpjRzMyaVNasp/8NvAP46DTHAnhVIYnMzKypzFgoIuId6fdXLl4cMzNrNlldT+dmvTAiNjY+jpmZNZusrqc1GccCWJRCIelo4P8AyyPivMX4TDMz+42srqe3LfTNJX0eOAvYHREvqms/E/gE0AFcExHrMnI8CFwg6eaF5jEzs7nL6np6c0RcL+kvpzseER/L8f7XAp8Erqt73w7gU8AZwE6SJcy/TlI0rpjy+rdHxO4cn2NmZgXJ6no6OP1+yHzfPCK+K2n1lOZTgB+mVwpI+jJwTkRcQXL1MS+SLgQuBFi1atV838bMzKbI6nr6bPr9Aw3+zC7g0brnO4GXznSypGcBHwR6JA2kBWUfEbGe9I7x3t7e7LsIzcwst1lXj00Hkz8B/A7JIPYdwKW1K4KiRcTPgYsW47PMzGxfeRYF/HvgRuA5wErgJmDDAj5zFHhu3fMj0zYzM2tCeQrFQRHxpYh4Kv26HjhwAZ95F/ACSUdJOgB4PfD1BbyfmZkVaMZCIelwSYcD/yDpMkmrJT1P0nuAb+Z5c0kbSLqquiXtlHRBRDwF/DkwBOwAboyI+xf+o4CkNZLW79mzpxFvZ2ZmZKweK+khkjEJTXM4IuLoIoMthFePNTObm3mtHhsRRxUXyczMWkWePbOR9CLgOOrGJiLiuplfYWZm7SLP9Nj3Aa8gKRTfBF4DfJ+6u63NzKx95Zn1dB5wOvDjdP2nE4HlhaaaJw9mm5k1Xp5CUY2Ip4GnJB0K7GbyfRBNIyJuiYgLly9vyjpmZtaS8oxRbJHUCVwNbAWeIJnyamZmS8CshSIi3pU+/Iyk24BDI+LeYmOZmVmzyDvr6VzgNJL7Kr4PuFCYmS0Rs45RSLqKZFG+YeA+4H9J+lTRwebDg9lmZo03453Ze0+QHgBeGOmJkvYD7o+IFy5CvnnxndlmZnOTdWd2nllPPwTqdwJ6btpmZmZLQNZWqLeQjEkcAuyQdGd66BTgzpleZ2Zm7SVrMPsji5bCzMyaVtaigN+pPZb0bOAl6dM7I2J30cHMzKw55Jn19DqSrqbXAq8D/l3SeUUHmw/PerKF2rR9lFPXbeaoy27l1HWb2bTdmy+a5Zn1dA9wRu0qQtIK4J8i4sRFyDcvnvVk87Fp+ygDG4epjk/sbass6+CKc49nbU9XicnMirfQWU/7Telq+nnO15m1lMGhkUlFAqA6PsHg0EhJicyaQ547s2+TNARsSJ+fT86tUM1aya6x6pzazZaKzEIhScCVJAPZp6XN6yPiq0UHM1tsKzsrjE5TFFZ2VkpIY9Y8MgtFRISkb0bE8cDGRcpkVor+vu5pxyj6+7pLTGVWvjxjDdskvWT208xa29qeLq4493i6OisI6OqseCDbjHxjFC8F3izpYeA/AZFcbJxQZDCzMqzt6XJhMJsiT6HoKzxFg0haA6w55phjyo5iZtY2Zu16iohHgGcB5wBnA89K25qOt0I1M2u8PHdm/zXwRZJicQTwBUmXFx3MzMyaQ56upzcBJ0bEkwCS1gF3A39bYC4zM2sSeWY97QIOrHv+DMAL4JiZLRF5rij2APdL+keS/SnOAO6UdCVARPxFgfnMzKxkeQrFV9OvmtuLiWJmZs1o1kIREV9cjCBmZtacvAqsmZllaqtC4Y2LzMwar60KhW+4MzNrvFnHKCT9FtAPPK/+/Ih4VYG5zMysSeSZ9XQT8BngamBilnPNzKzN5CkUT0XEpwtPYmZztmn7KINDI+waq7Kys0J/X7dXv7WGy1MobpH0LpJ7KX5da4yIXxSWysxmtWn76KSNlkbHqgxsHAZwsbCGylMo3pp+769rC+Doxscxs7wGh0Ym7cYHUB2fYHBoxIXCGirPDXdHLUYQM5ubXdPs753VbjZfeWY9LQPeCbw8bbod+GxEjBeYy8xmsbKzwug0RWFlZ6WENNbO8txH8WngZOCq9OvktM3MStTf101lWcektsqyDvr7uktKZO0qzxjFSyLixLrnmyXdU1QgM8unNg7hWU9WtDyFYkLS8yPiPwAkHY3vpzBrCmt7ulwYrHB5CkU/8G1JDwIiuUP7bYWmmidJa4A1xxxzTNlRzMzahiJi9pOkZwC1js+RiPh11vll6+3tjS1btpQdw8ysZUjaGhG90x2b8YpC0qsiYrOkc6ccOkYSEbGxoSnNzKwpZXU9/T6wGVgzzbEAXCjMzJaAGQtFRLwvffg3EfFQ/TFJvgnPzGyJyHMfxVemabu50UHMzKw5ZY1RHAv8NrB8yjjFocCBRQczM7PmkDVG0Q2cBXQyeZziceAdBWYyM7MmkjVG8TXga5JeFhF3LGImMzNrInnGKC6S1Fl7IukwSZ8vLpKZmTWTPIXihIgYqz2JiF8CPYUlMjOzppKnUOwn6bDaE0mHk2/pDzMzawN5/sH/KHCHpJtI1no6D/hgoanMzKxp5Nnh7jpJW4FXpk3nRsQPio1lZmbNIm8X0gPAL2vnS1oVET8qLJWZmTWNPFuhXgy8D/gJyT4UIlnr6YRio5mZWTPIc0VxCdAdET8vOoyZmTWfPLOeHgX2FB2kESStkbR+z56WiGtm1hLyXFE8CNwu6VZg74ZFEfGxwlLNU0TcAtzS29vrJUbMzBokT6H4Ufp1QPplZmZLSJ7psR9YjCBmZtacspYZ/3hEvFvSLSSznCaJiLMLTWZmZk0h64riS+n3jyxGEDMza05Zy4xvTR+eFBGfqD8m6RLgO0UGMzOz5pBneuxbp2n70wbnMDOzJpU1RvEG4I3AUZK+XnfoEOAXRQczM7PmkDVG8a/AY8ARJCvI1jwO3FtkKDMzax5ZYxSPAI9IumDqarGSXgHcXmgyMzNrCnnGKG6U9B4lKpL+Drii6GBmZtYc8hSKlwKrSLqi7gJ2AacWGcrMzJpHnkIxDlSBCnAg8FBEPF1oKjMzaxp51nq6C/ga8BKSge3PSPrjiHhtocnMrK1s2j7K4NAIu8aqrOys0N/XzdqerrJjWQ55CsUFEbElffwYcI6ktxSYyczazKbtowxsHKY6PgHA6FiVgY3DAC4WLSBP19NWSW+W9NeQbIMKjBQby8zayeDQyN4iUVMdn2BwyP+UtII8heIq4GXAG9LnjwOfKiyRmbWdXWPVObVbc8k16yki/gx4EiAifon3pTCzOVjZWZlTuzWXXLOeJHWQLjUuaQXgWU9mllt/XzeVZR2T2irLOujv6y4pkc1FnsHsK4GvAv9D0geB84DLC01lZm2lNmDtWU+tSRH77Em070nSscDpgIB/jogdRQdbiN7e3tiyZcvsJ5qZGQCStkZE73TH8lxREBEPAA80NJWZmbWEPGMUZma2hLlQmJlZplxdT2WStBb4Q+BQ4HMR8a1yE5mZLS2FXlFI+ryk3ZLum9J+pqQRST+UdFnWe0TEpoh4B3ARcH6Rec3MbF9FX1FcC3wSuK7WkN6T8SngDGAncFe61WoH++5z8faI2J0+vhzfEW5mtugKLRQR8V1Jq6c0nwL8MCIeBJD0ZeCciLgCOGvqe0gSsA74h4jYNtNnSboQuBBg1apVjfkBzMyslMHsLuDRuuc707aZXAy8GjhP0kUznRQR6yOiNyJ6V6xY0ZikZmbW/IPZEXElyd3hZmZWgjKuKEaB59Y9PzJtMzOzJlRGobgLeIGkoyQdALwe+HoJOczMLIeip8duAO4AuiXtlHRBRDwF/DkwBOwAboyI+xv0eWskrd+zZ08j3s7MzMi5KGCr8aKAZrZQS22P7wUvCmhmtpR4j+/JvNaTmdkU3uN7srYqFB6jMLNG8B7fk7VVoYiIWyLiwuXLl5cdxcxamPf4nqytCoWZWSN4j+/JPJhtZjaF9/iezIXCzGwaa3u6lmxhmMpdT2ZmlqmtCoVnPZmZNV5bFQrPejIza7y2KhRmZtZ4LhRmZpbJhcLMzDK5UJiZWaa2KhSe9WRm1nhtVSg868nMrPHaqlCYmVnjuVCYmVkmFwozM8vkQmFmZplcKMzMLJMLhZmZZWqrQuH7KMzMGq+tCoXvozAza7y2KhRmZtZ4LhRmZpbJhcLMzDK5UJiZWab9yw5gZmYLs2n7KINDI+waq7Kys0J/Xzdre7oa9v4uFGZmLWzT9lEGNg5THZ8AYHSsysDGYYCGFQt3PZmZtbDBoZG9RaKmOj7B4NBIwz6jrQqFb7gzs6Vm11h1Tu3z0VaFwjfcmdlSs7KzMqf2+WirQmFmttT093VTWdYxqa2yrIP+vu6GfYYHs83MWlhtwNqznszMbEZre7oaWhimcteTmZllcqEwM7NMLhRmZpbJhcLMzDK5UJiZWSZFRNkZGk7ST4FH5vnyI4CfNTBOkVopK7RW3lbKCq2Vt5WyQmvlXUjW50XEiukOtGWhWAhJWyKit+wcebRSVmitvK2UFVorbytlhdbKW1RWdz2ZmVkmFwozM8vkQrGv9WUHmINWygqtlbeVskJr5W2lrNBaeQvJ6jEKMzPL5CsKMzPL5EJhZmaZXChSkj4vabek+8rOMhtJz5X0bUk/kHS/pEvKzjQTSQdKulPSPWnWD5SdKQ9JHZK2S/pG2VmySHpY0rCkuyVtKTvPbCR1SrpZ0gOSdkh6WdmZpiOpO/0zrX39StK7y86VRdKl6f9j90naIOnAhr23xygSkl4OPAFcFxEvKjtPFknPAZ4TEdskHQJsBdZGxA9KjrYPSQIOjognJC0Dvg9cEhH/VnK0TJL+EugFDo2Is8rOMxNJDwO9EdESN4RJ+iLwvYi4RtIBwEERMVZyrEySOoBR4KURMd8beQslqYvk/63jIqIq6UbgmxFxbSPe31cUqYj4LvCLsnPkERGPRcS29PHjwA6guMXoFyAST6RPl6VfTf3biaQjgT8Erik7SzuRtBx4OfA5gIj4r2YvEqnTgf9o1iJRZ3+gIml/4CBgV6Pe2IWixUlaDfQA/15ylBml3Th3A7uBf4yIps2a+jjwHuDpknPkEcC3JG2VdGHZYWZxFPBT4Atpt941kg4uO1QOrwc2lB0iS0SMAh8BfgQ8BuyJiG816v1dKFqYpGcCXwHeHRG/KjvPTCJiIiJOAo4ETpHUtF17ks4CdkfE1rKz5HRaRLwYeA3wZ2kXarPaH3gx8OmI6AH+E7is3EjZ0u6xs4Gbys6SRdJhwDkkxXglcLCkNzfq/V0oWlTa3/8V4IaI2Fh2njzSboZvA2eWHCXLqcDZad//l4FXSbq+3EgzS3+TJCJ2A18FTik3UaadwM66K8qbSQpHM3sNsC0iflJ2kFm8GngoIn4aEePARuB3G/XmLhQtKB0g/hywIyI+VnaeLJJWSOpMH1eAM4AHSg2VISIGIuLIiFhN0uWwOSIa9ptZI0k6OJ3MQNqF8z+Bpp21FxE/Bh6V1J02nQ403QSMKd5Ak3c7pX4E/I6kg9J/H04nGbtsCBeKlKQNwB1At6Sdki4oO1OGU4G3kPy2W5u+9wdlh5rBc4BvS7oXuItkjKKpp5y2kGcD35d0D3AncGtE3FZyptlcDNyQ/n04CfhQuXFmlhbfM0h+O29q6VXazcA2YJjk3/aGLefh6bFmZpbJVxRmZpbJhcLMzDK5UJiZWSYXCjMzy+RCYWZmmVworKlJ+tcc51wj6bj08Xvn8fonZjunCJJWz7ZacXrOG+ue90q6ssBMfyPp1UW9v7UmT4+1tiLpiYh4ZtGvaYR0na5vZK1WLOkVwF818wq21v58RWFNrfbbvqRXSLq9bi+DG9I7UEnbeyWtI1k9825JN0x5/TMl/bOkben+DefM8rmr6z5nR/q5B6XHTk8XtRtWso/JM9L2hyV9OG2/U9Ixafu1ks6b+jNN83nfS/Ntk1RbfmEd8Hvpz3Rp+ufwjfQ1h0vaJOleSf8m6YS0/f1prtslPSjpL6b5vI40131p3kvrs6Z/nrWbOYclRXr8+ZJuSxch/J6kY3P/x7SW5UJhraQHeDdwHHA0yR3qe0XEZUA1Ik6KiDdNee2TwB+lC+i9EvhordBk6AauiogXAr8C3qVkM5hrgfMj4niShe7eWfeaPWn7J0lWoc1rN3BGmu98oNa9dBnJ/g0nRcT/nfKaDwDbI+IE4L3AdXXHjgX6SNZ+el+6Nli9k4CuiHhRmvcL9QcjYkv6mScBt5GsTArJ3b4XR8TJwF8BV83hZ7QW5UJhreTOiNgZEU8DdwOr5/BaAR9Kl474J5L9O549y2sejYh/SR9fD5xGUjweioj/l7Z/kWSPhZoNdd/nsnvbMuBqScMkK5Uel+M1pwFfAoiIzcCzJB2aHrs1In6dbmi0m31/1geBoyX9naQzSQrhPiSdT7Jw32VKViv+XeAmJcvGf5ZkiRZrc/uXHcBsDn5d93iCuf39fROwAjg5IsaVrA4721aRUwfw8gzoxTSPnyL9pUzSfsAB07zuUuAnwInpuU/m+KwsmX9WEfFLSSeSXHVcBLwOeHv9OUqWg38/8PKImEizj6VXGbaE+IrC2s34NN0sAMtJ9pkYl/RK4Hk53muVfrOn8xtJtpocAVbXxh9IFmf8Tt1rzq/7fkf6+GHg5PTx2SRXD9Pleyy9WnoL0JG2Pw4cMkO+75EUwNqg98/y7ksi6Qhgv4j4CnA5U5b7VrLi7wbgTyLipwDpez8k6bXpOUqLjbU5FwprN+uBe2uD2XVuAHrTrp0/Id9S5yMkmwHtAA4j2XDnSeBtJN0vwyS74H2m7jWHpd1bl5BcJQBcDfy+klVeX0ayYc9UVwFvTc85tu6ce4EJSffUBpzrvB84Of28dcBbc/xMNV3A7WkX0vXAwJTj55AU06trg9pp+5uAC9Kc96fnWZvz9FizaeSZujrNax4GetNxAbO24SsKMzPL5CsKMzPL5CsKMzPL5EJhZmaZXCjMzCyTC4WZmWVyoTAzs0z/DfXWqjFoJtbPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(N0_list, PE_list, 'o', label='data')\n",
    "plt.yscale('log')\n",
    "plt.xlabel('initial population size')\n",
    "plt.ylabel('extinction probability')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the figure it looks like the extinction probability decreases exponentially with initial population size. Indeed, it can be mathematically shown that the extinction probability is simply $Q = q^{N_0} = (\\delta/\\beta)^{N_0}$. This can be intuitively understood, because in our model the individuals are all independent of each other. Therefore an initial population of $N_0$ individuals can be thought of as $N_0$ populations of one individual each. Since each of these populations has a probability $q$ of going extinct, the probability that they all go extinct is $q^{N_0}$. This relation fits our data pretty well, as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXhTZfrG8e/T0kLZBXEBVBC1DLLKLrKoYEUQWVRE3HcFxnGpgvsuDurMbwRxUEBFFB1EkE1QkEUEWQQFRVQUtQUFlbKMBVp4f3+cwBRs09AmPUl6f64rV5uTnORurTw55z3v85pzDhERkYIk+B1ARESimwqFiIgEpUIhIiJBqVCIiEhQKhQiIhJUGb8DRMKRRx7p6tSp43cMEZGYsWLFil+dczXyeywuC0WdOnVYvny53zFERGKGmf1Q0GM69SQiIkGpUIiISFBxVSjM7HwzG7Vt2za/o4iIxI24GqNwzk0FprZo0eJ6v7OISPTIyckhIyODXbt2+R3Fd+XKlaN27dokJSWFvE9cFQoRkfxkZGRQqVIl6tSpg5n5Hcc3zjl+++03MjIyqFu3bsj7qVAETF6ZybBZ69iYlU3Nqimkp6XSs1ktv2OJSBjs2rWr1BcJADOjevXqbNmy5bD2U6HAKxJDJq0mO2cvAJlZ2QyZtBpAxUIkTpT2IrFfUX4PcTWYXVTDZq0jO2cvCew7sC07Zy/DZq3zMZWISHRQoQA2ZmUDjlFJz/BYmdFUZmee7SIi4ffQQw/x9NNPF/j45MmT+fLLL0swUcHiqlAU9fLYmlVTKMNefnDH0C9xLnPK3kmvhIXUrFIuQklFJJpNXplJu6FzqTt4Ou2GzmXyysySz6BCERnOuanOuRuqVKlyWPulp6WSlFSWR3Mv5/w9j5PhjuIfySOZXGkobNHpJ5HSZP+YZWZWNo7/jVmGo1g8/vjjpKam0rlzZ9at8/5tefHFF2nZsiVNmjShT58+/PHHH3z88ce8++67pKen07RpU9avX5/v80pKXBWKourZrBZP9m5EraoprHV1GJTyFCubPESNnetgZDv44GHYU3L/UUTEP/vHLPMKx5jlihUrmDBhAitXrmTSpEksW7YMgN69e7Ns2TI+++wz/vKXvzB69GhOP/10evTowbBhw1i1ahX16tXL93klRVc9BfRsVuuQK5w6w87L4P374aNnYc1EOO9pOCXNt4wiEnkFjU0Wd8xy4cKF9OrVi/LlywPQo0cPANasWcN9991HVlYWO3fuJC0t/39jQn1eJOiIIpiKNaDXC3DVdCiTAq9fDBP6w7YMv5OJSITUrJpyWNsPR36Xpl511VUMHz6c1atX8+CDDxY4ezzU50WCCkUo6pwBN30EZz8I386B4a1g0b9gb47fyUQkzNLTUklJSjxoW0pSIulpqcV63Q4dOvDOO++QnZ3Njh07mDp1KgA7duzg2GOPJScnh/Hjxx94fqVKldixY8eB+wU9rySoUISqTDK0vx0GfAJ1O3inpP7dEX5c4ncyEQmjvGOWBtSqmsKTvRsVe/LtaaedRt++fWnatCl9+vShffv2ADz66KO0bt2aLl26UL9+/QPPv+SSSxg2bBjNmjVj/fr1BT6vJJhzrkTfsCS0aNHCRXzhoq+mw4y7YHsGNLscujwC5atF9j1FpEjWrl3LX/7yF79jRI38fh9mtsI51yK/5+uIoqjqd4OBS6HdrfDZG/Bcc/h0HOzbV/i+IiIxJK4KRYmvR5FcwTuSuHEh1EiFdwfC2K7wyxcl8/4iIiUgrgpFUSfcFdvRDeCqGXDBCPj1a3ihPcy+H3bvLNkcIiIREFeFwlcJCdDsMhi0Apr1h4//BSNaw9ppEIfjQCJSeqhQhFv5atDjObhmFpSrAm/2hzcuga0/+J1MRKRIVCgi5fg2cON8OOcx+H6hd3Sx8BnI3eN3MhGRw6JCEUmJSXD6IO/qqJM7w5xH4IUzYMNHficTkRKUlZXF888/D8C8efPo3r27z4kOjwpFSahSG/q+Bpf+B3J3wcvd4J2bYOfhLUcoIrEpb6GIlNzc3Ii9tgpFSTrlHLhlCbS/A1ZPhOEtYPkYzb0QiXODBw9m/fr1NG3alPT0dHbu3MmFF15I/fr16d+/P/snPq9YsYKOHTvSvHlz0tLS2LRpEwCrVq2iTZs2NG7cmF69erF161YAOnXqxD333EPHjh15/PHHqVu3Ljk5Xmuh7du3U6dOnQP3i0PdY0tacnk4+wFo3Bem3wHTboOV46H7P+DYxn6nE4l/MwfDz6vD+5rHNIKuQwt8eOjQoaxZs4ZVq1Yxb948LrjgAr744gtq1qxJu3btWLRoEa1bt2bQoEFMmTKFGjVq8Oabb3LvvfcyZswYrrjiCp577jk6duzIAw88wMMPP8w///lPwDtamT9/PgAbNmxg+vTp9OzZkwkTJtCnTx+SkpKK/ePpiMIvNVLhyqnQaxRk/QCjOsJ7Q2D3jsL3FZGY1qpVK2rXrk1CQgJNmzZlw4YNrFu3jjVr1tClSxeaNm3KY489RkZGBtu2bSMrK4uOHTsCcOWVV7JgwYIDr9W3b98D31933XWMHTsWgLFjx3L11VeHJa+OKPxkBk36eqek5jwKS0bCF+/AuU9Cg57e4yISXkE++ZeUsmXLHvg+MTGR3NxcnHOceuqpLF68+KDnFtZpokKFCge+b9euHRs2bGD+/Pns3buXhg0bhiVvXB1RlHgLj3BJOQK6PwvXzYGKR8F/roLX+sDv3/mdTETC4NCW4flJTU1ly5YtBwpFTk4OX3zxBVWqVOGII45g4cKFAIwbN+7A0UV+rrjiCvr16xe2owmIs0LhWwuPcKndHK7/ELr+HX5aCiPawLynIHe338lEpBiqV69Ou3btaNiwIenp6fk+Jzk5mYkTJ3L33XfTpEkTmjZtyscffwzAK6+8Qnp6Oo0bN2bVqlU88MADBb5X//792bp1K/369QtbfrUZj1bbN8Gse+CLSVCtHnR7BuqdCXiLvw+btY6NWdnUrJpCelpqsXvli8Sz0tRmfOLEiUyZMoVx48YV+JzDbTOuMYpoVflYuGgsnHa5d3XUuJ7QsA8za/2VITN+PrD4e2ZWNkMmeVdwqFiIlG6DBg1i5syZzJgxI6yvG1ennuJSvbPg5sXQaQisnUb7WV25aN8MEvjf3IvsnL0Mm7XOx5AiEg2ee+45vv32W0455ZSwvq4KRSxIKgedBsMti/l0bz0eSXqFycn309jWH3jKxqxsHwOKRL94PM1eFEX5PahQxJLq9RhS/mEG7PkrR9tWJic/wCNlxlKZ/1Kzaorf6USiVrly5fjtt99KfbFwzvHbb79Rrly5w9pPYxQxJv3c+gyZlMP83Y25o8x/uCJxNl0Tl/JDg3vAnam5FyL5qF27NhkZGWzZov5q5cqVo3bt2oe1j656ikF5r3rqVHkjz5R/hWpZq6FuB+j2LBx5st8RRSTGBLvqSYUiHuzbCytehjkPQ042tLvVazyYpNNRIhKaYIVCYxTxICERWl4LA5fDqb1gwTB4vg18877fyUQkDqhQxJOKR0HvUV6zwcRkGH8hvHk5bMv0O5mIxDAVinhUtwPctAjOuh++mQ0jWsHiEbA3cgubiEj8UqGIV2WSocOdMOATOOF0rx3IqE5eDykRkcOgQhHvjqgDl74FF4+D7N9hdBd496/wx+9+JxORGBFXhSJm24xHmhk06AEDlkLbgbDyNW8Z1pXjIQ6vehOR8IqrQhHzbcYjrWxFSHscblwA1U+CKbfA2PNg81q/k4lIFIurQiEhOqYhXP0e9HgOtqyFF86A9x+APf/1O5mIRCEVitIqIQFOuwIGroDGl8Ci/4MRreGr8LYnFpHYp0JR2lWoDj1HeEcYyRVhQj94ox9k/eh3MhGJEioU4jmhLdy0ELo8At/N844uPvoH7M3xO5mI+EyFQv4nMcnrEzVgqbdg0gcPeeMXGxb5nUxEfKRCIX9W9Ti4ZDz0mwB7/oCXz4PJt8B/f/U7mYj4QIVCCpbaFQYsgTNug8/fhOeae11q9+0rdFcRiR8qFBJccgXo/JDXO+roU2HqrTAmDX5e7XcyESkhKhQSmqPqw1XToecL8Pt6+HdHeO8e2L3D72QiEmEqFBI6M2jaz1v34rTLYckIGN4KvpyiViAicazQQmFmy81sgJkdURKBJAaUrwbn/x9c+wGUrw5vXQHjL4Lfv/c7mYhEQChHFJcANYFlZjbBzNLMzCKcS2LBcS3hhnmQ9iT8uNhbVW/+MMjd7XcyEQmjQguFc+5b59y9wCnA68AY4Ecze9jMqkU6oES5xDLQ9hYYuAxOORc+fAxGtoPv5vudTETCJKQxCjNrDDwDDAPeBi4EtgNzIxdNYkrlmnDxK9D/bdiXC6/2gLevgx2/+J1MRIoplDGKFcA/gGVAY+fcX51znzjnngG+i3RAiTEnd4ZbFkPHu71B7uEtYemLsG+v38lEpIjMFXK1ipmd6Jz77pBtdZ1zUTty2aJFC7d8+XK/Y8iv38L02+H7+ay1k7hr11X8XuVU0tNS6dmslt/pRCQPM1vhnGuR32OhnHqaGOI2kYMdeRKTGz3PHfv+SvV9vzI5+X6u3zmSxyd9wuSVmX6nE5EQlSnoATOrD5wKVDGz3nkeqgyUi3QwiQ/DZn9N5p42zKYxd5R5iysS3+c8ljJixjX0bPqANzdDRKJasCOKVKA7UBU4P8/tNOD6yEeTeLAxKxuAHZTnodyruGDPo/zsjuDhnGdhXE/v9JSIRLVQxijaOucWl1CesNAYRfRoN3QumYFisV8C+xhQcT53JL4Jubu8poNn3A5JOlAV8UuRxijM7K7At5ea2b8OvUUkqcSd9LRUUpISD9pWNimJet1u81qBNLgA5j/lTdb79gOfUopIMAWOUQBrA1/10VyKbP/VTcNmrWNjVjY1q6YcfNVTn5eg2WUw/Q54rQ806AnnPunNyxCRqFDoqSe/mdmJwL1AFefchaHso1NPMSh3Nyz6P1jwNCQmw1n3QsvrvZnfIhJxwU49FVgozGwqUGAVcc71COGNx+ANiG92zjXMs/1c4P+AROAl59zQEF5rogpFKfD7dzAj3TsNdUwj6P5PqJ3v366IhFGwQhHs49rTYXjvl4HhwKt5wiQCI4AuQAZes8F38YrGk4fsf41zbnMYckisqHYi9J/ozep+bzC81BmaXwWdH4QUNTAW8UOBhcI5V+yubs65BWZW55DNrYBv98/2NrMJwAXOuSfxjj6KxMxuAG4AOP7444v6MhINzODUnnDS2fDhk/DJSFg7Fc55DJpcorkXIiUs2FVPbwW+rjazz/PcVpvZ58V4z1rAT3nuZwS2FZSjupm9ADQzsyEFPc85N8o518I516JGjRrFiCdRo2wlOPcJuGE+VKsLk2+CV86HLev8TiZSqgQ79XRr4GuRP+UXIL+Pg8HGQn4DbgpzBoklxzaGa2bDylfh/Qe9NuanD4IO6ZBc3u90InGvwCMK59ymwNcfgN1AE6AxsDuwragygOPy3K8NbCzG60lpkJDgjVUMWgGNLoKPnoXnW8O69/xOJhL3Qmkzfh2wFOiNtw7FEjO7phjvuQw42czqmlky3gp67xbj9Q4ws/PNbNS2bdvC8XISjSocCb1GwlUzIKk8vNEXJvSHrJ8K31dEiiSUFh7rgNMDp4Aws+rAx8651EJf3OwNoBNwJPAL8KBzbrSZnQf8E+9KpzHOuceL9VMcQpfHlhK5e2DJCJj3FFgCdBoMbW6GxCS/k4nEnCLNo8iz8xygq3NuT+B+MjDDOdc57EnDRIWilMn6EWbeDetmwFENoPs/4Pg2fqcSiSlF7fV0u5ndDmQCn5jZQ2b2ILAEUMtPiR5Vj4d+b8Alr8PuHTAmDaYMgP/+5ncykbgQbIyiUuC2HpjM/65MmgJsinCuItEYRSlXvxsM+ATa3QqfTYDhLeDTV2HfPr+TicS0qO/1VBQ69ST88qXXaPDHj+G4NtD9WTj6VL9TiUStYi2FamY1zGyYmc0ws7n7b+GPKRJGRzeAq2fABc/Db9/AC+1h9n2we6ffyURiTihrZo8HvgLqAg8DG/AucRWJbmbQrL+37kWz/vDxczCitdcOJA6PpEUiJZRCUd05NxrIcc7Nd85dA+iSEokd5atBj+e82d0pVeHNy+D1vrB1g9/JRGJCKIUiJ/B1k5l1M7NmeLOpRWLL8a29vlHnPA4/LIIRbbz1L3L3+J1MJKqFUigeM7MqwB3AncBLwG0RTVVEuupJCpVYBk4fCAOWwsldYO6j8EI7+H6h38lEopauepLS7Zv3vaujsn6Axpd4rcwrqvuwlD7FverpRDObama/mtlmM5sSWJ5UJPad3MWbe9EhHda8DcObw7LRmnshkkcop55eB94CjgFqAv8B3ohkKJESlZQCZ90HN3/Mlor1YfrtrHqoBVc98RKTV2b6nU7Ed6EUCnPOjXPO5QZurxFk/QiRWDU5owIdfrmNW/fcQi3bwujdd7L9nTuYtlQLJUnpFqzXUzUzqwZ8aGaDzayOmZ1gZncB00suokjJGDZrHdk5+5iy7wzO3v0M4/d25jJ7jzYz0mDNJM29kFKrwMFsM/se78gh3xXpnHNRN05hZucD55900knXf/PNN37HkRhTd/D0Px0qN7b1PJ40mkYJG6DeWXDe01C9nh/xRCKqWG3GY5GuepKiaDd0LplZ2X/aflyVZBae9R3MfQxyd0P726Hd3yCpnA8pRSKjuFc9JZnZX81sYuA20My0MozEnfS0VFKSEg/alpKUyB3nNoDWN8LAZfCX7jDvSRh5OqxXyzMpHUIZzB4JNAeeD9yaB7aJxJWezWrxZO9G1KqaggG1qqbwZO9G9GxWy3tCpWPgwjFw+Tve/XG9YOI1sONn3zKLlIRQVrj7zDnXpLBt0USnniTicnbBov+Dhc9AmbLe5bUtr4OExML3FYlCxTr1BOw1swOjd4HJdnvDFU4kJiWVg053wy2LoXZLmHkXvHgmZK7wO5lI2IVSKNLxLpGdZ2bzgbl4fZ9EpHo9uOxtuOhl2LkZXjzbawmSneV3MpGwKRPsQTNLALKBk4FUvEtlv3LO7S6BbIctz+WxfkeR0sQMTu0F9c6GD5+Apf+GL6d4XWobX+w9LhLDQhmjWOyca1tCecJCYxTiq02fwbTbIXM51GkP3Z6FGqf4nUokqOKOUcw2sz5m+lgkEpJjm8C170P3f8DPn3uX0s55FHL+PEdDJBaEUihux2sEuMfMtpvZDjPbHuFcIrEtIQFaXAMDV0CjC2Hh094yrF/P9juZyGErtFA45yo55xKcc0nOucqB+5VLIpxIzKtYA3q9AFdOgzLl4PWLvKVYt6krrcSOUI4oMLPeZvasmT1jZj0jHUok7tRtDzd9BGc/AN98AMNbwsfDYW+u38lEChVKC4/ngZuA1cAa4CYzGxHpYCJxp0wytL8DBizxCsfse2FUR/jxE7+TiQQVyhFFRyDNOTfWOTcWOA/oFNFUIvHsiDrQbwL0HQ/ZW2HMOfDuIPjjd7+TieQrlEKxDjg+z/3jgM8jE6d4zOx8Mxu1bds2v6OIBGfmNRgcsBROHwQrx8PwFrDytcNa92LyykzaDZ1L3cHTaTd0rlbkk4gIZR7FfKAlsDSwqSWwGPgDwDnXI5IBi0LzKCTm/PKFN/fipyVwfFtv7sXRDYLuMnllJkMmrSY7538ddVKSEg9uZCgSomDzKILOzA54IMx5RORQR58KV8+EVePh/Qfg3+2h7QDoeDckV8h3F29FvoPbrmXn7GXYrHUqFBJWhRYK59z8kggiUuolJMBpl0P9bl6xWPR/3hKsXZ/yth1iYz6LLAXbLlJUIV0eKyIlqHw1uGA4XDMLylaGCZfC65fA1h8OelrNqin57l7QdpGiUqEQiVbHt4Eb58M5j8H3C7yZ3Qufhdw9QMEr8qWnpfqRVuKYCoVINEtM8q6KGrgUTjob5jzsjV9s+KjwFflEwiSUq57aAQ8BJ+CNaRjgnHMnRjxdEemqJ4lbX8+CGXdC1o/QpB90edRrEyJSTMW96mk0cBuwAq1sJ+KvU9K81uULn4ZF/4J1M6HzQ3Dald5guEgEhPKXtc05N9M5t9k599v+W8STiUj+kst7PaNuXgTHNIJpf/Nmd2+KynmwEgdCKRQfmtkwM2trZqftv0U8mYgEVyMVrpwKvf4Nv3/v9Y16bwjs3uF3MokzoYxRfJjPZuecOysykYouz1Ko13/zzTd+xxEpOdlbYc4jsHwsVDoGzh0KDS7QMqwSsmBjFIUWilikwWwptTKWw7TbvJX1TuoM5w2DalF73YlEkWIthWpmVQJrUSwP3J4xsyrhjykixVa7BVz/IZz7lNe+/Pm2MP/vkLvb72QSw0IZoxgD7AAuDty2A2MjGUpEiiGxDLS5CQYug9Su8OHj3rrd383zO5nEqFAKRT3n3IPOue8Ct4cBHcuKRLvKx8JFL8Nlk2DfXnj1Aph4Lez4xe9kEmNCKRTZZnbG/juBCXjqOiYSK046G25ZAh0Hw9p3vXUvPhnlFQ+REIRSKG4GRpjZBjP7ARiOtzSqiMSKpHJw5hC4eTHUOg1mpsOLZ0Hmp34nkxhQaKFwzq1yzjUBGgONnHPNnHOfRT6aiITdkSfB5ZPhwjGw42evWEy/E7Kz/E4mUazAFh5mdplz7jUzu/2Q7QA4556NcDYRiQQzaNjHu3z2wydg6Sj4cgqkPQGNLtTcC/mTYEcU+5fVqpTPrWKEc4lIpJWr4i2KdP2HUKU2TLrOG/D+VZNV5WAhdY91zi0qbFs00YQ7kcO0by+sGAsfPAK52dDub9D+dkjSIkilRbEm3AHPhbhNRGJVQiK0vA4GLYdTe8GCv8PzbeCbD/xOJlEg2BhFW+B0oMYh4xSVgcT89xKRmFbxKOg9Cpr2h+l3wPg+Xs+oc4dC5Zp+pxOfBDuiSMYbiyjDweMT24ELIx9NRHxzYkevjflZ93mLJQ1vCYtHwN5cv5OJD0IZozjBOfdD0CdFGY1RiITR79/DzLvgm9lwdCPo/g84rqXfqSTMijtG8ZKZVc3zYkeY2aywpROR6FatLlz6Flw8DrJ/h9FdYOqt8MfvfieTEhJKoTjSOXdgNo5zbitwVOQiFZ2ZnW9mo7Zt2+Z3FJH4YgYNesCAT6DtAPh0nNcKZNXrEIdLFcjBQikU+8zs+P13zOwEICr/MpxzU51zN1Spoi7oIhFRthKkPQ43LoBq9WDyzfByN9i8ttBdJ6/MpN3QudQdPJ12Q+cyeWVmCQSWcAilUNwLfGRm48xsHLAAGBLZWCIS1Y5pCNfMgvP/BZu/hBfOgPcfhD3/zffpk1dmMmTSajKzsnFAZlY2QyatVrGIEaH0enoPOA14E3gLaO6c0xiFSGmXkADNr4SBy6HxJbDonzCiDXw1409PHTZrHdk5B3erzc7Zy7BZ60oqrRRDKEcUAGWB34FtQAMz6xC5SCISUyocCT1HwNUzIbkCTOgHb1wKWT8deMrGrPxXJihou0SXAifc7WdmTwF9gS+AfYHNDu8UlIiI54TT4aaF3nyL+U/BiFbQ8W5oO4CaVVPIzKco1KyqFiGxoNBCAfQEUp1zWnRXRIJLTIIz/gYNe8PMwfDBg/DZBIa2GMIN88sedPopJSmR9LRUH8NKqEI59fQdkBTpICISR6oeD/1eh34TYM9/af/RFbxX5w1OrbIHA2pVTeHJ3o3o2ayW30klBKEcUfwBrDKzOcCBowrn3F8jlkpE4kNqV6jbAeb/nRMWD2d62flw0UPQ7ApvMFxiQiiF4t3ATUTk8CVXgC4PQ5NLvEaDU2+FleO9ViDHNPQ7nYSg0F5PsUi9nkSilHPw2Rsw+z5v+dU2N0Onwd5EPvFVkXo9mdlbga+rzezzQ2+RCisiccwMml7qzb047XJYPByGt/KWYo3DD63xosAjCjM71jm3KdCy40+iuaOsjihEYsRPS2Ha7fDLajj5HOj6d68JoZS4Ih1ROOc2Bb69xTn3Q94bcEskgopIKXNcK7hhHqQ9AT987K2qt2AY5Opq/GgSymUHXfLZ1jXcQUSklEos43WkHbAUTkmDuY/ByHbwveb0RotgYxQ3m9lqoP4h4xPfA6tLLqKIlApVasHFr0L/ibAvB145H96+HnZu9jtZqRdsjKIKcATwJDA4z0M7nHNRvWKJxihEYlxONix81ms0WCYFOj8Aza+GhES/k8Wtoo5RbHPObQBGHzJG8buZXRmpsCIiJKXAWffCzR9Dzabe/IuXOsPGVX4nK5VCGaN4wMxGmlkFMzvazKYC50c6mIgIR54MV0yB3i/Btgx48UyYcRfs0iqWJSmUQtERWA+sAj4CXnfOXRjRVCIi+5lB44tg4DJocS0sHQXDW8LqiZp7UUJCKRRHAK3xisVu4AQzs4imEhE5VEpV6PY0XD8XKh0Lb18L43rCb+v9Thb3QikUS4CZzrlzgZZATWBRRFOJiBSk1mlesTjvacj81Jt78eETkLPL72RxK5RC0dk5NwbAOZcd6Bo7uJB9REQiJyERWl3vnY76Sw9voaSRbeHbOX4ni0uhFIqfzOwyM3sAwMyOB0qsdJtZTzN70cymmNk5JfW+IhIDKh0DF46GyyeDJcBrveE/V8H2TYXuKqELpVA8D7QF+gXu7wBGhPLiZjbGzDab2ZpDtp9rZuvM7FszC3p04pyb7Jy7HrgKb0lWEZGD1TvTu5T2zHvhqxneYPeSkbA31+9kcSGUQtHaOTeAwFGEc24rkBzi678MnJt3g5kl4hWarkADoJ+ZNTCzRmY27ZDbUXl2vY8QC5SIlEJlykLHu2DAEji+Nbw32LucNkOTb4srlEKRE/jH3QGYWQ1gXygv7pxbABw6i7sV8K1z7jvn3B5gAnCBc261c677IbfN5nkKb0D905B/MhEpnaqd6LUBuehl+O8Wb6LetNsge6vfyWJWKIXiX8A7wFFm9jjeXIonivGetYCf8tzPCGwryCCgM3Chmd1U0JPM7AYzW25my7ds2VKMeCIS88zg1F5eo8E2N8OKl73TUZ9N0NyLIghphTszqw+cDRgwxzm3NuQ3MKsDTHPONQzcvwhIc85dFzkY3HIAAA1/SURBVLh/OdDKOTfosNMXQL2eROQgmz6H6bdDxjKo0x66PQM1Uv1OFVWK1OspL+fcV865Ec654YdTJAqQARyX535tYGMxX1NEpGDHNoZrZkP3f8LPq7025h88DHv+8DtZTAipUITZMuBkM6trZsnAJcC7PuQQkdIkIQFaXO0tw9roIvjoWXi+NXw9y+9kUS+ihcLM3gAWA6lmlmFm1zrncoGBwCxgLfCWc+6LML3f+WY2ats2NQwTkQJUrAG9RsJVMyCpPLx+MUzo7zUdlHyFNEYRazRGISIhyd0Di4fD/L97E/Y6DfYGvxOTmLwyk2Gz1rExK5uaVVNIT0ulZ7Ng193EtmBjFCoUIiJbf4CZd8PXM+GoU1lwyhBunJ9Eds7eA09JSUrkyd6N4rZYFHswW0Qkrh1xAlw6AS55HXZto8NHl/GgG0lVdhx4SnbOXobNWudjSP/EVaHQGIWIFEv9bjBwKS/knk+fxIXMLXsHFyXOwwJzjDdmZfsc0B9xVSicc1OdczdUqVLF7ygiEquSKzCu4jV02/ME37paDEsaxVvJj5BqP1Kzaorf6XwRV4VCRCQc0tNS+alMHfruuZ/0nBuoZxuZnnwPL9eeCrt3+h2vxKlQiIgcomezWjzZuxE1q1Zg4t5O9C83gow6vTn52zEwojWsnVaqWoHoqicRkVD9uASm3Q6bv4BTzoWuf/cGwuNAqbnqSYPZIhJRx7eBG+fDOY/B9wu9o4uFz3jzMeKYjihERIpiW4a35sXaqXBkqtdosG57v1MVWak5ohARKTFVakPf1+DS/0DuLnilO0y6EXbG3zIHKhQiIsVxyjlwyxJofyeseRuGN4flY2BfSOu7xQQVChGR4kouD2ff763bfUxjb0W90V1g02d+JwsLFQoRkXCpcQpcORV6vwhZP8CoTjBzMOza7neyYomrQqGrnkTEd2bQ+GIYuAyaXw2fvAAjWsGaSTE79yKuCoVaeIhI1Eg5Aro/C9fNgYpHwcSr4bU+8Nt6v5MdtrgqFCIiUad2c7j+Q29y3k9L4fm2MG8o5OzyO1nIVChERCItIRFa3+idjqrfDeY9CSNPh/Vz/U4WEhUKEZGSUvlYuGgsXP4O4GBcL5h4Dez42e9kQalQiIiUtHpnwc2LodMQr8Hg8Jbwyb9h397C9/WBCoWIiB+SynlrdN+yGGq3gJl3wYtnQuYKv5P9SVwVCl0eKyIxp3o9uGwSXDgWdvwCL54N0++A7Cy/kx0QV4VCl8eKSEwyg4a9vcHu1jd6LUCGt4TP34qKuRdxVShERGJaucrQ9Snvctqqx8Gk6+HVHrDla19jqVCIiESbmk3h2veh27Nev6iRp8OcRyEn25c4KhQiItEoIRFaXgsDl3unpRY+7S2U9PXsko9S4u8oIiKhq3gU9B7lNRssUxZevwjevAy2ZZZYBBUKEZFYULcD3LQIzrofvnnfG+z+eDjszY34W6tQiIjEijLJ0OFOGPAJ1GkHs++FUR3hx08i+rYqFCIiseaIOnDpW3DxOMjeCmPOgXcHwR+/R+Tt4qpQaMKdiJQaZtCgBwxYCm0HwsrxMLwFbN8Y9reKq0KhCXciUuqUrQhpj8ONC7yFkirXDPtblAn7K4qISMk7pqF3i4C4OqIQEZHwU6EQEZGgVChERCQoFQoREQlKhUJERILSVU8iIjFu8spMhs1ax8asbGpWTSE9LZWezWqF7fVVKEREYtjklZkMmbSa7Bxvve3MrGyGTFoNELZioVNPIiIxbNisdQeKxH7ZOXsZNmtd2N4jrgqFWniISGmzMSv/xYwK2l4UcVUo1MJDREqbmlVTDmt7UcRVoRARKW3S01JJSUo8aFtKUiLpaalhew8NZouIxLD9A9a66klERArUs1mtsBaGQ+nUk4iIBKVCISIiQalQiIhIUCoUIiISlAqFiIgEZc45vzOEnZltAX4o4u5HAr+GMU4kxVJWiK28sZQVYitvLGWF2MpbnKwnOOdq5PdAXBaK4jCz5c65Fn7nCEUsZYXYyhtLWSG28sZSVoitvJHKqlNPIiISlAqFiIgEpULxZ6P8DnAYYikrxFbeWMoKsZU3lrJCbOWNSFaNUYiISFA6ohARkaBUKEREJCgVigAzG2Nmm81sjd9ZCmNmx5nZh2a21sy+MLNb/c5UEDMrZ2ZLzeyzQNaH/c4UCjNLNLOVZjbN7yzBmNkGM1ttZqvMbLnfeQpjZlXNbKKZfRX4+23rd6b8mFlq4He6/7bdzP7md65gzOy2wP9ja8zsDTMrF7bX1hiFx8w6ADuBV51zDf3OE4yZHQsc65z71MwqASuAns65L32O9idmZkAF59xOM0sCPgJudc4t8TlaUGZ2O9ACqOyc6+53noKY2QaghXMuJiaEmdkrwELn3EtmlgyUd85l+Z0rGDNLBDKB1s65ok7kjSgzq4X3/1YD51y2mb0FzHDOvRyO19cRRYBzbgHwu985QuGc2+Sc+zTw/Q5gLRC5ZvTF4Dw7A3eTAreo/nRiZrWBbsBLfmeJJ2ZWGegAjAZwzu2J9iIRcDawPlqLRB5lgBQzKwOUBzaG64VVKGKcmdUBmgGf+JukYIHTOKuAzcD7zrmozRrwT+AuYJ/fQULggNlmtsLMbvA7TCFOBLYAYwOn9V4yswp+hwrBJcAbfocIxjmXCTwN/AhsArY552aH6/VVKGKYmVUE3gb+5pzb7neegjjn9jrnmgK1gVZmFrWn9sysO7DZObfC7ywhauecOw3oCgwInEKNVmWA04CRzrlmwH+Bwf5GCi5weqwH8B+/swRjZkcAFwB1gZpABTO7LFyvr0IRowLn+98GxjvnJvmdJxSB0wzzgHN9jhJMO6BH4Nz/BOAsM3vN30gFc85tDHzdDLwDtPI3UVAZQEaeI8qJeIUjmnUFPnXO/eJ3kEJ0Br53zm1xzuUAk4DTw/XiKhQxKDBAPBpY65x71u88wZhZDTOrGvg+Be8P+it/UxXMOTfEOVfbOVcH75TDXOdc2D6ZhZOZVQhczEDgFM45QNReteec+xn4ycxSA5vOBqLuAoxD9CPKTzsF/Ai0MbPygX8fzsYbuwwLFYoAM3sDWAykmlmGmV3rd6Yg2gGX433a3X/53nl+hyrAscCHZvY5sAxvjCKqLzmNIUcDH5nZZ8BSYLpz7j2fMxVmEDA+8PfQFHjC5zwFMrPyQBe8T+dRLXCUNhH4FFiN92972Np56PJYEREJSkcUIiISlAqFiIgEpUIhIiJBqVCIiEhQKhQiIhKUCoVENTP7OITnvGRmDQLf31OE/XcW9pxIMLM6hXUrDjzn0jz3W5jZvyKY6REz6xyp15fYpMtjJa6Y2U7nXMVI7xMOgT5d04J1KzazTsCd0dzBVuKfjigkqu3/tG9mncxsXp61DMYHZqAS2N7CzIbidc9cZWbjD9m/opnNMbNPA+s3XFDI+9YJvM8rZvZ54H3LBx47O9DUbrV565iUDWzfYGZPmbf+xlIzOymw/WUzu/DQnymf91sYyPepme1vvzAUaB/4mW4L/B6mBfapZmaTA/mWmFnjwPaHArnmmdl3ZvbXfN4vMZBrTeDnuC1v1sDvc/9kztVm5gKP1zOz9wJNCBeaWf2Q/2NKzFKhkFjSDPgb0ACvE2m7vA865wYD2c65ps65/ofsuwvoFWigdybwzP5CE0QqMMo51xjYDtxi3mIwLwN9nXON8Brd3Zxnn+3OuVbAcLwutKHaDHQJ5OsL7D+9NBhv/Yamzrl/HLLPw8DKQL57gFfzPFYfSMPr/fRgoDdYXk2BWs65hoGfY2zeB51zywPv2RR4D68zKXizfQc555oDdwLPH8bPKDFKhUJiyVLnXIZzbh+wCqhzGPsa8ESgdcQHeOt3HF3IPj855xYFvn8NOAOveHzvnPs6sP0VvDUW9nsjz9fDWb0tCXjRzFbjdSptEMI+ZwDjAJxzc4HqZlYl8Nh059zuwIJGm/nzz/odcKKZPWdm5+IVwj8xs4vxGvcNNq9b8enAf8xrG/9vvBYtEufK+B1A5DDszvP9Xg7v77c/UANo7pzLMa87bGFLRR46gOfwCk6o++z/PpfAh7LAUUxyPvvdBvwCNAk8d1ch70MBWfa/Z9DflXNuq5k1wTvqGABcDFxz0IubnYp31NLBObfXzBKArMBRhpQiOqKQeJOTz2kWgCp460zkmNmZwAkhvNbx9r81nfvhLTX5FVBn//gDXnPG+Xn26Zvn6+LA9xuA5oHvL8A7esgv36bA0dLlQGJg+w6gUgH5FuAVwP2D3r+Gui6JmR0JJDjn3gbu55B234EjkwnAFc65LQCB1/7ezC4KPMcCxUbinI4oJN6MAj43s08PGacYD0w1s+V4p61CaXW+FrjSzP4NfIO34M4uM7sa7/RLGbyOuC/k2aesmX2C9yGsX2Dbi8AUM1sKzMFbsOdQzwNvB/4R/jDPcz4Hcs3rEPsysDLPPg/hrRb3OfAHcGUIP9N+tQL77v+wOOSQx3viFdMX9w/lBI4k+gMjzew+vII3AfjsMN5XYpAujxXJRyiXruazzwagRWBcQCRu6NSTiIgEpSMKEREJSkcUIiISlAqFiIgEpUIhIiJBqVCIiEhQKhQiIhLU/wN9TZHorbmmxQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(N0_list, PE_list, 'o', label='data')    # plot simulation results\n",
    "plt.plot(N0_list, (delta/beta)**N0_list, label='theory')    # plot theoretical curve\n",
    "plt.yscale('log')\n",
    "plt.xlabel('initial population size')\n",
    "plt.ylabel('extinction probability')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It also explains why extinction becomes extremely unlikely when the population reaches a decent size, such as $N > 100$. In that case, the probability that the population will go extinct in future is $Q < q^{100} \\sim 10^{-30}$!"
   ]
  }
 ],
 "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": 4
}