{ "cells": [ { "cell_type": "markdown", "id": "a886faa5", "metadata": {}, "source": [ "# Modelling evolution in fluctuating environments\n", "\n", "In this notebook we will explain the method for modelling evolution when the population dynamics fluctuate developed by Ferris & Best (2018). In particular, this code contains the workings to produce the results for the two case studies given in the main manuscript.\n", "\n", "## Case study 1: seasonal forcing in competition\n", "\n", "We first look at the evolution of host avoidance to parasitism (reduced infection) when the degree of competition for resources oscillates seasonally.\n", "\n", "The population dynamics of our system are given by a standard susceptible-infected-susceptiblke framework as follows,\n", "\\begin{align}\n", "\\frac{dS}{dt} &= (a-q(t)(S+I))S-bS-\\beta SI + \\gamma I\\\\\n", "\\frac{dI}{dt} &= \\beta SI - (b+\\alpha+\\gamma) I.\n", "\\end{align}\n", "\n" ] }, { "cell_type": "markdown", "id": "f3dca7d3", "metadata": {}, "source": [ "In the first code cell below we run the resident dynamics of the system for some chosen parameter values, demonstrating the resulting yearly limit cycle." ] }, { "cell_type": "code", "execution_count": 4, "id": "492c9518", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.0, 30.0)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAELCAYAAAAP/iu7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydd3hUZfbHvyeQQkIxEBJC7yCCIAQERYo0FQXs9Sd23dUVO4q6gq66WNa2uLZVbCvqoojiKiBVAekoVVoIEEILPSH1/P4485LJZMrtE8L7eZ55Jrlz537fe+fec95z3kbMDI1Go9Fo3CQm2gXQaDQaTdVHOxuNRqPRuI52NhqNRqNxHe1sNBqNRuM62tloNBqNxnW0s9FoNBqN62hnozkpIaLmRMRENNHEd5iI5rhXKvchopt853FTtMviNkTUz3euY6NdFo19tLPRmMLPyPu/Cokoi4g+JqLTo13GygQRZRJRZrTLodFEm+rRLoDmpGUDgEm+v2sDOAfADQBGENHZzLzWZf2dAE4HcMhlHU30WAz5jfdFuyAa+2hno7HKemYe67+BiCYA+DOARwHc6KY4MxcBWO+mhia6MHMe9G9cZdBpNI2TTPS9dwv8gIhaEdEHRLTDl3bbQURvElH9IPsOJKIZRJRDRMd9+/5IRMP89gnZZkNEVxLRCt93dxLRy0RUI1ShiagBEb1ORFuIqICIdhPRJ0TUIsi+TERziCiNiD4ion1ElOfb1jWwfACaAWgWkHa8KfxlBIgoyVfubCLKJ6LlRHR5kP2aE1EJEU0JcZyWRFRKRN/4bZvjK0csEY0jom2+815DRNcGOUZbInqRiFYS0QHfdV1DRE8QUWyQ/TN9r9OI6B3f9Tzq+w3b+PZpT0RTfcc75Ls3agUcJ2SbDRF1JaLPiWiXr+w7iGgKEZ3nt08CET1CRL8T0REiOkxEG3xajSL8BBqH0ZGNxg2K/P8hol4AfgCQAGAqgK0A2gO4C8AQIurOzLm+fS/27ZMD4BsAuQAaADgbwDDfZyEholsA/BvAAQDvA8gHcBmAdiH2bwNgjk/jfwAmA2gC4Cpf2Xoy8+aAr50G4BcA+wF8BHEolwH4iYhOZ+YcAAcBjANwn+87r/p9f2WEc4gB8B2AfgCW+zQaAvgPgBn++zJzJhHNADCUiBr4tP25BQD5rkkgkwB0BTANYguuBfAfIjrIzP/z2+8yADcDmAVgJoB4AH0BPAMgA8CIIMeO85W1uq/cLQAMBzCDiIYCmA/gV8hvdB6Am3zfuznEZTkBEV0N4GMApQCmQO6nBr7jXO47Nnz7XAH5rd4BwL5yXOa7HjsjaWkchJn1S78MvwA0hzy0U4J89pbvszf9tsUB2AZxGqcH7H+lb/9/+m37CkABgPpBjl8vSDkm+m2rA+AwpB2nhd/2WgDW+PafE3DMhT69PgHbe0Gc5ncB29n3eg0A+W1/yrf9sYD9MwFkmrzGt6hrDCDGb/sAP/2b/LZf7tv2SMBxYgBsB7ALQHW/7XN8+y8AUNNve1/f9h8DjtMIQFzANgLwrm//3kHOmQF8BqCa3/Y3fNsPAPiz3/bqAFb4rnea3/Z+vv3H+m1rAOCY7xiB9xMBaOh3L5QC+CrI9U0AkBTtZ+lUe+k0msYq7YlorO/1MhH9CuBOAFsA/M1vv4sBNAXwHDOv8z8AM38JYBmAawKOXYSA6Mi3//4IZRoOcSxvM/NWv+8dAfBs4M6+tFdPAO8y87wArYWQyOpCIqoT8NVjAMawz3L5mOh7z4hQRiPc4Hsfw8ylfmX6CQGRjY+pAPZAnJQ/QwA0BvAhMxcH+d4YZj7qd/y5EEdR7hyYeSczFwZsYwD/8v07MMR5jGbmEr//P/e97/P7LnxlmwxxOh1CHEsxEkAigPFB7idm5mz1L8T55AcegJmPM/OxCDoah9FpNI1V2kFq8/5sBXAul0/lnO177xgs9w6gBoB6RJTCzPsgBulSAKuJaBKkFv4zMx80UKbOvvf5QT77Ocg2VbbGIcqWDokO2gBY6rd9YxBjpVIypxkoZyQ6A8jl4D36fgYwyH8DMxf52q4eIaJzmfkX30e3+t7fD6GzIsi2nZDKwQl8ab1bIamuMyC9D8lvl/QgxznAzFkB29R98XuAo/b/rGGIsiq6+96nh9uJmQ8T0Q8AriOiJpAocT6A5QEOUOMR2tlorPINM48AACJKA3A3gCcBfEFE5/vVpOv63kdGOF4SgH3M/DkRFQN4AMD9AB4EUExE3wK4L4gB80dFIHuCfLY7yDZVtuG+V7iy+VOhuzUzFxMRAFQLcxyj1AGwMcRnwc4DAN4D8DDEKfxCRCkALgEwj5n/CPYFZg7WbbwYFTsOvQHpZbgNkubMAVAIcayjIG04gRwOcexIn1XocBCA+o2zw+4lXAngCQDXAXjZt20fEb0C4O/+UaPGfXQaTWMbZt7NzH+F5PDPgxgghTIsg5iZwry2+R1vMjOfCyAF4gT+C4l2ppLPoodAGc/UIJ+lBdmmynZ7hLLNjXgRnOUQgp8DEPw8wMwbAcwFcBUR1QTwf5D2smAdAwzjq0j8CcAqSBvJLcw8hqXb+2d2jm0RFeFGioDAzEeZ+VFmbgqJxP8MYC8kpToq7Jc1jqOdjcZJxgA4CuBxIqrt27bY997T7MGY+QAzT2XmawH8BEkvNQ3zlVW+9/OCfNY7yDbLZTNBCcxHO6sA1CWiYO0Xwc5D8S4kCrsa0n5zGOKo7dACkjKbycyB7R/n2jy2FZb43geb+RIz/8HM//L73rBw+2ucRzsbjWP42lwmAEhGWc1xCqRH1KNE1CPwO0RUg4jO9vv/fCKKD9inOspSXsfDFGEqgCMA7vAfI+Or6T8epLy/QhzOzUR0SZCyxRJROONuhFwAKYHnFIFPfO/P+tpLVHkGIKC9JoDJkO7YTwPoCOAzloGRdlBpy17+USURtQXwmM1jW+EjAHkARlPA1EgkpPv+rk9EwTprqMiwQscBjbvoNhuN07wM4B4A9xPRa76G2ishY1gWEdF0AGsh915zSHfbhQAu8H3/H5AG+zmQnlHVIL2dOgL4hJlDtVmAmQ8S0X2Q1NEyXwcDNc5mDYL3dLoOwGxIim4+ZAxMMWTszHkQZ9He0pUQZkN6d31DRL9AetlNY+bfw3xnIiQNNgLAEt84moaQiGUagKHBvsTMBUT0McrG9rxno9zqmNlE9DUkjbmEiGb7yjIMwI+Qbteewcw5vrFUnwBY4SvbVkjasQ+A7yHn38hX3t8gHSF2QbpNXwqJNl/zstwa7Ww0DsPMe4noXwAeAnAvgL8x869E1AXAIwAuBNAfUjvdAampfuR3iOchBqwbgIsgzmIzZABoxPYHZn6fiI5CUnq3QGr6kyANxRVq+cy8mYjO8pV3OIDbIQ5hJ4BvIQMS7fA3APUgDmIQJJuQAyCks2HmUl+k9TRkoOUoyLQt10G6dgd1Nj6Us/mNmZeG2c8MIyERzqUA/gL5PR6HXB9PnQ0A+DqRbIFMizQQ0jtuDyRK/dK3WyaAsZCxSUMgkfFuSDr2BV9Uq/EQqtgD0UUxoiEARkNqmMmQxroFkEFba/32SwbwIqRmVwNS870/Qm1Qoznl8ZtBYRQzvx7t8mg0Cq+dzbWQ6TF+hTiappDaSRMAnZh5my8vPA/SMPkwZKTwY5D+/V2YeYdnBdZoTiKIqBpkkGwbAI2Z+UCUi6TRnMBTZxO0AETtICmCh5j5ZSIaDmlUPp+ZZ/v2qQPJy37CzPdGr7QaTeWDiDpBxtScB2n7Gs/Mj0a3VBpNeSpDbzQ1BYmanmQYgGzlaIATg8++RfiBdxrNqUo3yNiRsyGdAsZGtTQaTRCi4myIqBoRxflm3H0b0mCqFuI6A8DqIF9bA6CprxurRqPxwcwTfYNP6zLz7cwcrnu4RhMVotUb7VeUrXmyCZIyU1OM1IX0JAkk1/eeDBk4WA4iugPAHQCQlJTUrX17O71VNRqNRmOFZcuW7WPmCutURcvZ/B+ku2JLSJfTGUTUm5kzIaOVgzUkhZumBMz8DmTNCmRkZPDSpU71+tRoNBqNUYhoW7DtUUmjMfM6Zv6VmT+D9IOvCemVBkgEUzfI15J977qHjUaj0ZxkRL2DgG/q+E0AWvs2rYG02wTSAUCW//obGo1Gozk5iLqz8c0q2x4yKhmQ+a0aEVFfv31qQ7p2hl0SWKPRaDSVE0/bbHzzGC0H8BtkRtq2kDVLilG23sRUyIwBnxCR/6BOAvCCl+XVaDQajTN43UFgEYCrIAtixUFmA54D4Hlf5wA1L9TFAF4C8CZkvfCFAPoz83aPy6vRaDQaB4j6DAJuoHujaTQaTXQgomXMXGF5h1N21ueCggLk5ubiyJEjKCnRS5Jr7BEXF4eUlBTUqVMn8s4azSnIKelsCgoKkJWVheTkZDRv3hyxsbEIv9qwRhMaZkZ+fj527NiB+Ph4JCQkRLtIGk2lI+q90aJBbm4ukpOTkZKSgri4OO1oNLYgIiQmJiIlJQV79+6NdnE0mkrJKelsjhw5gtq1a0e7GJoqRq1atXD8uJ6WTKMJxinpbEpKShAbGxvtYmiqGNWrV0dxcXG0i6HRVEpOSWcDQKfONI6j7ymNJjSnrLPRaDQajXdoZ6PRaDQa19HOpgoxZcoU9OnTB6mpqahRowaaNWuGESNG4Icffoh20Rzh4MGDGDt2LJYvX17hs379+qFfv34n/p8zZw6ICDNnzox4XCLC2LFjHSypRqMJRDubKsLrr7+OSy+9FG3atMG///1vTJs2DU888QQAYNasWVEunTMcPHgQ48aNC+ps3nzzTbz55ptRKJVGozHCKTmosyry0ksvYcSIEfj3v/99Ytv555+P22+/HaWlpVEsmTd06NAh2kXQaDRh0JFNFSE3NxcNGjQI+llMTNnPPHbs2KC9pm666SY0b978xP/FxcV48skn0apVKyQkJCAlJQW9e/fGzz//XO577777Lrp27YoaNWogOTkZffv2xYIFC058npeXh9GjR6NFixaIi4tDixYt8Oyzz5ZzgCrlNXnyZNx0001ITk5G7dq1cf3112P//v0AgMzMTLRo0QIAcPvtt4OIQESYOHEigIppNMWhQ4dCHjMcq1atwrBhw5CcnIwaNWrg3HPPxfz58yN+T6PRBEdHNn7cdx+wcmV0y9ClC/Dqq+a/16NHD3z44Ydo2bIlhg8fjrZt29oqx/jx4/HKK6/g2WefRZcuXXD48GEsXboUubm5J/Z56KGH8PLLL+PWW2/FuHHjEBMTg0WLFiErKwvnnHMOiouLMWTIEKxduxZPPvkkOnXqhEWLFuGZZ55Bbm4uXn755XKa9913HwYOHIjPPvsMGzduxJgxY5CdnY3Zs2cjPT0dX331FS677DI89thjGDZsGACgVatWYc8j3DFDsXz5cpx33nk466yz8O677yIxMRFvvfUWBg4ciAULFqBbt242rqxGc2qinU0V4a233sIVV1yBRx55BI888gjq1auHQYMG4eabb8bgwYNNH2/hwoUYPHgwRo0adWLbJZdccuLvTZs24ZVXXsH999+Pf/zjHye2Dx069MTfn332GX7++WfMnTsXffr0AQAMGDAAADBu3DiMHj0aqampJ/Y/44wz8MEHHwAALrjgAtStWxc33HADfvrpJwwYMABnnXUWAKBly5bo2bOnofOIdMxgPPzww2jatClmzZqFuLg4AMCQIUPQsWNHPPPMM5gyZYohbY1GU4Z2Nn5YiSgqC23btsWKFSvwyy+/YPr06Vi0aBG+/vprTJo0Cc8888yJzgJG6d69O55//nk8/vjjuPDCC9GjR48ThhcAZs6cidLSUtxxxx0hj/HDDz+gWbNmJ6IcxeDBg/HEE09g0aJFJyIUALjqqqvKff/KK6/EjTfeiIULF4Z0DJEwe8z8/HzMnTsXY8aMQUxMTLlyDxw4EJ9++qmlcmg0pzra2VQhqlWrhj59+pyIIrKzs3HBBRdg3LhxuPvuu5GcnGz4WGPGjEFCQgI++eQTPPfcc6hZsyauuOIKvPjii0hJSTnR7tG4ceOQx9izZw+2bdsWcmqgwLaTtLS0cv/HxcUhOTkZO3fuNFzuQMweMzc3FyUlJXjmmWfwzDPPBN2ntLS0XDuYRqOJjHY2VZiGDRvitttuw6hRo7Bx40b06NHjxPT3hYWF5SKVQMMfGxuL0aNHY/To0cjJycF3332HBx54AHl5efj888+RkpICANi5cyfatWsXVL9evXpo0aIFvvjii6Cf+3dIAIDdu3eX+7+wsBAHDhxAo0aNTJ23nWOedtppiImJwd13340bb7wx6D7a0Wg05tHOpoqwfft2NGnSpML29evXA8CJnmrNmjUDAKxevRpdu3YFIONXFixYgFq1agU9doMGDXDbbbfh+++/x+rVqwFISikmJgbvvPNOhYZ+xQUXXIDJkyejZs2aaN++fcRz+OKLL3DLLbec+P/LL79EaWkpevXqBQCIj48HIKkuo0Q6ZiBJSUk477zzsGrVKnTt2lU7Fo3GIbSzqSJ07NgR/fv3x6WXXooWLVrg8OHD+P777/HWW2/hqquuQtOmTQEAF154IerUqYPbb78d48aNQ0FBAV544QXUrFmz3PGGDx+Ozp07o2vXrkhOTsaKFSvwww8/4M477wQgvcBU54AjR45g2LBhqFatGhYvXoz27dvj6quvxvXXX48PPvgAAwYMwIMPPojOnTujsLAQmzdvxtSpUzFlyhQkJiae0FyzZg1uvvlmXHPNNfjjjz/w+OOPo2/fvifaVtLS0lCvXj1MmjQJZ555JpKSktCiRQvUq1cv5HWJdMxg/OMf/0CfPn0wZMgQ3HrrrUhPT8e+ffuwfPlylJSU4O9//7vl30mjOWVh5ir36tatG4dj7dq1YT8/GfnXv/7Fl1xyCTdt2pTj4+M5MTGRu3TpwuPHj+eCgoJy+86fP58zMjK4Ro0a3KZNG/7444955MiR3KxZsxP7vPTSS3z22Wdz3bp1OSEhgdu2bctPPfUUFxYWVtDt1KkTx8XFcXJyMvft25cXLFhw4vP8/Hx+6qmnuF27dif2ycjI4KeeeoqLioqYmXn27NkMgCdPnswjR47kOnXqcM2aNfnaa6/lvXv3ltP7+uuv+fTTT+fq1aszAP7ggw+Ymblv377ct2/fE/uZOSYAfuqpp8ptW7t2LV999dVcv359jouL40aNGvEll1zC06ZNC/s7VMV7S6MxA4ClHMQuk3xWtcjIyOClS5eG/HzdunU4/fTTPSyRJhxz5sxB//79MWPGDAwcODDaxbGFvrc0pzpEtIyZMwK364S0RqPRaFxHOxuNRqPRuI7uIKCJOv369UNVTOdqNJoydGSj0Wg0GtfRzkaj0Wg0rqOdjUYTRaKdPSwoiJ42M3D0aHT1o3n+xcWAifHJjlNUBBw/7p2edjaaqFJQED2De+QIsH69vEeDf/4TSE4GvvkmOvoPPCD606dHR/+664DUVOCXX7zXZgaGDQMaNADWrPFev6gIOOccoGlTIDPTe/28PFnOpFUrIGBGJ9fQzkYTNfbsAX7/HcjKio7+jh1Ss96+3Xvt0lLg6aeBQ4eAEPN9usq+fcBrr0nN+sUXvdffuROYNEn0ozHb+rp1wHffAQcPAq+/7r3+jBnAkiXyO7z7rvf6P/4IrF0LZGcDH3/sjaZ2Nqc4JSXRiyz27pX3ffukHF5SVAQcOwbExEgtr6jIW/0VK+T827UDli0To+cl8+eLw+vVS/4uLPRW/4cf5D0jA5gzx/t78Pvvy+t7zezZQFwc0KOH/O01//sfUKsW0KaNd+evnc0pzJEjsjJpNML4wkKp1daqJYbm2DFv9VXqrGFDefdaXxmYxx+X9+XLvdVfsgSoXh244w5JZXqdSpo7F0hLA26+WSobXke3c+YAp58OjBgB/PGHRJhesmSJpLHOPVeeQb9lkzxh7lygf3+pbCxb5o2mdjZVhIkTJ4KIsGnTJsPf2b1bDP0HH7yP1q3bIC4uDqeddpqj5Vq5ciXGjh1bbjlpoKxh2DcZNfLyzB1XnW+mRU+Znw8QAb6VEkzr2+X338XRXXSR/O/VA69YuhTo2BHo3bvsfy9Zuxbo3Bno3j06+qtXi7HP8E2q4qWzLy2V37t7d6BbN7kX163zTv/4cWDTJrn+GRlATo6k09xGO5tTFGap3R89mo3nnrsDGRnnYNasWZg5c6ajOitXrsS4ceMqOBvVC6dWLSA+3vvIIj9fdKtXBxISvHc2a9YAZ5wB1KsHNGkC/Pabt/rK2LZqJb/B7797p11aKsa1QwegUydx+l7qHz0KbNsm+p07yzYv9bOypAydO8tv4LX+hg3yG5xxhrf6egaBU5S8PGknOXx4I0pKSjB8+Ej0VtVcDzh+XIx9TAxQo4a3XTABcTZqdYOEBG+7wCpjq1bUbt0a2LzZO/38fGDXLnE0RPLupX5Wltx/HTrItW/c2Ft9FUV06CCpvKQkYMsW7/RV8qF1a6BlS/nby/NXKdMOHaSy45W+jmyqKP369UPv3r0xc+ZMdO3aFYmJiejYsSOmTJkCQAzO2LE3YfjwfgCA664bACLCTTfddOIY7777Ljp37oyEhASkpKTg1ltvrRChFBcXY/z48ejQoQMSEhJQv359XHDBBVi/fj0mTpyIm2++GQDQpk0bENGJ1Fd+PhAbW4znn38eF1zQHmedFY+GDRviwQcfxPEAz7NlyxYMHToUiYmJqF+/PkaNGoUCG96htFSci2/RUsTHe9sFOztbjK1a4LRlS2+Nnco8tmhRpu+lsfvjD3lX5++1s9u4sUyfyPvzV1qtW0tFq1Gj6Jx/27ZAerqUwQt9Hdn4c9990loXTbp0cawv6ObNmzFq1Cg89thjSElJwcsvv4wrrrgC69evR0JCa9x++5MYPLgbRo26F48+OgHDhnVFamp9AMCjjz6Kl19+Gffeey9efPFF7Ny5E0888QRWr16NBQsWoFq1agCAa665BlOmTMF9992HgQMH4vjx45g3bx527dqFoUOH4oknnsDf/vY3fPnll2jcuDEAoEGDdKxZA4wZcwN++ulb3HPPaDRteg4KC9dh3LgnkZmZicmTJwOQZZwHDRqE/Px8TJgwAampqXj77bfx1VdfWb4uqjFWrYodHy8OqLgYiI21fFjD7Nol72plajXW4ehRIGANO1fYulXelbNp1Uq6AZeWSqTpNqozgG/RWLRqBXz7rfu64fQ3bPBOf9Mmuef8f3+vI8v0dCkD4J2z1c6mCrNv3z7MmzcPbdq0AQB07doV6enp+OKLL3DllWPQunUrJCXJ2ivNm3dA9+49Ub06kJmZiRdffBFPPfUU/vrXv544Xtu2bdG7d298++23GDFiBGbNmoXJkyfjtddew7333ntivxEjRpz4u1WrVgCALl26oHXr1gCkJ9ry5fPx7bef48MPP8Tw4Tdi40agXbuBSE2tixtuuAErV65Ely5d8OGHH2LLli1YuHAhevbsCUBWG+3UqZPl66K6OSvHoh66ggJvnU16urz7LhG2bAHOPNN9feVsVAqnVSv5TXbulPYjt9m+XZya6gnYqpWMufLK2W7fLoNZlVarVtIV2Ctnu3mzOHql1apVWVdwL8jKksGkCq+cnXY2/kRjdJmLtGnT5oSjAYDU1FSkpqYiKysLBQUSPvtz/Lg8gDNmzEBpaSmuv/56FPv1yTz77LNRu3ZtzJs3DyNGjMD06dNBRLj99ttNlauwEFi48AfExcXh8ssvB3Mxioulk8DgwYMBAPPmzUOXLl2wcOFCNGnS5ISjAYCYmBhcddVVGDt2rPmLgorORkU4BQXeGDvlbFRPPFXDzsryztnUqCHtFf7627d742yyssTRqOvvr+/FunOBxrZZM/nt9+2TGQ3cJju7/HVu1kx6hBUVeVPZycoq6xgByLWYO9d9Xd1mU4WpW7duhW3x8fE4fvw4CgvLjKxCNYPs2bMHANC6dWvExsaWex0+fBj79+8HAOzfvx9169ZFjUCvFYHCQiA3dw8KCwtRs2ZN1KoVi169YpGeHotU39OuNHbt2oU0ZRX9CLbNKKGcjVcDO3Ny5F2dgkqneNH9FJAIqnlzaa/w19+50xv9rKzyxlZFOF6df6Cz8Vo/J6esoqH0mcsqIW7CXLFS0aiRjDNyu0eojmxOQZglZRBYi1LOpp6vi8r06dORnJxc4fvq85SUFOTm5iI/P9+UwyksBOrUqYeEhATMnz8fgOTMa9cuSy019FmA9PR0rAky4nC3jQmdlFOp7rv7q1WTlIZXo+h37ZLxPer6N2gght8rY79tmzgbhTK2Xjqbbt3K/vfa2W3fXja+KFBfdQV2C+aKzsa/suHvBN1g/37pHOSv46/vlwhxHB3ZnIKUlsq7v7OpVq3MCA8aNAgxMTHIyspCRkZGhVcLX8vy4MGDwcx47733QmrF+xpE8v2mty0sBM499wIcP34chw4dQkZGBrp0yUD79mUaytn06tUL27dvx6JFi/zKX4ovvvjC8vkXFYmj8c/Px8V5G9kopwpIWdLSvDO2u3aVGRgAqFtX2q28qNmrmnW0IosjR4ADB6Knf/CgVOoCIxvAm99fdY4Idv5u6+vI5hQkmLOpXr2sZt+qVSuMHj0a99xzDzZs2IC+ffsiISEB27dvx4wZM3Dbbbehf//+6N+/Py6//HI88MAD2L59O84//3wUFRVh3rx5GDp0KPr164cOHToAACZMmICRI0ciNjYWiYln4pxz+uHaa6/FFVdcgQceeACNGvVAaWkM5szJxPfff4/x48ejbdu2GDlyJP7+97/jsssuw3PPPYfU1FS89dZbOHz4sOXzD5Ybj431NrLxdzaAGH8vjE1pqTTG+2chicTgeKG/d68YW39jl5QE1Knjjb6adNVfPz3du8hSpVCDRTbRcjZe6WtncwoSydkAwHPPPYfTTz8dEyZMwIQJE0BEaNKkCQYMGFCu08GkSZMwfvx4fPjhh3j11VdRp04ddO/eHbfddhsAoHPnzhg7dizeeecdvPvuuygtLcX06VvRrFlzfPLJJ3jjjTfw/vvvY/36ZxEXF4+WLZtjyJAhJ9pk4uLiMGPGDNxzzz3486z9AiAAACAASURBVJ//jKSkJFx33XUYOnQo7rrrLkvnH8zZxMV5t9TArl1A+/bltzVq5M0cdfv3y2Bef2On9L009r5e8OX0vYgsgunHxkrHAC+djX9lo149KUO0zt+zyI6Zq9yrW7duHI61a9eG/byqk5PDvGQJc1FR2bbMTOYVK7zRX7mSeevW8tt27JAylZS4r79qFfOWLeW3bd/OvHQpc2mpvWNHurdKS5ljY5lHjy6//U9/Yq5Xz562EX77jRlg/vLL8tuvvpq5TRv39b/7TvQXLSq/feBA5rPPdl9/4kTR37Sp/PazzmK+6CL39f/zH9EPvE2aNWO+4Qb39ceMYa5WreJzVrMm86hRzmgAWMpB7LJus4ky27bJvFheTpdSVCRpA9+4TABSsyouLot63II5dGQBuD/7bSj92Fj5zG393FzRD0yjNWwoUYfb0/YE9oTz19+50/1ZFIKlkQDvIivVryTw/L3SD3X+DRt6E9ns3i1RXOB4Ii/SqNrZRJH8fMlhFxaW3YReoIyt6voKlBl7t9stlDEPZuy90Ffr94Rydm53EggcY6PwqvtzOGOfl+f+VPuhnF2jRnJt3K7s5ORIG1HgeCqv0ng5OXKvBU6u7qWzC/ztlb7b56+dTRRRD3ZSkreLZ4WLLNw2toFjXCqLvlfOLljOHijLm7s91kLV7IPVrL3Sr1OnbF46RXq6VAT27XNXP5SxTU+Xip/b95/S96/oKX0vxtns3h36/N3W184mihw7Jl1O69aVm9yr3lCh0kiA+2WIZOyj7WzcTqOFimzUyHXfeFrXyMmR2a4Da/Ze6gczdl7p795dMaoCyrZFy9mlpQGHD3uTRg11/m5f+1PW2XC01kL2Iy9Ppg1JSir73wvCORuvjH31gH6Q6v9oORsn9I3cU6EiG2UAvDD2aWkVa9bRdjZenn84Z2djrLAj+mqpdDdgDh3ZpKZK5dfNWQROSWcTFxdXbpBhNCgpkU4BiYllc5R54WxCzW5crZoYILeNfag2GyIx+G5HFqGcjZpFwI5+fn4+YiNMbrVrl1QuatUqv12tGOpFzT7axj5YzdpLYx9Ovyo7uwMH5P6P1u9/SjqblJQU7NixA7m5uSgqKopKlKPC5Ro1xNDFx3uzgFg4Y696pLlJUZEYdf+ecIrq1b2JbELpx8Za02dm5OXlYefOnSfmdgvFrl3BH/a4OJmJOFo163r15B5wWz9czRpw19jJnHzRM7bFxRK5BEa1XumH6pwBeHP9PR3USURXALgWQAaAVABZAL4C8BwzH/HbLxnAiwBGAKgBYCGA+5nZkcVL69Spg/j4eOzduxf79+8vN7OxVxw7Jvnh2Fi5CfbvLxtd7SZqdluiivnpffvkYXTT6e3bJ2UItub6vn1lXZPdQvX+C6W/f7+13yA2NhZpaWmoXbt22P0Cp6rxx4u8eU5O+XnBFNWqSXTlpn5enrRLBDP2yclS2XDT2alzi1ZksWeP3N/RcraheiICZQ7IzfP3egaBhyAOZgyAHQDOAjAWQH8iOoeZS4mIAEwF0ALAXwAcAPAYgNlE1IWZdzhRkISEBDTxYj71EDz3HPD447KGR1ISMGEC8PHH0istMJ/uJFOnAsOHA0uWVJzO/eGHpfvlihXu6Q8YIMb8558rfjZuHLBsWdlKgm7wpz9JCtM3/2c5Hn1Uxj25uX7erl2hlxFITXXX2BQViUMNZmwAMThuGptQY1wAiTbr149ezb52bYkuo2XsvXB2oXoi+utXpTTaJcx8FTN/ysxzmflVAPcCOBtAP98+wwD0BvB/zPwZM//g2xYD4BGPy+saW7bIj646B7RpI7U+NxsIgYoLd/njtrEBQufMlb4XNftQkUVqqvvnHyqN5oW+urfC6bt5/cMZu2jrE7l//4VzNjVrSvtttNNobt5/njobZg5mSpf43tU8tMMAZDPzbL/vHQLwLYDh7pbQO7ZsKVspEShbrVGtougWu3bJgxWsaUE9bG4OrAvVZgBImdzu/hnO2KeliUF26/xVGilaabRwxg5w39hH0ne7smPk/KOt7/b1DzagFJBxT7VrV63IJhh9fe8qi34GgNVB9lsDoCkRebCWovsELmCk/lYT5blF4Foq/qSlSYopN9cd7XANtEofcO+Gj2TsU1PdPX8jxkZNZ+OmfrjI0gtjG0rfK2dXWfW9SGMG6/aucPv8o+psiKgRgKcBzGTmpb7NdSHtNIEoE1BxNS851h1EtJSIlu51OxflAIHpHK+cTaQGasC9G17dyOEedv/9nCbUGJfKpu/W7WskjXXkiEyj5KZ+qA57XqTRgs1eoPDC2dapU3E5doUXzi7Ubw+4f/5Rcza+COUbAMUAbvb/CECwvshhm82Z+R1mzmDmjPr16ztXUBc4elRe/kYnOVlytl5ENpGcjdvGNlJk49YNH2r0fqC+W+cfrr3MC30jNWu39UNF1Uo/L0+eDbf0wxlbZezdGglhVN8tQnU790o/Ks6GiBIgPc5aAhgS0MMsFxLdBKIimmBRz0lFsBoukUQ3anEjt4jUZgG4Z+yN1KyB6EcW0XJ2Xpx/7drha9Zu64dydJVFv7BQUq1u6UeKLNxsM410/lUusiGiWACTAfQAcFGQsTNrIO02gXQAkMXMLtV7vCOU0WnSxN3IprRUbqZoGVsjOXM39StDZFGtmnTxDYYX1z9SzRZw7/xDzUtWWfS9iKxD3XtAWZvhAReq0yUlFVdoDaa/f797A7s9dTZEFAPgUwADAAxn5kVBdpsKoBER9fX7Xm0Al/g+O+kJVcN229mEWktFUbeuGMNoOZukJHm5mcZSgxeDUbeujPdwUz8treJaIgq3nV2kNIrbxtZIzT6a+l5EVtHS379fKpuRrj+z7OsGppwNEQ2yqTcBwJUAXgJwjIh6+r3UQqVTITMGfEJE1xDREN82AvCCTf1KQbjIJifHvZmXI6VxYmLc7f4ZqYEWcDeUV2mEUMZeDSyMlrGrVUumLYpWGklFXG7oh5sEUuGmsc3PDz17gRf6qp3WiLN1Qz9SeyngfmRtNrL5kYg2EdHDRGSlFf5C3/vjEIfi/7oNAJi5FMDFAGYAeBPA1wBKAPRnZpebz71h1y5pJK0b0DLVpIk8lG4tYhSpzQJw39iHu9kBdxspI6UxAHfHukTSV+OfonX93Ywsjx6Vxv9oObtwsxco3IysIrVXAu4aezPn79b9b9bZnA8ZhPkMgO1E9B//dFckmLk5M1OI11i//XKZ+RZmrsvMicw8gJlXmSxrpSVUDVut1ujWin2R2iwAd41tpJy50q/Kzi6SvlvX//hxWazPiH60jF1CgkS+0dJ3c+ZtI5GFm8bejLOrFM6Gmecw87WQ0f5PQibUnE1E64holG8CTU0EQtVw3V4aOFIaDYi+sY9mZAG4F1moBloj+tEyNm7qGzG20dZX2QY3fn8j+m62GUZqL/X/rLKk0QAAzLyfmV9k5rYABgHYB+AfAHYS0UQi6uRkIasa4dYBB9yNbGrWrLhKoz/K2bgx1sBIZJOaKoMaS0qc1VbGPlrOTk2DEy1nZ8TYuKlvJLKoDPpu/f5GnI3qvOLW+SckVFxHyZ86daSD0t13O68P2OyNRkQXQSbS7AlgD4CPINPPLCeiP9kvXtUkVA27bl1pIHbL2YSbPUCRliazMjs91sBMGqe01PkpY8wYe9W+4CRGokqgzNg57eyNRhZuGVujkZXb+hGWG3I1soqJCd0TUuGms2vQIPyM8kRA48Yyf5obmHY2RNSAiB4noq0AvgNwGoAbADRh5rsAtAbwNoC/OlrSKoJaQCnUzLONGrkb2RhxNoDztSszNVs39I20VwHu5c2NdM4A3BtYaCaN5sZkpDk5cn9HMrZuGvvkZKnMRdJ3K7JMSwu+aF+gvlvONtKz5zZmuz5PBrANMtX/9wA6MXNfZv6cmYsBgJlLAPwHQJRPrXKiaq2hjE7Dhu622RipWQLuOZto6RuNLNx2dkbPXzknp1DHi1SzV5Fl4MJ6dtm9W3qbVY+wglZammg7PRlppG7X/vpuOBsjz57Sd/q3B4yfv5uYjWzaALgPQCNmvpuZ14TY73cA/W2VrIoSqYbtZmRjNI0GOP/AGW0zcFs/2ucfTf169ULPS6ZQBsmNyoaRmrXax+nJSM3oHzrk/DIXRjrHKH03I6toYtbZXAzgvWBTxhBRdSJqCgDMfISZ5zpRwKpGpNy5cjZO5+yPHZMZfSMZOzeNDRA9Z2M0slCfO1273LVL1hEJN6AVcNfZmDH2bkRWRvTduv5m9d1Ioxp1NseOycspSkokWjzZnM1WyFLOwejs+1wTBiORTX6+LA/thm6kG75ePckrRyuNc9pp0kDphrMxYuzdTKMZMTZuOXsj7XVu6htN40Rb3w1nq+YkNKPv5PmrNriTLY0Wbpr/WAAurvFYNYiUTmrYUN6dbrcxmsapVk0MrtPOZvduYw20anleN5ydEWMbF+fOWAuj+m46ezPOxkl9NVVNtCKbvDyJ6s1Edk7+/rm50jHIjLN18vyNZhXcJkJzHUBEp6H8lP+NiKhlwG41AIwE4ELTVtVi166yLs7B8B9rc0awua9t6ALGDY4bxs5ozapBg+jV7AF3nN2uXUDPnpH3c2N+Nmbj51+rlkR/Tp6/WpAtWpGF0c4p/vpOXn+j3c7d0j9pnA2AUQCegixoxgD+G2I/8u2nCUMko+vWwE6jaTTAHWNrputlWhqwY0fk/cyQk2PM2Ct9p419NJ3twYMydsqIPpHz+maMXWKiOLxo6VdFZ2NG302MOJspADIhzuR9AH8DsDlgnwIAa5n5N0dLVwWJVMNUaTQ3nE316pKmiUSDBsDq1c7q5+QAGRnG9k1LA5Ytc07bTM0ekPN3Uv/IEUnlRCuyMhPVAs5HtkZ7IlYGfTU/mxv6Rq6/moz0lEyj+SbAXAUARMQApjGzw73wTx1ycoBzzgn9eUKCpNncaLNp0CD09Pr+qJptaamx/Y1gJrJp0KBsxUIn9A8fNp7GAdwz9mb0161zTt+MsVP6mwOrkzYwk8ZS+tFKoyn9aEUWsbHOT1mTkyMRY7hpqrzA7EScH2pHYx2jNWw3xtqYrdkXFTm3YqBqoDXzsJeUOLeIkxVjq6IRJ7AaWTjV/d2KfjTbDNzSj9QT0i39XbtkKW6jxt4NZ5uWFn6qGi8w0kFgFoA/M/N639/hYGYe4EzRqh5qsFgko+uWs2nWzNi+/j1ijKTdImHF2Cj9UEsom0FFiSpFaVR/926gRQv7+srYG9VPS5Mpaw4dku7aTumbaTNSo/gjDQI1gtF5wfz1Z860r+uvb2RAqyItDfjNwQaB7Gx5po0ae6cjq8owVQ1gLLLxv0Qxvv9DvTxdZvpkw2gN2w1nY7TrK+B890uzDZRON5KadTZVTV/VrGvXNq7P7NwofjVVTaR5wfz1Dx50bhS/WWPrdGSRnW38twecj6zMdE5xEyNtNv39/u7nammqOEbTGY0ayc1WXBx5LikjFBWFnvwzGE47G7ORTVUz9tnZkjM3auz9r3+7dvb1VUXDaM3aP7IzYyRDYdbY+4/ib9rUe33/KWsiDQI2Qna28c4xSt/pNNq55zp3PKvoSMRDjNbwGzYs6y7rBGp9GtWtOhJuRTZW0mhOkJ0t+fJwa3m4rd+wobk0CuBsZGM0qgWcP/+dO83X7CuDvhNT1qhl3s3op6VJe+HRCpOCmaegQCqaTlQa7GJ21ufhRHSz3//NiGghER0hov8SUZT7O1RuzEQ2gHOpNLM1+9q1nR3Yl50tOXujzqZ2bRn06mRkYeZhc3rKGivGxkl9o1PlBOo7aeyNVnQAZ52NMvZm9J28/ocPi+OI1u+vbI6Z83cLs5HNEwD8m2z/AaAxgHcA9AEw1pliVU1ycsSI1qkTfj+nnY06jtEbXg3sc9LYNGhgPCXotL5ZYx8bKw3K0dJ3esoas5GNk8auqEiOEy1jv2+fdLawou/E9VcVvWhFlurZPxmdTSsAvwEAEdUAcBGAB5j5QQBjAFzqbPGqFuqhj5ROcSuyMVu7jFbNFnC2R45ZY++kvurubkY/Jsa5Rbzy86X9wYyxS0qSlKMTv7/qwm3m91eRZbSMrZOTgZrNKgDOOlsrz75bmHU2CQDyfX+fA+lgMN33/wYAlSAzWHkx2iskJUVq104N7MzOlpqymW7E0XY2TulbyZkDzjXSHjki08VHS9/qVCVO6VsxtvHxMrA5Ws7GjcgmWm1WJ3Nkkwmgt+/v4QCWMfMh3/+pAA4F+5JGMJrOiImR/ZyMbNLTzY3Gj7azSU8vyzfbQXWhNWvsGzVyxtlbMTZO6psd0KlwqvutVWPnVGRpNoUMiLOrV8/Z39/M9U9NlcqhE8//zp1yPsnJ9o9lF7PO5m0AY4loKYA/A/i332e9AKx1qmBVETMNtU6OtTHbGwcoP7DPDvn5MhOBWWPTqJH0oikosKdvx9jv3ClT5kRL34nJSK0YW8C5yoZVZ+OkPpF5Z+vU9c/OlpSk0Z6QgDgapyqbqqIX7dkDAPPT1bwG4CYACwHcwszv+n1cC8BEx0pWxSgokHUtjN70Tjobs71xAOe6f1o1No0by7vd2qVVY9+4sTjafTYnZ7Kj74Sz3b5d3s2OV3HS2Kv5vszqOxHZ7twpkYLZmRAaN3bm+TPbXqdw6vm3klVwC9PjbJj5U2b+CzN/FLD9zsBtmjLMTgbYsKGzbTZWaraAfYNj1dmo/e3WLu1EFk7qW6lZ+3/fKllZ0uBvdtqbhg0lBWl3eWIVVZudUFWlEe3OD2eloqX0nYosrTibxo2d0z9pnY2CiFKJqGngy8nCVSXMTgbZqJE0Lh85Yk9XpbHM3vCqnHZrV3YjG7v66oG1auzt6m/fLl3dzc6469T5b98uUY3ZNEqTJvJu1+BZNXaNG0tbm93JWO3o79kj3abtsH27dWdn97dnPomdDRHVJqIPiCgPwC4AW4O8NEEwOxmiU8bOas1eGRuVhrFKtCObrCzphZeYaO57Thn7bduMT4Dqj5Pnr35LMzj5+1up2Tupb9XZAPYiy+Ji0bfy+zduLBXNw4et6x84IA67sjgbszNvTQBwOaRjwO+QRdM0BjA7ktff2bRvb13Xaj/7tDTJczvxsJttIAVkFoGaNZ0x9lbm10pLk4ZaJ4y9FX0nI5suXcx/zwljr2rWF15oT/+ss6zpq8jIamQByO/fvLk1fdXBxE5lY+dO43PqBdP3P1a0MetshgB4mJknuFGYqoyassXomhpO5eyVsTB7w8XEiMGLVs2SyJm89bZtQIcO5r/nVI+gbduA3r0j7xdI7drS1mLn/AsKpK3QSmSjfjM7v/+BA9LmY8XZOuHssrLkPVrOfts2ebca2Sj900+3pm/n/N3AbJsNQQZvakySnV1WWzaCU8tD27nhmzSx72wyM61pA/YbaZmtp7Gc0D98WBrZregrZ2vn91dlt2Js4uOlYmTn/Lf6kupW1gRSPcjs3H+ZmfJuJTJxIo2pjH200qh2rr8bmHU2kwBc4kZBqjpme4QlJUnDshPOpl49OZ5ZmjQpe2CskplpPQ1h19ju2ycdJKw6G7v6dhw9YL+RWBlqK5GN+l60jH1MjH1nr/StGNs6deSZceL3t+LsnWizzcyUCXWNZlPcxmwabTqAV4moFoDvAeQG7sDMkVbzPCXJzjZ/0znRI2XbNuvGvmlT0S8pMR6R+XP0qBh8q/qq+6tVfTs1S6U/fXrk/UJhx9go/blzreur87fjbDZutK5vx9grfTvObutWmfzVSgcFIvvObts2MfQ1apj/bkKCVBLtOtvmzSvHgE7AvLP5xvfeAjK4U8GQFBsDsGAWqj7Z2UDPnua+48SUJdu2Wc/5NmkiPWqsLqJl19g0biyOxqq+XWPv3yPISiOt3cimcWN7ztZuZNO4MTDLRtVx61aJEKwubd2kCbBwoXX9zEz57a1cO8B+m6XVzilO6W/dar2i5wZm02j9Q7zO93vXBFBYaG0Bo4YN7UU2zPbaTOw20qqcsdUbXpVbGW2z2DX2qtxbLXbo37YNiIuzvv57ixZl3WetkJlpvWYNyO9/+LD17rd2UqhAWQcRq1MGZWbaa69o0cL6bw9IZGn13lP6qsJmBbvn7zRmp6uZG+nlVkFPZtSATiuj2HftkpqtFey2WShnY7Xdxm5k07KlvG/ZYu3727ZJ3r1uXXv6dpxN06bmR88H6ls9/02bgNatrX0XsD+w066xa9KkbElzK9it2bdsKVF1Xp7579rtnKL0t2yxNovC4cMyPdbJHNkAAIgohYguJqKRRFTXty2BiPQy00GwM2VKSYn1+clUzd5Omw1gL7JJTDS3tIE/qtx2nE2zZtZz1k44O7vGxo6+XWejym7F2TLbN/Z29PPzxVHYcXZ2Kht2K3pKX52HWex0znALszMIEBG9CGAHgKkA3gfQ3PfxNwAed7R0VQS783NZbbexm0Y67TQZWGk1lLfbQJmQINfAqrHdvLnMYFghOVmugR1j36qVdf0mTaS9wYp+fr5EJHacjfru5s3mv7tvn0QEdoy9HX27FS2grOxWrr/qWGHn97ejf9I7GwCPAbgHwNMAzoZ0ClB8C+Bih8pVpbAz8y8QvTYLInngrTzsgDMNlCqVYJbSUjH2bdpER//AARm9bke/enX57azoq9q4HWdTv77M/LBpk/nvOmHsWrSQe9CKvt32QsBeZKmcjZ3f345+ZRtjA5h3NrcBeJqZnwOwPOCzTZBlozUBZGeL4TA7zbqqFVk19pmZYiys9gYCxFhZ6f7qRBoFsG7ss7Oldh8tZ+OEsQHEWFjRVwbaTs1aVTasGHsnavYJCRLdWdF34vxTUiSyt/r7V6tmz9jb6aCycaM8+2ZtjpuYdTaNACwK8VkhAAtDB6s+O3aYXykTECeRkmLtYQPkhmvTxl4/+zZt5GYvLjb3vX37gEOHnDH2O3fKPFdmcMrYt2wp52+2R5ST+nacjZ3IRn3fSmVnwwa5353Qt3L/r18v3a6t9gQE5LlRv79ZNm4UZ2F2HR1/7KSRN2wA2rWrPGNsAPPOZieAjiE+6ww963NQ7HSBtPqwAWXOxg6tW4ujMdsjbf16eW/Xzp5+y5ZlPXvM4KSxLyw03262aVOZsbKrv3ev+aUmNm2SNierPfEUrVtbq2xs2CC1+vh4+/pW7n+njK3VyNKJZ8+Ovjr/yoRZZ/MlgL8S0bl+25iI2gJ4EDKdjSYAqzP/AnLDWnnYCgvFSLRta01XoWqmZsuwwTeDnhPOBjBfu9y4UQyd1QGNgfpmH/iNG+U3T0hwRt/s+dvtiaZo3Vq6H5vtkbh+vTPGrnXrsijZDE4ZWyvdj5nl93fi+luJbI8dk9/rZHc2YwGsBzAPgMrkfwlZbmAjgL87VrIqQklJ2QJWVmjdWr5vNo2kUj9ORDaA+XabDRvE2Nvp+gmUGVuzzu6PPyRfb3WMi0Ll/P/4w9z3Nmyw7+iBsutvVv+PP5ypWVupbJSWir5TzgYwl8o7dkxS107p5+ebG1i7a5dM1eTU779jh7kVU9WzelI7G2bOB9APwEgACwDMBLAEwB0ABjGzzXXtqh45OZKCsJNGYzZfu1HGye4Nn54uY2WsRDZt2lifKkSRliZtV2vXmvve6tVAx1AJXxM0ayYj8NetM/6d0lJgzRpn9Nu3l1SQGf3DhyXt6IS+crZmfv8dO8RAO+lszOire98JfbU8hZnrv3q1vHfq5Jz+BhNz7TuVVXAas+NsEgD0giyaNgXAOAA3M/OHzGwyq3tqYHdNCatpLPXA2a3dWu2R5FQahUgeODPO5tgxcc5OGNtq1cTgm9HfulWMrRP6iYnS0GxGX+3rhH56uszCoNrgzOjbWfRP0aqV3ANW9K3OCeiPMvZr1hj/zu+/y7sT19+K/tq1EtE7EVk5iSFnQ0TxRPQaZJbnuZC2mc8h6bT9RPQSEcW5V8yTF7uTQSpnYdbYr1snPdnsNhCrMpip2eXnS9rDyqJlwTjjjOgZW0DOw8zDrmq2TuqbOX8n9WNi5PorA2qEVavkvXNn+/qJieJwzOivXCkpXCcqO6mp8hyZvf4NGjjT7bh1a+nRZkZ/5UpxNFbnxHOLiM6GiAjAd5DBnD8AuBPAhQAu8v09A8D9kEhHE4DdyCY5WaYaN1OzA+ThdCKMB+Q4W7YYzxuvXi2pJCvLEQejQwfpkWV0jiwna5ZKf/t24xNSKn2nnG2HDvL7G+0Rtnq1RCN228sUnTrJORltJF+5UrTtjO8Kpm+UVavEQVY3O6d9CKw4e6fuvdhYcRxm9Fetcu7ZcxIjkc0VkBmdr2Dmy5j5PWaezsw/+v4eAeBKAIOJ6DJXS3sSkpUlD53VdcQBuXFVbdUIqs3AKWdz5pliaIzW7p2s2QJlRtvoA7d6tdTq7HY7DtQ3Gt2tXi1dVmvWdE6/sNB4I/nvv8t37HaOUHTqJD3CjM7RtWqVc7+90t+4USLmaOiryNaIsy0pca69LlDfCAcPSjbFyfN3CiO347UAvmDmr0PtwMxfQXqlXe9UwaoKdro9Kzp1EgNmtGaZmSlRiJPOBgB++83Y/itXyuhlp6bKUA+ucmKRWL5cvmO3c4JCnb9R/aVLna1ZmtEvLQWWLQO6dnVe30h0kZ8vDdROGrszz5TzMlLZyMmRiWud1O/cWYy4kbFm69bJNXDy+nfuLBUNI5G1ekZPVmdzFoBpBvb7DoCDl7hqYHfmX0AM55Ejxgc2Op1GatFC0jJGnc2qVWIgnKpZp6dLDnzp0sj7lpTIfmef7Yw2IOefnGxMf/9+MQxO6nfsKOviLFkSed+NG2VMSo8ezumrSosRZ/fbb86mUM3qL/dNouWkfkaGvBu5/osXy7uT11/pL1sWeV+1z8maRqsPwMj48SwAEVe7JqLGRPQGES0kojwiYiJqHmS/ZCJ6j4j2EdExIppJRA7V1b1BdVm2Oz+YetiMtwgV2gAAIABJREFU5q2VUzjjDHu6ipgYKYORh724GFixAjjrLGe0AemN1L27MWO/bp1EdU4+7ETywBsxNmofJ/Xj4qSmauT8f/1V3p10dikpcg+rY4djwQJ5N7sqbThat5ZUtDLkkfSrVSsz0E7QqZO0nRi5/osXS1mdGOOk6NZN3o3oL1woldv0dOf0ncKIs0mEdHWORCEAI+OlWwO4CsABAPOD7eDrlDAVwAUA/gLgcgCxAGYTUWMDGpWC3btlcJfdG09FKEYji8WLpdtprVr2dP1Rxj5SI/Vvv4mxP/fc8PtZ0V+/PvK0LW7ULJX+6tWR2w0WLxbnpAyEk/rLlkWeo+3XX6WtyIlux/707GlsieYFC8QxWVnGOxQxMeI8jej/8otUdJIcnKUxPl4idaORTffuzkX1QJmzj6TPLOfv9LPnFEYvSSMiahnuBcCoE5jHzGnMfBGknScYwwD0BvB/zPwZM//g2xYD4BGDOlFHdVe262xq15ZjGKnZMMsN72TNFgB69ZL1SSI5PFWzPeccZ/UzMuTcIl2DhQtlAkYna5ZKv7hY2qPC8csv0qBrp0NIKP0jRyL3SvzlF3G0TrVXKXr1klH04VbtdNPY9ewpzj5cZaOoSO59p+89oKyyFW7V3MOH5flw+tlT+pEiy6wsmcPPjfN3AqPO5r+Q6WjCvUI5jnIws5H5c4cByGbm2X7fOwRZM2e4wTJbx8o6rEFwajJIQB62RYsiF23bNmkgdfqGVzdwpNrlL7/ITLV25yQLpGdPiRjmzQu/308/Af36OVuzBIDeveV9zpzQ+xw/DsyfD5x/vrPaANCnj7zPnh16nz17JNU5YIDz+r16ybuqTARjyxaZqsUNY9erl0R14VJpK1ZIhcgN/T59xJksD1xYxY+5c8UZuXH9+/QRZxJuJpH5vjzRyexsbgZwi4GX2s8JzgAQrLPvGgBNicihTqUVWfr3mfgjsQtyNx+wfayNG8sWwLLL2WdLT5tIPWLcyNkD0qMuPT28sWEWZ9C7t/NTmycnSw+fmTND77N5s4zeHzjQWW1AFhI788zw+gsXSppt0CDn9Vu2lFTKTz+F3mfWLHl34/y7dJG0bDj9//1P3gcPdl6/Vy95lsJd///9T+47N5y9Oma4858xQwahKsfsJMqBRbr+9etXzp5ogAFn45uKxvDLoXLVhbTpBJLre08O/ICI7iCipUS0dK/R0X9BoNT6aHv8N+we+y/Lx1Bs2iQ9mZwYXKYaXCOF0vPmSb7aqW7PCiKJGGbODN1usGqVhPEXXOCstmLgQInujh4N/rkyRG4Ye6X/yy+h221mzJD0Vb9+zmsTif7s2aFTOTNmSArR6fYiQBrIBw0Cvv8+dHT9/fcSxTsx23EgtWtLJUY5tGBMmyaVrPr1nddPS5O200jOpk8f+8sqBKN9e6nshdIvKQF++AG48ELno3qnqKTFAgEIdkuHrC8z8zvMnMHMGfVt3G0dru2MH+gCNP7qNeOjyELg1JoWgNSqk5LCp3EAueH79bO3aFMohg6VVE2odpNpvg7ybjqb4uLQD9yUKVL7d2tOqEGDgIKC4PrMwOTJYmyc7JgRqH/wIPDzzxU/KyoCvvlGjI3T7TWKiy6SNptgA4yPHRNHOHSoO9pKf9Wq4DMw5+RIA7qb+oMHS2Xu4MGKn61eLe1pbukTif6PP8o9GMivvwK5uXKNKiuV1dnkQqKbQFREYz/HFYIaNYBvTx+NWnl7gIkTLR/HyTUtAHEe558vN1sotm0TTbdq9hdcILWmaSFGXU2dKg3ZDRq4o9+nj8z1NinIqkl794qjveYa91YnPP980f/ss4qfLV8uk59ee6072oAYssTE4PozZsgYn+uuc0//oovk2n7+ecXPvvxS2qwuc3EOkYsvlvcvvqj42aefyrub+ldfLTM5TJ5c8bPPPhMnf+WV7uofPBg8uvvoI7FdblX0HIGZo/YCcBskgmkesP19ADuC7D8RwLZIx+3WrRvb4aEHS/lX6sElbdoyl5ZaOkZ2NjPA/MYbtopSjn/+U465cWPwz996Sz5fs8Y5zUD69mVu06biZVmzRrRfesk9bWbmu+5irlGD+ciR8tvfeEP0f/vNXf0772ROTKyoP2oUc2ws8/797upfdx1z3brMx4+X337VVczJycwFBe7qX3QRc3o6c1FR+e29ezO3tf64GKZnT+bTTy+vU1rK3L49c69e7mqXljK3bi3PgD9FRcxNmzIPHuyuflERc/36zJddVn770aPMtWoxjxzprr5RACzlYPY+2EavXmGczQjf9r5+22oD2A/gjUjHtetspkxhvgEfyeWZOdPSMWbOtPX1oGzaJMd85ZXgn/fvz9yunbsP/McfSxl++qn89vvvF2O7Z4972szMCxaI/uuvl20rLhYHmJHhrjYz88KFFX+D/fuZa9Zkvv569/VnzBD9d98t27ZpE3NMDPNDD7mv/803oj9pUtm2RYtk2wsvuK///vuiNX162bbvv5dt//63+/rPPy9ay5aVbVPPxDffuK8/ejQzUfkK5auviv7PP7uvb4RK5Wwgk3teAeBfPqfyJ9//fX2fx0AWZ9sO4BoAQwDMgaTXmkQ6vl1ns2cPczzy+VhivYrVCIO88opc3ZwcW0WpQLduzF26VNy+Y4fchOPGOasXSH4+c716zBdeWLZt927mpCTma691V5tZHGnv3lK7PnBAtqmIb/Jk9/WZmfv1Y05JkfNmZr7tNjH2v//uvnZpKXOPHnL++/eLox0yRK5/drb7+sXFzJ06MbdowXzwIHNhoZSnfn3mw4fd18/PZ27WjLljR+a8POZjxySqadHC/aiOWc45OZn5vPPkWuzbx9yoEfNZZzGXlLivv3evVGyGDpV7YccOiXTPP999baNUNmfDIV5z/Pap60un5QLIA/ATgM5Gjm/X2TBLhPDflg8zV6smv6hJbr1VHkCnUemiFSvKb3/6adn+xx/Oawby4oui9fnn8sBdfjlz9erM69e7r83MvGSJ/Cx9+0pZ4uIkheF2CkexejVzfDxzhw7M11wj12L0aG+0maVWHRvLfMYZzAMGiP6ECd7pz58v179zZ+Zzzy27F7xi2jSpWHXvzty1q/ztH+m4zcSJcs59+4qji49nXrzYO30VyQwZwtyypVQ0vHr2jFCpnI3bLyecza23Mnet48tbjR1r+vs9erhT29i/v2IUcfCg1LSHDnVeLxgFBZIfr1ZNctiAGH0v+fhjabsBxOC53VYSyI8/Sp4+Pp75wQfF6XrJ99/LtU9NDZ1WdZOvv5ZoIj29fErPKz79VCKc5s29dXSKN96Qc2/fXu4FLyktZX7uOanMduxYedJnCu1sTPLBB3J1jpw7hLlJE1PWpKREHMK999ouRlDGjJGyTZsmN95NN0ntbulSd/SCceAA81/+IjXrDz7wTtef3Fzmdeu8SV8Eo7Q0etoaTWUllLNxaC27qoea32nBGbdj8C9XSJ9jg53YN2+WcQdqHRCneeIJ4LvvgBEjZEzJmjWyzY3BfKE47TTg9de90wtGcrK8ogWRe92sNZqqRmUdZxN1WreWEbsfHbhEFiJ/5x3D31XzN3Xv7k7ZatSQqUluvllmhJ0wAXj6aXe0NBqNxgl0ZBMCIqB/f2DmT3HgkTeB/vGyzDJoYKGIxYtl8J1Ta9AHo1494O233Tu+RqPROImObMLQv7+sSbO5/20y+dAHHxj63uLFktJyYk40jUajqQpoZxOG/v3lfUZmG5lw7L33Iq5eVVgoU507vXiXRqPRnMxoZxOGli1lXZbZswHccYfMX6/mcQ/B4sUyUV5lXS1Po9FoooF2NmFQ7TZz5gClwy+VWRgjdBSYNUu+17evN2XUaDSakwHtbCIwZIjMKLz4twTgxhtlHvsw6+XMni1roNcNNme1RqPRnKJoZxOBiy6Shv5vvgFw++2ycMiHwdeIO3pUVrJUbT0ajUajEbSzicBpp0lKbMoUSF/mc86RjgJccW23H36QDgKXXOJ9OTUajaYyo52NAUaMkFX41q+HdBTYsAGYP7/Cfl99JYMse/f2vowajUZTmdHOxgCXXSYrVH74IWQpvjp1KnQUOHZMppAZNsy9ZXk1Go3mZEU7GwM0bChtNx98ABTFJgLXXw/897+y6LePzz4DjhyRKWQ0Go1GUx7tbAxy550ym8CkSZCOAgUFwCefAJDmmwkTgI4d9fgajUajCYZ2Nga56CKgc2dg3Dig6IwuQEYG8O67ADO+/BJYuRJ44AE9C7BGo9EEQzsbg8TEAM8+K8sHjB0L4K67gNWrceSbWXj4YaBTJxmGo9FoNJqKaGdjgqFDgVtuAZ5/Hhi/43oU1muAtTeNR3a29BfQHQM0Go0mOHpeYpP885/AoUPAo2MTkIv7MB6P4scXlqFnTw9XLtNoNJqTDB3ZmKRGDemItngxcM6Hd6G0Vm2cv/SFaBdLo9FoKjXa2Vike3dg+I11EHP3n4EvvwR+/z3aRdJoNJpKi3Y2dnn4YRnk+dhj0S6JRqPRVFq0s7FL3briaKZNA+bOde64QeZe02g0mpMV7Wyc4C9/ARo3Bh56SJaPtsOsWbIqaHw8kJgIXHwx8PPPjhRTo9FoooV2Nk5QowYwfjywdCnw5pvWjsEMPPccMGAAkJkJ3Huv9LNetgw47zzg/vtlSmmNRqM5CdHOximuvVZWWhszBsjKMv/9v/4VePxxmXdt/XrgpZekn/XmzRI5vfqqDPQ5dsz5sgejpEQme9NoNBoH0M7GKYiAt96SCOX662WRNaO8/Tbwt78Bt94KfPwxkJBQ9lliIvD668D770uKbcgQGejjFlu3AtdcAyQlAbVryyyko0cDBw64pxlISQmwZIksELR5s3e6/hQXA3l50dHWaKoizFzlXt26deOo8emnzADz/fcb2//bb5ljYpgvvJC5sDD8vl98wVy9OvPZZzMfPGi/rIHMns1csyZzrVrMd9/NPH4886WXMlerxly/PvOUKc5rBvL118xNm8o1VK/u3Zl/+sl9bWbm9euZL7uMOSFBtJs0YX78ceY9e7zR37KFeeRIud41ajBnZDBPmMCcl+eN/oYNot+wIXNyMnO/fswffcRcXOyN/po1zNdey9ygAXNKCvOQIcxffcVcUuKN/pIlzJdfzpyWJq9LLmGeNo25tNQb/XnzmC++WH7/9HQpy7x53mgzyzN++LCtQwBYykHsctQdgxuvqDobZuZ77pFL++qr4ff79VfmxETmbt2YjxwxduwpU8Th9OzJfOiQ/bIqfv5ZjFuHDszbtpX/bOVK5q5d5ZyefNK9B/+FF0SjSxfm//xHyvTqq8zNm8v2O+9kzs93R5tZrm2NGsx16jD/5S/Mzz0nDz4Rc926zJ995p42M/PUqcxJSXJP3HijVFi6dZNzb9WKee5cd/W//FLOPymJ+frrmf/0J+b27UW/c2fm3393V/+jj+TerlOH+f/+j/n228t++/POY960yV39V16RilW9esw33cR8yy3MjRqJ/qBBzDt2uKddWsr81FOi1aAB8623yj2QmirbLruMed8+9/RLSphHjRKtp5+2dSjtbLyksFAiAvXDBTPO8+czn3Yac8uWzDk55o7/9dfyUPbq5YzDyc6WG7x1a+bdu4Pvk58vDx8gD0GkKMwsEybIsa++mrmgoPxneXnMDz8sn3ftyrxzp7PazMwzZsg17dFDroc/q1fLtQaY773X+XNnZp4+XQxdt27MWVll20tLpWwtW0oE/PLL7tSyf/6ZOTZWztP//EtKmD//XIxefDzz++87r80sxwWYBwxg3ru3bHtREfN774kDSkpyL7p+6SXRv/TS8lmDwkLm118X7Xr1mH/80R39N98U/ZEjmY8dK9uelyeVnthYiXQWLnReu7SU+a67RH/UKNv3t3Y2XlNQwHzDDXKJ+/ZlnjlTbqJt25j/+lfmuDjmdu2Yt261dvyvvhLjeM459sLewkKpNdaowfzbb+H3LS0V5wkwDx1a/qGww4IFci4XXyzGJRRTp0qar2lTcQBOsWYNc+3azB07hk5PFhVJpAEwn38+c26uc/qrV4t+p06hKw9HjkjtFpD7yskILytLUkatW4c+r927mQcOFP1HH3U2uv31V3keBg5kPn48+D7bt0s6lUiMr5MO99tv5byuuir0/bd+vfw+1aoxv/22c9rMzHPmyP0/dGjodOWKFRLdxsczT5rkrP7rr8v5jx7tyHXVziYalJYyv/OO5J792yBUDd5uWPzf/8rNb8fhPPiglOeTT4x/56235KHv1Yt5/35ruorduyVV0bIl84EDkfdfvlyisDp1mGfNsqfNLNetXTsxtoHpw2BMnCi1zLZtmf/4w77+oUPMbdrIOflHNMEoKSlz9j17mo+Ig5H3/+3de5BU5ZnH8e8PZwBZFNBgKkFx1pAyJaJGjasBVpmNBokMmrDspuKuypqbrKWlrokb/8iuICYmWF7iajaWlPeAqCRVogneFyUGFFZBjECIGMAQxeDAwAwzz/7xnN4Zxx7s2+nTzDyfqq6eOTPd73P6cp73dt6z01tTgwd70t2b1lbvygQfS6hEZePtt80OPdTs8MM/+vuwc6eP5+Ra1z0lpmK88YZ/lo4//qPHxbZvN5s0ycu/8srKJNwNG/z4cOSRHz0Ou3Wr2bhxXv4111Qm4S5Z4olu8uSKVSAi2WRpxw5vicya5bWI116r3HPnEs7YscUnnNxkhhkzii93/vzO1tm6dcU/3sxrkY2NPhj/8suFP27DBh9bqq83u/PO0so28y/rtGnePfX004U/7tlnvUtl2DCfVFFO+VOn+vtXzCDwgw96S/Sww3w8rZzyv/Y1/wwsXFj4Y+bM8crGCSeU16XZ1mY2YYK//8uXF15+LuGOH//BLrdi7djhrZWDDiq8h6GtzeyiizoTbnNz6eU3N/v45JAh3nIqxK5dnT0m555bXsLdssUnghRa0StQJJvebP58P2CNGVP4IOqyZf4lHz/+w2MkhXrmGT/gDh9eWl9yrluqlISxbVtnt86VV5Y2W+rmm/3xs2cX/9i1a33wvK7O7I47in+8mU9+AJ8YUazly71FWM44xvXXd9aSi5WbzDBihLc2S3HZZV7+3LnFP/b++71L6VOfKvxA3VVHhx+sJbPHHiv+sXPmeCXluOMKaxHne45p07z8Rx8t/rEzZ/pr9/nP9zzOujetrd69X2xFrwCRbHq7xx/3GtrQod5i2VsTe9kyTxIjR5b2Qe1qzRqvGQ0Y4IOchTbtcwPCF19cetmtrT5jCrwboJguvcWLvWV01lmldx9s2+azlMAnMBST8BYv9kQ1ZUrp3SGbNnWOY1x3XXHPs3ChHyynTi29/BUrvHU1aFDxM/Xuuqv89//5572iM3SoT08uRm6copyZV48+6mNthxxS/PTk3MyzH/yg9PLnzfNk0dDgU7aLMWOGFd19XqBINn3B+vU+myo3iP3kkx88kO7e7bX5QYM80ZTa/dXd1q1mEyd6uU1N3s21N3fe6QfI008vf2ZXR4fvU329dwkUMltoyRI/l+joo8vvPmht7ZzJc8ophdWyly6tXPk7d/r4H3jiKqRb66mn/CD1uc+V1w1kZrZ5s9euc+MohezPI4/4+3XaaeW//+vXmx1zTGfiKmR/7r7bP39NTeWPU6xe7RMr+vXziROFdGvlWtTnn1/+uMuLL/qYV12dt9D3NsEmJ5foLr+8vLJ7EMmmr9izx+yWW7zGBz7w/YUvmJ1xhrdmwH/vPr23XO3t3rWw//5+ILv00g8feP/4R7MLL+yMoVKz2cy8Wyl3TkhTk7fe8sV4223eChs1ymc4VUJHh9cQhw3z577iivwngXZ0eJfRwIHeGqzUeRsdHT51d+BAr2nPnp1//K693WdS1df7a1WpE1Xb2vwA1q+fD3bfemv+wfY9e7zrsNInJre0+OcNvFtv7tz8XcOtrZ0H2gkTKnei7Pbtfl4M+Odq3rz8rdyWls6uw8mTKzeF/p13vIUKZqNH++y6fEm0ublzgsf06amdL9dTspH/rXc58cQTbdmyZVmHka2WFliwABYtgvXrfQmY0aN9KZozzvDlddLw5pu+xtsDD/iSLw0NcNhhvs7aypWw336+qOjMmdC/f2XLbmmBG2/0BU3ffx+OPRZOPRU+/nHYssWXv3njDWhshHnz4OCDK1v+5s3w3e/CPff4fk6cCOPG+WUotmzx92PFCt/20EMwfHhly1+71hdwXbTIlxpqaoKTT/afN270fV650hd7nT8fhg2rbPkvvwyXXALPPef7fM45cNJJMHiwx/bzn8Pq1TB5si/LNGRIZctfssT3/6WX/D3/8pfh+ON9yafXX4f77vM4zjvPF8wdNKiy5T/+uK/8/uqrMGIETJ3qn8H+/WHVKv9cbNwIF13kn9O6usqVbQaPPAKXX+5LTo0aBVOmwJgx/l1/5RV/zd9+25efuvZa6JfOamWSlpvZiXlizL4lUulbn27Z1IrNm70WO22ad5eceabXKtM+C9zMu3Juvrnz/CHwqb0TJvgJimkvPfLaa966yZ39nrt99rM+mSDtpV9+8xuzCy7wGXPdy587N92lXzo6fGbftGk+y6r7skMPPJDu69/ebrZokdnZZ/t73rX8sWN9YkOa5e/ZY7ZggU+Rzn32wCfwNDb6WF2adu/28bPGRp8tmiu/rs6/g0uWpFu+Rcsm9FVmfmmG+vrUanJ7tW0bNDd7Lf7AA6tbtpm3tnbu9FZcpVsyH6W9HTZt8pXKR4yAAw6obvltbV7+zp3euh48uPrlv/mm3zc0fHCB3WrYtcv3f88eGDmyauX31LKJZBNCCKFieko2cYmBEEIIqYtkE0IIIXWRbEIIIaQukk0IIYTURbIJIYSQukg2IYQQUhfJJoQQQuoi2YQQQkhdJJsQQgipi2QTQgghdTWbbCQdJulBSX+RtF3SQ5JGZh1XCCGE4tVkspE0CHgS+AxwHvBPwKeBpyT9VZaxhRBCKF4FL6hQUV8HjgCONLO1AJL+F3gD+CYwJ8PYQgghFKkmWzZAE7A0l2gAzOz3wBJgSmZRhRBCKEmtJpvRwKt5tq8CjqpyLCGEEMpUq91oBwHb8mx/F8h7BShJ3wC+kfzaLOn1Msr/GPDnMh6/r4v9j/2P/e+7yt3/w/NtrNVkA5Dvqm7q8Z/Nfgr8tBIFS1qW7+I/fUXsf+x/7H/sf6Wft1a70bbhrZvuhpG/xRNCCKGG1WqyWYWP23R3FLC6yrGEEEIoU60mm18AJ0s6IrdBUgMwNvlb2irSHbcPi/3v22L/+7ZU9l9m+YZGspWcuLkSaAGuxsdvrgEOAI4xs+YMwwshhFCkmmzZmNkOoBH4HXA3cC/we6AxEk0IIex7arJlE0IIoXepyZZNFvrywp+SpkpaIOkPklokvS5ptqQDso4tK5Iek2SSZmYdS7VImiTpWUnNyXdgmaTGrOOqBkljJf1K0p+SfX9J0vSs46o0SYdKulnSC5J2Jp/xhjz/N0zSzyT9WdIOSYsljSmn7Eg2xMKfwBVAO/DvwETgv4BvA7+W1Oc+I5K+ChybdRzVJOmbwEJgOXAO8PfAfGBQlnFVg6RjgMVAPb4u41eA3wJ3SPp2lrGlYBQwDT+F5Ll8/yBJ+ESsicDF+OtRjx8PDy25ZDPr8zfgEvxgO6rLtr8G9gCXZR1fFfZ/eJ5t/4xPzGjMOr4qvxZDgS3AV5P9n5l1TFXY5wZ8Ms6lWceS0f5fC7QCg7ttXwq8kHV8Fd7Xfl1+vjD5jDd0+58pyfYJXbYNwVdwuanUsvtcrbUHfXrhTzPbmmfzb5P7EdWMpQb8EFhlZvdnHUgVTQc6gNuyDiQj/YE2POF29R69rPfHzDoK+LcmYJOZPdXlcX8BfkkZx8Ne9UKWIRb+/LBTk/vXMo2iiiSNw1t0F2UdS5WNA9YA/yhpnaQ9ktZKmpF1YFUyN7m/SdInJQ2V9HXg74AbsgsrM3s7Ho6UNLiUJ63ltdGqqeiFP3szSSOA/wQWm9myrOOpBkn1wO3Aj8ysnEVc90WfTG7X4+N26/Axm1sk1ZnZjVkGlzYze1XSacDDdFY02oBvmdkDmQWWnYOADXm2v5vcDwOKPgUlkk2nohb+7K2SWstCfLzqgozDqabvAPsDs7IOJAP98BOmzzezh5JtTyazlK6SdJMlHfe9kaRPAwvwmvu38O60KcBtknaZ2b1ZxpcBkcLxMJKNi4U/AUkD8VkoRwCnmtlbGYdUFckU9+/hA6YDJA3o8ucBkoYC75tZeyYBpu8dfPblr7tt/xU+I+kTwKZqB1VF1+ItmbPMrC3Z9oSkg4EbJd1f4FhHb/EuPR8PocRjYozZuD6/8GfSjbQAOAmYZGavZBxSNR0BDATuwb9IuRv4tPBtQFnnGNS4VT1sz9Vke/uBdgywskuiyXkROBg4pPohZWpvx8M3rcRVXCLZuKwX/sxUci7NvfiA6BQzW5pxSNW2ApiQ5waegCYAa/M/tFd4OLn/YrftXwTeMrMtVY6n2rYAx0nq32373wC76Byr6Ct+AYyQlJskhKQDgcmUcTyMbjT338C/AgsldV34cyM+aNzb/QQfEJ4F7JB0cpe/vdXbu9PM7D3g6e7b/dw2/mBmH/pbL/Mo8BRwu6SPAeuBqcAZ9I1xu1vwE1h/KelWfMymCT/X6gYza80yuEqTNDX58YTk/kxJW4GtZvYMnlBeAO6R9G94y/4qvKX7w5LL7cXjfkVJ+u1vAE7HX9Qn8JPcNmQZVzVI2kAPl3IF/sPMvl+9aGqHJANmmdnVWceStqTmOhtPMsPwqdDXmdl9mQZWJZLOxCeJjMa7VNfhS+3f3tvG6pLPdT7PmNlpyf8cBPwIOBt/PV7AT3BfWXK5kWxCCCGkLcZsQgghpC6STQghhNRFsgkhhJC6SDYhhBBSF8kmhBBC6iLZhBBCSF0kmxDKlFwkOtBuAAAB2ElEQVRa96NuGyQ1JD+fn3XMIVRbrCAQQvlO6fb7w8BK4Ptdtu0GNif/u646YYVQO+KkzhAqLFmR4X/M7NysYwmhVkQ3WghVkq8bTdJcSW9JOlHS85JaJL0u6UvJ3y9LuuC2S1ooaXi356yTdJWkNZJ2S9ok6cfJ5SJCqBmRbELI3oHAXcDPgHOAPwELJP0YX3F6BnBp8vNPuj32HuBq4D7gS/j6Zv+Cr+IdQs2IMZsQsncAfgniZwEkbcLHfM4CjsotBCnpaOBiSfuZWbuk8cA/AOeZ2V3Jcy2W9C6+Yu9xZrai6nsTQh7RsgkheztyiSaxJrlf3G3F4TV4BfETye8TgVa8FVSXu+FX2AT42zSDDqEY0bIJIXvvdf3FzFqTa+l0v/xu7roqufGYQ4D+QE9XTjy4UgGGUK5INiHsu97BryQ5voe/b6piLCHsVSSbEPZdj+EX/BpiZk9kHUwIexPJJoR9lJk9Lel+4EFJc4AXgQ6gAZgEfMfMfpdhiCH8v0g2IezbzgUuBqYD38NXKtgAPA68nV1YIXxQrCAQQgghdTH1OYQQQuoi2YQQQkhdJJsQQgipi2QTQgghdZFsQgghpC6STQghhNRFsgkhhJC6SDYhhBBS93/EUcecUazvVQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Import libraries\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.integrate import solve_ivp\n", "plt.rcParams.update({'font.size':16})\n", "\n", "# Parameter values\n", "alpha=1\n", "b=1\n", "gamma=1\n", "q0=0.5\n", "delta=0.5\n", "beta=0.2\n", "a=10\n", "\n", "def popdyn(t,x):\n", " \"\"\"Function to run population dynamics\"\"\"\n", " q=q0*(1+delta*np.sin(2*np.pi*t))\n", " sdot=(a-q*(x[0]+x[1]))*x[0]-b*x[0]-beta*x[0]*x[1]+gamma*x[1]\n", " idot=beta*x[0]*x[1]-(b+alpha+gamma)*x[1]\n", " return [sdot,idot]\n", "\n", "# Run the model\n", "ts=np.linspace(0,10,1001)\n", "sol=solve_ivp(popdyn,[ts[0],ts[-1]],[10,10],t_eval=ts,method='Radau')\n", "\n", "# Plotting commands\n", "plt.plot(sol.t,sol.y[0],color='b',label='Susceptible')\n", "plt.plot(sol.t,sol.y[1],color='r',label='Infected')\n", "plt.legend()\n", "plt.xlabel('Time')\n", "plt.ylabel('Density')\n", "plt.title('Resident dynamics')\n", "plt.tight_layout()\n", "plt.ylim([0,30])\n", "#plt.savefig('resident.png')" ] }, { "cell_type": "markdown", "id": "f0ead1ca", "metadata": {}, "source": [ "The next code cell focusses on the evolution of host avoidance. In particular it produces a pairwise invasion plot to show the evolutionary outcome when there is a trade-off between avoidance and reproduction in the host for a given parameter set." ] }, { "cell_type": "code", "execution_count": 5, "id": "e46e8a1d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Resident strain 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 , 10 , 11 , 12 , 13 , 14 , 15 , 16 , 17 , 18 , 19 , 20 , 21 , 22 , 23 , 24 , 25 , 26 , 27 , 28 , 29 , 30 , 31 , 32 , 33 , 34 , 35 , 36 , 37 , 38 , 39 , 40 , 41 " ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAELCAYAAAAP/iu7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de7ylY/3/8dfbKDkkW04ljCHjVIqpSIeh/DoPQjmEkSil0HFEmhBTSScpozIRpYQZ1JcYI/mmTCqMcpoZjL4yGKfBOH1+f1z3Msuatfdaa6/73vc6vJ+Px3rsve51H651773Xe1+H+7oVEZiZmRVpubILYGZmvc9hY2ZmhXPYmJlZ4Rw2ZmZWOIeNmZkVzmFjZmaFc9hYriTNlzS/14/Zy7r1fEoaLykkTS67LLYsh00fkjQ6+6OsfizJPmR+JmmjssvYbSRNy87j68ouSy+pOq+Vx3OSHpb0v5IOkqScjzcxO87EPPdrsHzZBbBS3QL8Kvt+VWA8cACwq6Q3RsRtw9jnO3IqW6cfs5d14vk8DbgXGAWMBnYDpgKvAz5VXrGsWQ6b/vbviJhceZL9l3gGsD9wFDCx1R1GxB15Fa6Tj9nLOvR8/jgi/lF5ImkKcB1wiKSTImJeeUWzZrgZzZ4Xae6iU7On4wAkbSLpW5L+IWmRpCclzZF0tKQX1e6jXnt/VVPIRpK+JOlWSU9JmizpiOy1d9Vsc162/PSa5Xtkyw9scMwBSSdI+rekx7Oy3yTpVEmr1Ky7qqTjs3WflPSApAslbdX6WXzBfivNldMkbZzt82FJj0q6qLq5UtJK2fKbBtnXS7Jtb6la1urPZhNJZ2Xna4mk/0q6VtIXa9ar22eTvZ9pkv4v+/ndKen7ktaos25ImiVpbUlnSro/+znMkrR1i6dyGRExB5gFCNim0fqSdpF0dXaOF0u6rvp3KFtnGumfLYAzqpru5rdbXnPNxpZV2wb+QVLT2kzgcmAF4O3AcaRA2qWFfZ9C+mC4BJgOzAVuzF7bAbgUnq9hva1qebXx2ddZg76BtP2lWfkuA2YALwY2ItXWTgAey9ZdA/gjsBlwZVa2l5OaaXaS9M6I+HML77Ge0cC1wD+B04HXAO8HtpC0RUQ8ERGPS7oQ+IikrSLinzX7eD+pqfPbVcua/tlIWhf4K/Ai0rmfD6wObAF8DPjmUG9A0ljgT6RzcyGpCXZr4NPA+yRtGxELazZbDbgGeAA4E9ggK/MVkjaLiHuHOmYTmuqvkfQ54CRgIfBz4ClgV+An2bn+TLbqhVmZdyado0pN6qE2y2kAEeFHnz1IH34BXFizXMC07LUzsmXrAi+us97p2XpvqXltPjC/Zllln/OAV9a8thywCLi2atmW2fozs6+vqnptDnD3UMcEXpttd3Kd975q9fsBfpmtu1fNehsDDwM3NnlOK+/xdXXOcwCfq1n/jNrjAu/Oln2zzv4vyF7bqGpZ0z8b4DPZsgl19v3yJn6GV2bb71uz/Jhs+c9qllfe9/cAVS3/arb8yOGe12z5FsDjwHPAhtmy8dm6k6vW2wh4GrgHWKdq+Sqk8A/gbVXLJ2bLJpbxt9nLDzej9bdNs6asyZJOBmaT+msWkf77JyLuiYinqjeK9Ff5o+zpO1s43kkR8Z+afT0HXA1sI+ml2eLx2dfJ2dcdACStBWzOELWaGk/ULoiIRyrvJ6vVfAi4JCJ+WbPe7aQP7S0lbdnk8QYzF/hOzbJp2ddxVcv+APwX2DOrnZGVcwB4LymQn+9PGebPpt45eWCowktan/Qz+XtEnFXz8jezMu8l6cU1ry0GvpyVqWJa9nUcrflE9nt6rKSzSP01KwI/jKH7a/YhteB8I6pqUhHxGCn4IP3OW8HcjNbfxrL0D+5p4D/Az4DjK3/AkpYDDiT9x7cFqWZQ3XzxihaON3uQ5bOADwBvAX5P+mC7JSL+KGlB9vwsmmhCy9wM3AQcqTQU+RJSoN1U88H3BlLNahXVvzZjs+zrptn+huuGLFSr3ZN9Xa2yICKelXQuqRbyNuCq7KXdSc2Av6jeQYs/m4uAE4ELJf2aFGx/ioi7mih/ZTj3rNoXIuJJSdeSmp7GsrRZFOC2iFhcs8ky77tJH68cEngUuJ70u3rGoFskg5adVFurXscK5LDpb9MjolGfyw+ATwJ3AueThp8+RfqwOIzUT9Cs+wZZPiv7Ol7S/5A+aH+bLbuKpSEzvmb9uiLiGUk7AseS+gjem710t6TjI2Jq9nz17Ovbs8dgVh7qeE14uM6yZ7Kvo2qWn00Km71ZGjb7ZOufW7Nu0z+biJgn6c3A14APk400lHQdqYnv6iHKv2r29b+DvH5vzXoVy7zv7GcDy77vRl4fVaPRWjBo2SPiYUlLWLbcVgCHjQ1K0trAIaS27e0i4omq195E+kBrxWB36vsH6YNpB9J/6GuyNFBmAftkTTk7AAuiiaG5kTqrD5H0KVIf0E7A4cBpkhZGxAXAI9nqX4+Io1t8L4WIiL9Kug3YQ9KngbVJ4fu7iLi/st5wfjaRBh3sImkF4I2k2uShwO8kbR4Rdw9SrMp5WnuQ19euWa+TVJf9BYEj6WWkQO7Ecvcc99nYUDYkNctcXv1hltk+r4NU9dtsTfoAhBeGDaT/xjel+f6a5/cdETdExLeBvbLFE7Kv15ECcNvhlLtA5wADpAEDe5F+Br+oWWfYP5uIWBIRV0fEF0l9c6sAOw6xSaVG8bbaF7LgehPwJGmEWqcZtOwsrc1W15iezb62WvOyBhw2NpRKe/52NR3WmwBH5nysWaQ/8COAmyPiv/B8R/0C4AtV6w1J0oaSNq3zUuU/8Ceyfd8LnAe8Q9IhdfaznKShmteKUgmWfbLHo6ShuNVa+tlIekO962GoOSf1ZP06V5EGcXy45uXPk/qGflU7WKFDnEMKkC9IWrOyUNLKpCZFSMOyKx7Mvq47MsXrH25Gs0FFxH8kXUC6JuE6SVcCryTVDC4lXYuSl0r/xJrAr+u8tk/2/awm9rUVcEHWcT2H1Fe0Iem6k8eBH1etewipxnSqpI+RrkV5DFgf2A5YC3hJi++lLRFxu6S/ks77i4Aza2svw/jZ7ENqVrwSuJ00Uuz1pObFW4CLGxTrENJ1NudI2gO4lVQTfRdpSPuXhvl2C5Wdyy8D3wBulPQbll5nsyFwakRcVbXJtaRa2mGSVgXuBx6OiB9hbXHNxhrZn3StxJqkC/heR5rK5gtDbTQMf2dph/Ksmtcqz5vqryGNevsG6fd7AvBZUlPPr4FxEXFDZcVs2O92pPe0HLAf6YN1a9IFiXu3/lZycTYpaGDZJrSKVn42vyQNO14P2Jc0n9irSM1o20fE40MVJiL+RRq9dzZp1ODnScPQfwhsGxGDDf4oXUR8kxS+t5Mugj2ENLz/4Ij4VM26D5CabOeTztGJdGiQdhu9cCRoCQWQ1iNdg7ATWRs0cHgzQzIlDVb44Y5cMTOzApQaNpJWIo2mWQIcTeqsPR5YCXhtnTH6tdsH6b+102peuqHRf2pmZjZyyu6zOQgYA4zNOoKRdANwG+kirpOb2Mc9EXFtcUU0M7N2ld1nM4E0BcftlQXZlevXkK5INjOzHlB22GxB/WlA5pA6H5txSDZd+uOSZkp6a37FMzOzPJTdjLY6aVRIrQdJF7U18gvSkM3/kKYv/wIwU9JOETGr3gaSDgYOBhg1atQ2L33pS+utZi1Yd11fkmDWS+65J01h99BDD90fEWs2WL0pZQ8QeAr4dkQcWbP868CXIqKlMMxmDb6JNAX9WxqtPzAwEDvsUHu7FGvVlClTyi6CmeVk0qRJz39/wQUX/C0iWp2hu66yazaLWDoZYrUB6td4hhQRj0q6hDQTro0AB41Zb6gOmSKUHTZzSP02tTYnTRM/HGLwCR/NzKxK0SFTUXbYzABOkjQmIuZCus85aSLBls9ANr3E+4C/5FhGM7OeM1IhU1H2aLTTSdNCTJe0s6QJpAkH76bqQk1JG0h6RtIxVcs+L+l0SXtLGi9pf9KQ6XVIF4iamVkdIx00UHLNJiIWZze5+g7pTowCriBNV/NY1aoizQhcHY63kCbT2xV4GemeFNcAB0bEX0eg+GZmXaWMkKkouxmtMn35kLMHR8R8Xni7WyLiItKtbs3MbAhlhkxF6WFjZmbF6ISQqSi7z8bMzArQSUEDrtlYG3yNjVnn6bSQqXDYmJn1gE4NmQqHjZlZF+v0kKlwn42ZWZfqlqAB12zMzLpON4VMhcPGzKxLdGPIVDhszMw6XDeHTIX7bMzMOlgvBA24ZmNm1pF6JWQqHDY2LL6g06wYvRYyFQ4bM7MO0KshU+E+GzOzkvV60IBrNmZmpemHkKlw2JiZjbB+CpkKh42Z2Qjpx5CpcNiYmRWsn0OmwgMEzMwK5KBJXLMxMyuAQ+aFHDbWMl/QaTY4h0x9Dhszsxw4ZIbmPhszszY5aBpzzcbMbJgcMs1z2JiZtcgh0zqHjZlZkxwyw+c+GzOzJjho2uOajZnZEBwy+XDYmJnV4ZDJl8PGzKyKQ6YY7rOxlnj2AOtlDpriuGZjZn3PIVM8h42Z9S2HzMhx2JhZ33HIjLzS+2wkrSfpPEkPS3pE0vmS1h/Gfo6UFJL+VEQ5zaw3OGjKUWrNRtJKwExgCbA/EMDxwJWSXhsRi5vczxjgKOC+ospqZt3NIVOuspvRDgLGAGMj4nYASTcAtwEfB05ucj8/As4GxlL+ezKzDuKQ6QxlfzBPAK6tBA1ARMyTdA2wM02EjaS9ga2BvYDziyqomXUXh0xnKbvPZgvgpjrL5wCbN9pY0gDwHeCLEfFgzmUzsy7loOk8ZddsVgcW1Vn+IDDQxPbfAm4FpjV7QEkHAwcDrLjiis1uZviCTut8DpnOVXbYQBoUUEuNNpL0VmA/YOuIqLeP+geLmApMBRgYGGh6OzPrXA6Zzld22Cwi1W5qDVC/xlPtNOCnwAJJq2XLlgdGZc+fiIgluZXUzDqOQ6Z7lB02c0j9NrU2B25usO1m2eMTdV5bBBwBfLet0plZx3LQdJeyw2YGcJKkMRExF0DSaGB7oNFv0g51ln0XGAV8Gri9zutm1uUcMt2p7LA5HTgUmC7paFL/zXHA3aRmMgAkbQDcARwbEccCRMSs2p1JeghYvt5rZtbdHDLdrdSwiYjFknYkDV8+izQw4Arg8Ih4rGpVkWosZQ/VNrMR5pDpDWXXbIiIu4DdGqwznyZGqEXE+HxKZWadwEHTO3IJG0lrA49ExBN57M/M+ptDpve0HTaSzgQ+QBpyfBcwG5gdEae0u2/rHL6g00aCQ6Z35VGzeTOwDvA0aSjyNtnDzKwpDpnel0fY/At4NiKeI103Mwc4M4f9mlkfcND0hzzC5mjgNElfiYj/5LA/M+sDDpn+kkfYTCY1pe0iaRFL+2xOymHfZtZjHDL9KY+weS3wqoh4Orv6/w3AuBz2a2Y9xCHT3/IImzlk18Bk18PMB36Tw37NrAc4ZAzyuSJ/MfBrSa/OYV9m1kMcNFaRR83mVtJQ5z9KEkv7bCbnsG8z60IOGavVdthExFcr30taF19nY9a3HDI2mKbCJguRlYHbhrorZkTcA9xDunWA9QjPHmCNOGSskYZ9NpIOJHX6/wu4V9JHs+UnS7pH0t8kTZK0UrFFNbNO5KCxZjRTs/kCcC7wI2B3YKqknYAPAb8k3YPmS8BeknaIiAeLKqyZdQ6HjLWimbBZD/hERFwDXCPpcdJdNKdGxCHw/KzPlwNHksLJzHqUQ8aGo5mwuQ94VdXzM0ihckllQUT8V9IU4Cs4bMx6kkPG2tHMdTa/ByZL2jR7fifwa2BuzXoLgA1yLJuZdQgHjbWrmZrNMaQpaG6SNBv4M3AR8EzNev8PeDjf4plZmRwylpeGYRMR90vaHtgD2IU0SOAwICQ9AlwP3EsaMOBbC5j1AIeM5a2p62wi4ing7OyBpLVIk21WHuOBUcABknYF/gZcFxFfLqDMZlYQh4wVZVgzCETEfcDvsgcAkl7JCwPoo4DDxqxLOGisSHnMjQZAduO0GXj2ALOu4pCxkZBb2Fhv8lQ1vcshYyPJYWPWZxwyVoY87mdjZl3CQWNlcc3GrA84ZKxsLYWNpLnArhHxzzqvbQnMiIgxeRXOzNrjkLFO0WrNZjSwwiCvvQRPV2PWERwy1mmG02cz2M3TxgEPtVEWM8uBg8Y6UcOajaQjgCOypwFcJOmpmtVWBFYHfpVv8cysWQ4Z62TNNKPNBa7Ivt8fmA0srFlnCXAz8JP8imZmzXDIWDdoZiLO6cB0AEkAx0bEvILLZWYNOGSsm7Q0QCAiDiiqIGbWPAeNdZuWr7ORtD+wF7A+aQRatYiIjfIomJktyyFj3arV62y+AnwNuAn4B6mvpi2S1gO+A+wECLgcODwi7mqw3QbA94HXAWsBi7NyfSMift9uucw6iUPGul2rNZsDge9FxBEN12yCpJWAmaTQ2p802u144EpJr42IxUNsvgpwP3A06ZbUqwIHAb+TtFtEnJ9HGc3K5JCxXtFq2LycdEvovBwEjAHGRsTtAJJuAG4DPg6cPNiGETGHFH7Pk3QJMA84AHDYtMkzPpfLQWO9pNWwuQrYilQbycME4NpK0ABExDxJ1wA7M0TY1BMRz0h6GHg6p/KZjTiHjPWiVsPmcOB8SQ+Q7tL5YO0KEfFcC/vbgmxYdY05wB7N7EDScqSZENYg1ZQ2AQ5roQxmHcEhY72s1bC5Nft6xiCvR4v7XB1YVGf5g8BAk/v4JvC57PvHgD0j4orBVpZ0MHAwwIorrth8Sc0K4pCxftBq2BzL4HOjDVe9/amF7b9LmiZnHWA/4BxJu0fExXUPFjEVmAowMDCQ93sxa5pDxvpJqxd1Ts75+ItItZtaA9Sv8SwjIhaQRqMBXCxpFnASUDdszDqBg8b6Tdk3T5tD6reptTlprrXhmE3qWzLrOA4Z61fDChtJWwFjWXYGASLizBZ2NQM4SdKYiJib7Xs0sD3Q8l9lNljgLcAdrW5rViSHjPW7VmcQWA24BNi2sij7Wt330UrYnA4cCkyXdHS2n+OAu4HTqo67ASlAjo2IY7Nlk0lNcNcA95L6bA4E3gjs3cr7MiuKQ8YsafXmaSeQLux8GylodgV2BM4m3Yrgja3sLJshYEfSKLezsv3MA3aMiMeqVhUwqqa81wNbAj8ALiONSnsSeGtE+L46VjoHjdlSrTajvYs0N9q12fMFEfE3YJakH5Gub9mvlR1mc6Dt1mCd+dSMUIuIGaRmOLOO4pAxW1arYfMKYG5EPCvpSeClVa+dj+/UaX3MIWM2uFbD5l5gtez7O4HtgFnZ841zKpNZV3HImDXWatj8iRQwF5P6WL6ajR57hjRrs5u1rK84aMya02rYfA14Zfb9t0iDBT4MrEQKmk/nVzQrk2d8HppDxqw1rc4gcAfZNSwR8TRpTrLPDbmRWQ9xyJgNT9kzCJh1BYeMWXtaus5G0txs9oB6r20paW4+xTLrHA4as/a1WrMZDawwyGsvATZoqzRmHcQhY5af4TSjDTYt/zjgoTbKYtYRHDJm+WsYNpKOAI7IngZwkaSnalZbkTRPmS/qtK7lkDErTjM1m7lA5c6X+5Om8F9Ys84S0i0BfpJf0cxGjoPGrFgNwyYipgPTASQBHFe5HYBZt3PImI2MVq+zOaCogpiNJIeM2chq9X42MxusEhHxjjbKY1Yoh4xZOVodjbYcy45Geznprp0LSfelMetIDhqz8rTajDa+3nJJGwEXkm6uZtZRHDJm5ctlupqIuEPSFNLknK/PY59m7XLImHWOPOdGWwhskuP+zIbFIWPWeVqaG20wklYHPks2I7RZWRw0Zp2p1dFo81h2gMCLgbWz73fLo1BWrm68l41DxqyztdqMdhXLhs2TpFtE/ya7343ZiHHImHWHVkejTSyoHGYtcciYdZdmJuLcsZUdRkSjCz/N2uKgMes+zdRsLmdp05kGWSey1wIYlUO5zJbhkDHrXs02oz0K/DZ7LC6uOGbLcsiYdb9mwmYHYD/SSLM9gAuAn7u5zIrmkDHrHQ2vs4mIqyLiQGAd4BPAWsClku6SdKKkzYoupPUfB41Zb2l6NFpEPAmcA5wj6RXA3qQazxcl/SgiDi2ojNZHHDJmvWm409U8AMzPHlsAAzmVx/qUQ8ast7U6g8D2wL6kvpsVSHfwfB/wh/yLZv3AIWPWH5q5zmZjUsB8BBgN/BH4PGnGgMcKLZ31LIeMWX9ppmZzK/AIcD7wMdLUNABrSVqrduWImJtf8awXOWjM+k+zzWirAhOB/ZtY1xd1Wl0OGbP+1UzYHFB4KaynOWTMrGHYRMTPiyyApPWA7wA7kaa8uRw4PCLuarDdOOBg4G3A+sD9wNXA0RExr8gyW3McMmZWkcvN04ZL0krATGBTUhPdvsCrgSslrdxg8z1Jw66/D7wHmARsDczOAsxK5KAxs2p53hZ6OA4CxgBjI+J2AEk3ALcBHwdOHmLbb0TEwuoFkq4B5mX7PaaQEtuQHDJmVk/ZYTMBuLYSNAARMS8LjZ0ZImxqgyZbdqekhcC6RRTWBueQMbOhlB02W5AuDK01h3ThaEuyedrWAv7VZrmsSQ4ZM2tGqX02wOrAojrLH6TFKXAkLQ/8GFgI/HSI9Q6WNFvS7CVLlrRyiL7RbIA4aMysWWWHDSy9MVu1wW7SNpRTgDcDH4mIegGWDhYxNSLGRcS4FVZYYRiH6Q9DBcmkSZMcNGbWkrKb0RaRaje1Bqhf46lL0omkYdD7R8RlOZWt702aNIkpU6a84LmZ2XCUHTZzSP02tTYHbm5mB5KOIg17/kxEnJVj2QwHjJnlo+xmtBnAtpLGVBZIGg1sn702JEmfAY4HjoqIHxRURjMza1PZYXM66Z440yXtLGkCaXTa3cBplZUkbSDpGUnHVC3bE/gu8D/ATEnbVj02H9F3YWZmQyq1GS0iFkvakTRdzVmkgQFXkKarqb59gUgTfFaH47uz5e/OHtWuAsYXVGwzM2tR2X02ZHOg7dZgnfnUjFCLiImkmajNzKzDld2MZmZmfcBhY2ZmhXPYmJlZ4Rw2faz6gk0zsyI5bPrQlClTng8aB46ZjQSHTZ9xuJhZGRw2fcRBY2ZlKf06GyuOw8XMOoVrNj3KQWNmncQ1mx7jkDGzTuSaTQ9x0JhZp3LNpgc4ZMys0zlsuphDxsy6hcOmCzlkzKzbOGy6iEPGzLqVBwh0CQeNmXUz12w6nEPGzHqBw6ZDOWTMrJc4bDqMQ8bMepHDpkM4ZMyslzlsSuaQMbN+4NFoJXLQmFm/cM2mBA4ZM+s3DpsR1IkhM2nSpLKLYGZ9wGEzAjoxZMzMRpLDpkAOGTOzxAMECuKgMTNbyjWbnDlkzMyW5bDJiUPGzGxwDps2OWTMzBpz2AyTQ8bMrHkeIDAMDhozs9a4ZtMCh4yZ2fA4bJrgkDEza0/pYSNpPeA7wE6AgMuBwyPiria2PQEYB2wDrA4cEBHT8iqbQ8bMLB+l9tlIWgmYCWwK7A/sC7wauFLSyk3s4tPAisDFeZfNQWNmlp+yazYHAWOAsRFxO4CkG4DbgI8DJzfY/mUR8ZykjYH98iiQQ8bMLH9lh80E4NpK0ABExDxJ1wA70yBsIuK5vArikDEzK07ZYbMFML3O8jnAHiNRAIeMmVnxyg6b1YFFdZY/CAwUcUBJBwMHZ0+XjB079qYijtNn1gDuL7sQPcDnMT8+l/kYm9eOyg4bgKizTIUdLGIqMBVA0uyIGFfUsfqFz2M+fB7z43OZD0mz89pX2TMILCLVbmoNUL/GY2ZmXajssJlD6reptTlw8wiXxczMClJ22MwAtpU0prJA0mhg++y1ok0dgWP0A5/HfPg85sfnMh+5nUdF1OsyGRnZhZv/BJ4Ajib13xwHvBR4bUQ8lq23AXAHcGxEHFu1/duBNYF1gB8APwRmAUTEeSP2RszMbEilhg2ApPV54XQ1V5Cmq5lftc5oYB7wtYiYXLV8FvD2evuNiMIGGZiZWWtKDxszM+t9ZffZ5E7SepLOk/SwpEcknZ/VnprZ9gRJl0l6QFJImlhwcTvWcM+jpHGSpkr6t6THJd0l6WxJG45EuTtRG+dyA0nTJd0p6QlJ90uaJek9I1HuTtPO33bNfo7M/r7/VEQ5O12bn5ExyON1jbbtqbDp5Ik9u0mb53FP0gjD7wPvASYBWwOzsxm++0qb53IV0oWJRwPvBQ4EHgN+J+mDhRW6A+Xwt13ZzxjgKOC+IsrZ6XI6j9OA7WoetzbcKiJ65gEcBjwLbFy1bEPgGeCzTWy/XPZ1Y9JghYllv6duO4/AmnWWbQA8RxrgUfr765ZzOcj+lgfuBi4q+71143kELgVOIw0k+lPZ76vbzmP2uXj8cI7dUzUbBpnYE6hM7DmkyHFizy437PMYEQvrLLsTWAism3M5u0Fbv5O1IuIZ4GHg6dxK2B3aPo+S9ibVso8spITdIdffx1b0WthsAdSb62wO6UJRa06u51HSZsBawL/aLFc3avtcSlpO0vKS1pH0FWAT0jD/ftLWeZQ0QBr1+sWIeDDnsnWTPP62D5G0JOuTnSnprc1s1GthM+ITe/ao3M6jpOWBH5NqNj9tv2hdJ49z+U1STeb/gC8Ce0bEFfkUr2u0ex6/RepXmJZjmbpRu+fxF8AngXeSJjR+OTBT0vhGG3bCRJx5G9GJPXtYXufxFODNwPsiol/nu2v3XH4X+BXp4uX9gHMk7R4R/TaQZVjnMfvPez9g68g6HvrcsH8fI2LfqqdXS5pOqikdD7xlqG17rWbjiT3zkct5lHQi6b+fj0bEZTmVrdu0fS4jYkFEzI6IiyPiQ8C1wEk5lrEbtHMeTyPVqhdIWk3SaqR/tEdlz1fIt6gdLdfPyIh4FLgEeEOjdXstbDyxZz7aPo+SjiINez4sIs7KsWzdpojfydmkEZP9pJ3zuBnwCdKHaeWxPbBt9v0h+RWz4xXx+yjq15ZeoNfCpuyJPXtFW+dR0mdI1eqjIuIHBZWxW+T6OylpOVJzxR05la9btHMed6jz+Cep+WcHoJ/mUcz793FV4H3AXxquXPa475zHkK8M3A7cSBrGN4H0SwZQDG4AAAWASURBVDUXWKVqvQ1I48qPqdn+7cDuwKGkpD4le7572e+tW84j6aLO54Dfk/5zrH5sXvZ767JzOZl0ceyHs9/NDwOXZed3z7LfW7ecx0H2N4v+vM6mnd/HzwOnA3sD40kXhd4IPAW8teGxy37zBZzM9YHfAo8AjwIXAqNr1hmdhcnkmuWzsuXLPMp+X91yHkmjfeqeQ2BW2e+ry87lBNLV3vcBS4A7Sf99bl/2e+qm8zjIvvoybNo5j8AHSNfj3E8aHflA9vv4xmaO64k4zcyscL3WZ2NmZh3IYWNmZoVz2JiZWeEcNmZmVjiHjZmZFc5hY2ZmhXPYmJlZ4Rw2ZmZWOIeN9QVJEyVF1eMpSXdIOkHSSwo65mRJDa+aljRL0qwiypDtfxdJnx3GdmtKOlXS/OxGWbdI+lgRZbTe57CxfrMHsB1p8sBLSbcI/lZBx/pJdqyy7QK0FDbZtPu/A94NfBV4P2m26dMkrZd7Ca3n9eLN08yG8o9Yev/1P0h6NXCgpMMi4rk8DxQRC4AFee5zBO0ObAO8JiLmAEh6ijQJ48plFsy6k2s21u+uB1YE1qheKGkrSTMkLZL0hKRrqu+1LmkTSRdIuk/Sk5LukvSb7DbYlXWWaUaTtKekf2f3cJ8jadd6hWp0/Or9S3q1pEskPSbpTknHZLciQNI00uy861Y1Ic5v4rzsCsyuBE1mN9IkjLc1sb3ZCzhsrN+NBh4mzWALgKStgf8l3dHwINKH7APA5ZK2yVa7GFiXdOOtd5FuFLeEIf6mJL0TOIf0Yf1BUvPd94CxNes1c/xqF5Bmh96FNIPv10gBA3AcqTlsIalJbztSkDSyLfBnSS+WtLGkrwOHAV+OiGeb2N7shcqe7toPP0biAUwkTZk+ltR8PAB8lHTPjkNr1r0C+Bfw4qplo7JlF5JqQQFMaHDMyVTdnoI0PfvNwHJVy95Eze0XGh2/dv/AATXHvRG4rOr5NGBBC+dq7Wy/+wLnsvQWEedWl90PP1p5uM/G+s2/a56fGhGnVJ5IWpF0o7ITgOeqm8WAy4F9SLWMucAUSWuTgmLIpiVJo0j3aZ8SVX1DEfGX6matJo9f65Ka5zcBrx+qPA2My75eB/wdOJs0oOIg4P+Aw9vYt/UpN6NZv9mV9KH/XtKH9ycl7Vf1+uqkWsRXSDeIqn4cSqoRCdiJNDrrROBWSXMlDXUv+zWAFwH/rfNa9bKGx6/0x1R5sOb5EqCd4dzjSDfWuiUiboqIGRHxceCXpLuFmrXMNRvrNzdFNhpN0kzgBuBbkn4bEYuBh0i3Xf4hcGa9HWQ1k7nAfpIEbEUKglMlzY+I39fZrHJ3w7XrvLY26S6ctHD8Io0Dro+I2muEAri34GNbj3LNxvpWRCwBvgCsBXwyW7YYuJoUINdHxOzaR80+IiL+wdLrWLYc5FjPkpqldq+umUh6E2mQQmW9lo7fpCWkEXfN2oaaUMmaC98PnDeM45u5ZmP9LSJmSLoO+LykUyLiCVJw/BG4VNJPSf0UawBbk5q4ziGNIjsXuD1bNpE02GDmEIf7KnAZcKGk04A1SSPHamsLQx4/Iia1+DZvBlbPmvlmA09GxI31VpS0LvAK4EOS5mXlXZd08esC0vs2a5lrNmZwNKl28wmAiLie1K/zAPB90gfu94DXkELgXuAuUijMIPVlvBJ4f0T8bbCDRESlg38scD6pVnU4cEvNeo2O36qfAL8iDTr4K3DREOtWBgecCXyMNMtCJSTHR8Rjwzi+GVq2WdbM+pWk44BPAwN1+mzMhs01GzOrNg74m4PG8uawMbNq25D6dcxy5WY0MzMrnGs2ZmZWOIeNmZkVzmFjZmaFc9iYmVnhHDZmZlY4h42ZmRXOYWNmZoX7/xqTfBCIbAbHAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Parameter values\n", "alpha=1\n", "b=1\n", "gamma=1\n", "q0=0.5\n", "delta=0.5\n", "\n", "# Additional variables/parameters needed\n", "strains=41\n", "IC=[20,5]\n", "fit=[]\n", "\n", "# Trade-off \n", "curve=3\n", "I_c=15*(9-q0*15)/(15*(q0+0.2)-1)\n", "grad=(2)/(3)*I_c\n", "\n", "def resdyn(t,x):\n", " \"\"\"Function to run population dynamics\"\"\"\n", " q=q0*(1+delta*np.sin(2*np.pi*t))\n", " sdot=(a-q*(x[0]+x[1]))*x[0]-b*x[0]-beta*x[0]*x[1]+gamma*x[1]\n", " idot=beta*x[0]*x[1]-(b+alpha+gamma)*x[1]\n", " return [sdot,idot]\n", "\n", "def mutdyn(t,x):\n", " \"\"\"Function to run resident-mutant dynamics\"\"\"\n", " q=q0*(1+delta*np.sin(2*np.pi*t))\n", " sdot=(a-q*(x[0]+x[1]))*x[0]-b*x[0]-beta*x[0]*x[1]+gamma*x[1]\n", " idot=beta*x[0]*x[1]-(b+alpha+gamma)*x[1]\n", " smdot=(am-q*(x[0]+x[1]))*x[2]-b*x[2]-betam*x[2]*x[1]+gamma*x[3]\n", " imdot=betam*x[2]*x[1]-(b+alpha+gamma)*x[3]\n", " return [sdot,idot,smdot,imdot]\n", "\n", "# Loop through resident strains\n", "for res in range(strains):\n", " if res==0:\n", " print(\"Resident strain 1\", end=\" \")\n", " else:\n", " print(\", %i\" %int(res+1), end=\" \")\n", " beta=0.1+0.01*res\n", " a=10 -(grad)**2/(curve)*(1 - np.exp((beta - 0.2)*(curve)/grad))\n", " solr=solve_ivp(resdyn,[ts[0],ts[-1]],IC,t_eval=ts,method='Radau')\n", " \n", " ICm1=[solr.y[0,-1],solr.y[1,-1],1,0]\n", " ICm2=[solr.y[0,-1],solr.y[1,-1],0,1]\n", " \n", " # Loop through mutant strains\n", " for mut in range(strains):\n", " betam=0.1+0.01*mut\n", " am=10 -(grad)**2/(curve)*(1 - np.exp((betam - 0.2)*(curve)/grad))\n", " solm1=solve_ivp(mutdyn,[0,2],ICm1,method='Radau')\n", " solm2=solve_ivp(mutdyn,[0,2],ICm2,method='Radau')\n", " C=[[solm1.y[2,-1],solm1.y[3,-1]],[solm2.y[2,-1],solm2.y[3,-1]]]\n", " eigs=np.linalg.eig(C)\n", " fit.append(np.log(max(eigs[0])))\n", "\n", "# Plotting commands\n", "fita=np.array(fit).reshape(-1, strains)\n", "bplot=np.linspace(0.1,0.5,41)\n", "plt.contourf(bplot,bplot,np.transpose(fita),0,cmap='Greys')\n", "plt.xlabel('Resident $\\\\beta$')\n", "plt.ylabel('Mutant $\\\\beta_m$')\n", "plt.title('Pairwise Invasion Plot')\n", "plt.tight_layout()\n", "#plt.savefig('pip2.png')" ] }, { "cell_type": "markdown", "id": "3f05f500", "metadata": {}, "source": [ "The next code cell builds on the previous cell to find a CSS point as a parameter is varied by considering the invasion fitness of a nearby mutant." ] }, { "cell_type": "code", "execution_count": null, "id": "c0aa13ce", "metadata": {}, "outputs": [], "source": [ "# Parameter values\n", "alpha=1\n", "b=1\n", "gamma=1\n", "q0=0.5\n", "delta=0.5\n", "\n", "# Additional variables needed\n", "strains=41\n", "IC=[20,5]\n", "fit=[]\n", "\n", "def resdyn(t,x):\n", " \"\"\"Function to run population dynamics\"\"\"\n", " q=q0*(1+delta*np.sin(2*np.pi*t))\n", " sdot=(a-q*(x[0]+x[1]))*x[0]-b*x[0]-beta*x[0]*x[1]+gamma*x[1]\n", " idot=beta*x[0]*x[1]-(b+alpha+gamma)*x[1]\n", " return [sdot,idot]\n", "\n", "def mutdyn(t,x):\n", " \"\"\"Function to run resident-mutant dynamics\"\"\"\n", " q=q0*(1+delta*np.sin(2*np.pi*t))\n", " sdot=(a-q*(x[0]+x[1]))*x[0]-b*x[0]-beta*x[0]*x[1]+gamma*x[1]\n", " idot=beta*x[0]*x[1]-(b+alpha+gamma)*x[1]\n", " smdot=(am-q*(x[0]+x[1]))*x[2]-b*x[2]-betam*x[2]*x[1]+gamma*x[3]\n", " imdot=betam*x[2]*x[1]-(b+alpha+gamma)*x[3]\n", " return [sdot,idot,smdot,imdot]\n", "\n", "# Vary the baseline competition\n", "solstore=[]\n", "for q0 in [0.1,0.5]:\n", " # Trade-off\n", " curve=-3\n", " I_c=15*(9-q0*15)/(15*(q0+0.2)-1)\n", " grad=(2)/(3)*I_c\n", " \n", " for d in range(21):\n", " delta=d*0.05\n", " if d==0:\n", " print(\"Delta = 0\", end=\" \")\n", " else:\n", " print(\", %.2f\" %(delta), end=\" \")\n", " \n", " for res in range(strains):\n", " beta=0.1+0.01*res\n", " a=10 -(grad)**2/(curve)*(1 - np.exp((beta - 0.2)*(curve)/grad))\n", " solr=solve_ivp(resdyn,[ts[0],ts[-1]],IC,t_eval=ts,method='Radau')\n", "\n", " ICm1=[solr.y[0,-1],solr.y[1,-1],1,0]\n", " ICm2=[solr.y[0,-1],solr.y[1,-1],0,1]\n", "\n", " betam=beta+0.01\n", " am=10 -(grad)**2/(curve)*(1 - np.exp((betam - 0.2)*(curve)/grad))\n", " solm1=solve_ivp(mutdyn,[0,2],ICm1,method='Radau')\n", " solm2=solve_ivp(mutdyn,[0,2],ICm2,method='Radau')\n", " C=[[solm1.y[2,-1],solm1.y[3,-1]],[solm2.y[2,-1],solm2.y[3,-1]]]\n", " eigs=np.linalg.eig(C)\n", " if(np.log(max(eigs[0])))<0:\n", " solstore.append(beta)\n", " break\n", "\n", "#Plotting commands\n", "delt=np.linspace(0,1,21)\n", "plt.scatter(delt,solstore[0:21],color='b',label='$q=0.1$')\n", "plt.scatter(delt,solstore[21:42],color='r',label='$q=0.5$')\n", "plt.xlabel('Amplitude, $\\\\delta$')\n", "plt.ylabel('CSS Transmission, $\\\\beta$')\n", "plt.ylim([0.1, 0.5])\n", "plt.tight_layout()\n", "plt.legend()\n", "#plt.savefig('varyalpha.png')" ] }, { "cell_type": "markdown", "id": "6e51e295", "metadata": {}, "source": [ "## Case study 2\n", "\n", "Now we move on to a model with a free-living parasite stage, that also leads to limit cycles for some (but not all) parameter values. Noitice the dynamics have bene log-transformed to prevent the model failing to run due to stiffness." ] }, { "cell_type": "code", "execution_count": 1, "id": "1249917d", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:20: RuntimeWarning: overflow encountered in exp\n", "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:20: RuntimeWarning: invalid value encountered in double_scalars\n", "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:21: RuntimeWarning: overflow encountered in exp\n", "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:22: RuntimeWarning: overflow encountered in exp\n" ] } ], "source": [ "# Import libraries\n", "from scipy.signal import find_peaks\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.integrate import solve_ivp\n", "plt.rcParams.update({'font.size':16})\n", "\n", "# Parameter values\n", "alpha=1\n", "b=1\n", "gamma=0.1\n", "q=0.1\n", "mu=0.1\n", "theta=5\n", "beta=0.1\n", "a=10\n", "\n", "def logdyn(t,x):\n", " \"\"\"Function to run population dynamics\"\"\"\n", " sdot=a-q*(np.exp(x[0])+np.exp(x[1]))-b-beta*np.exp(x[2])+gamma*np.exp(x[1]-x[0])\n", " idot=beta*np.exp(x[0]+x[2]-x[1])-(b+alpha+gamma)\n", " pdot=theta*np.exp(x[1]-x[2])-mu\n", " return [sdot,idot,pdot]\n", "\n", "# Plot commands\n", "ts=np.linspace(0,500,50001)\n", "sol=solve_ivp(logdyn,[ts[0],ts[-1]],[-15,-10,5],t_eval=ts,method='Radau')\n", "plt.plot(sol.t,np.exp(sol.y[1]),color='r',label='Infected')\n", "plt.legend()\n", "plt.xlabel('Time')\n", "plt.ylabel('Density')\n", "plt.title('Resident dynamics')\n", "\n", "# Add crosses for peaks\n", "peakcheck=np.exp(sol.y[1])\n", "peaks,_=find_peaks(peakcheck,height=1)\n", "plt.plot(ts[peaks], (peakcheck[peaks]), \"x\")\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "id": "c75f31f4", "metadata": {}, "source": [ "Now we allow the host to evolve avoidance but with a trade-off with reproduction. Again we demonstrate the results trhough a pairwise invasion plot. *Note: this system is slower to run than the previous one: have patience!*" ] }, { "cell_type": "code", "execution_count": null, "id": "f43ae790", "metadata": {}, "outputs": [], "source": [ "# Parameter values\n", "alpha=1\n", "b=1\n", "gamma=0.1\n", "q=0.1\n", "mu=0.1\n", "theta=5\n", "\n", "strains=41\n", "IC=[-5,-3.5,2]\n", "fit=[]\n", "ts=np.linspace(0,5000,50000)\n", "\n", "# Trade-off\n", "curve=-400\n", "grad=75\n", "\n", "def resdyn(t,x):\n", " \"\"\"Function to run population dynamics\"\"\"\n", " sdot=a-q*(np.exp(x[0])+np.exp(x[1]))-b-beta*np.exp(x[2])+gamma*np.exp(x[1]-x[0])\n", " idot=beta*np.exp(x[0]+x[2]-x[1])-(b+alpha+gamma)\n", " pdot=theta*np.exp(x[1]-x[2])-mu\n", " return [sdot,idot,pdot]\n", "\n", "def mutdyn(t,x):\n", " \"\"\"Function to run mutant-residentpopulation dynamics\"\"\"\n", " sdot=a-q*(np.exp(x[0])+np.exp(x[1]))-b-beta*np.exp(x[2])+gamma*np.exp(x[1]-x[0])\n", " idot=beta*np.exp(x[0]+x[2]-x[1])-(b+alpha+gamma)\n", " pdot=theta*np.exp(x[1]-x[2])-mu\n", " \n", " smdot=am-q*(np.exp(x[0])+np.exp(x[1]))-b-betam*np.exp(x[2])+gamma*np.exp(x[4]-x[3])\n", " imdot=betam*np.exp(x[3]+x[2]-x[4])-(b+alpha+gamma)\n", " return [sdot,idot,pdot,smdot,imdot]\n", "\n", "# Loop through resident strains\n", "for res in range(strains):\n", " beta=0.05+0.0025*res\n", " a=10 -(grad)**2/(curve)*(1 - np.exp((beta - 0.1)*(curve)/grad))\n", " solr=solve_ivp(resdyn,[ts[0],ts[-1]],IC,t_eval=ts,method='Radau')\n", " if res==0:\n", " print(\"Resident strain 1\", end=\" \")\n", " else:\n", " print(\", %i\" %int(res+1), end=\" \")\n", " \n", " # Use peaks to find period\n", " peakcheck=np.exp(solr.y[1])\n", " peaks,_=find_peaks(peakcheck,height=0.2)\n", " if len(peaks)<1 or ts[peaks[-1]]-ts[peaks[-2]]<1:\n", " period=1\n", " else:\n", " period=ts[peaks[-1]]-ts[peaks[-2]]\n", " \n", " ICm1=[solr.y[0,peaks[-2]],solr.y[1,peaks[-2]],solr.y[2,peaks[-2]],0,-100]\n", " ICm2=[solr.y[0,peaks[-2]],solr.y[1,peaks[-2]],solr.y[2,peaks[-2]],-100,0]\n", " \n", " # Loop through mutant strains\n", " for mut in range(strains):\n", " betam=0.05+0.0025*mut\n", " am=10 -(grad)**2/(curve)*(1 - np.exp((betam - 0.1)*(curve)/grad))\n", " solm1=solve_ivp(mutdyn,[0,period*2],ICm1,method='Radau')\n", " solm2=solve_ivp(mutdyn,[0,period*2],ICm2,method='Radau')\n", " C=[[np.exp(solm1.y[3,-1]),np.exp(solm1.y[4,-1])],[np.exp(solm2.y[3,-1]),np.exp(solm2.y[4,-1])]]\n", " eigs=np.linalg.eig(C)\n", " fit.append(np.log(max(eigs[0])))\n", "\n", "# Plotting commands\n", "fita=np.array(fit).reshape(-1, strains)\n", "from numpy import inf\n", "fita[fita == -inf] = -100\n", "bplot=np.linspace(0.05,0.15,41)\n", "plt.contourf(bplot,bplot,np.transpose(fita),0,cmap='Greys')\n", "plt.xlabel('Resident $\\\\beta$')\n", "plt.ylabel('Mutant $\\\\beta_m$')\n", "plt.title('Pairwise Invasion Plot')\n", "#plt.savefig('pipsip.png')" ] }, { "cell_type": "markdown", "id": "7b1b283a", "metadata": {}, "source": [ "Again, our next stage is to vary a parameter value and then identify the CSS." ] }, { "cell_type": "code", "execution_count": null, "id": "67dde3b0", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Parameter values\n", "alpha=1\n", "b=1\n", "gamma=0.1\n", "q=0.1\n", "mu=0.1\n", "theta=4\n", "\n", "strains=41\n", "IC=[-5,-3.5,2]\n", "fit=[]\n", "ts=np.linspace(0,1000,50000)\n", "\n", "# Trade-off\n", "curve=-400\n", "grad=75\n", "\n", "def resdyn(t,x):\n", " \"\"\"Function to run population dynamics\"\"\"\n", " sdot=a-q*(np.exp(x[0])+np.exp(x[1]))-b-beta*np.exp(x[2])+gamma*np.exp(x[1]-x[0])\n", " idot=beta*np.exp(x[0]+x[2]-x[1])-(b+alpha+gamma)\n", " pdot=theta*np.exp(x[1]-x[2])-mu\n", " return [sdot,idot,pdot]\n", "\n", "def mutdyn(t,x):\n", " \"\"\"Function to run mutant-residentpopulation dynamics\"\"\"\n", " sdot=a-q*(np.exp(x[0])+np.exp(x[1]))-b-beta*np.exp(x[2])+gamma*np.exp(x[1]-x[0])\n", " idot=beta*np.exp(x[0]+x[2]-x[1])-(b+alpha+gamma)\n", " pdot=theta*np.exp(x[1]-x[2])-mu\n", " \n", " smdot=am-q*(np.exp(x[0])+np.exp(x[1]))-b-betam*np.exp(x[2])+gamma*np.exp(x[4]-x[3])\n", " imdot=betam*np.exp(x[3]+x[2]-x[4])-(b+alpha+gamma)\n", " return [sdot,idot,pdot,smdot,imdot]\n", "\n", "cssstore=[]\n", "for d in range(21):\n", " theta=4.0+d*0.2\n", " if d==0:\n", " print(\"Theta = 4.0\", end=\" \")\n", " else:\n", " print(\", %.1f\" %(theta), end=\" \")\n", " for res in range(strains):\n", " beta=0.05+0.0025*res\n", " a=10 -(grad)**2/(curve)*(1 - np.exp((beta - 0.1)*(curve)/grad))\n", " solr=solve_ivp(resdyn,[ts[0],ts[-1]],IC,t_eval=ts,method='Radau')\n", " peaks,_=find_peaks(np.exp(solr.y[1]),height=0.2)\n", " if len(peaks)<1 or ts[peaks[-1]]-ts[peaks[-2]]<1:\n", " period=1\n", " else:\n", " period=ts[peaks[-1]]-ts[peaks[-2]]\n", " \n", " ICm1=[solr.y[0,-1],solr.y[1,-1],solr.y[2,-1],0,-100]\n", " ICm2=[solr.y[0,-1],solr.y[1,-1],solr.y[2,-1],-100,0]\n", "\n", " betam=beta+0.005\n", " am=10 -(grad)**2/(curve)*(1 - np.exp((betam - 0.1)*(curve)/grad))\n", " solm1=solve_ivp(mutdyn,[0,period*2],ICm1,method='Radau')\n", " solm2=solve_ivp(mutdyn,[0,period*2],ICm2,method='Radau')\n", " C=[[np.exp(solm1.y[3,-1]),np.exp(solm1.y[4,-1])],[np.exp(solm2.y[3,-1]),np.exp(solm2.y[4,-1])]]\n", " eigs=np.linalg.eig(C)\n", " if(np.log(max(eigs[0])))<0:\n", " cssstore.append(beta)\n", " break\n", " \n", "# Plotting commands\n", "alpha=np.linspace(3,5,21)\n", "plt.scatter(alpha,cssstore[0:21],color='b')\n", "plt.xlabel('Virulence, $\\\\alpha$')\n", "plt.ylabel('CSS Transmission, $\\\\beta$')\n", "plt.ylim([0.05, 0.15])\n", "plt.tight_layout()\n", "#plt.savefig('varyalpha.png')" ] }, { "cell_type": "markdown", "id": "af12d588", "metadata": {}, "source": [ "For this model we also want to know the period of the cycles. We calculate this below. We note that sometimes an equilibrium can show up as being a very slow cycle. Some trial and error is needed to spot these values and account for them." ] }, { "cell_type": "code", "execution_count": null, "id": "1a675131", "metadata": {}, "outputs": [], "source": [ "# Parameter values\n", "alpha=1\n", "b=1\n", "gamma=0.1\n", "q=0.1\n", "mu=0.1\n", "theta=5\n", "\n", "strains=41\n", "IC=[-5,-3.5,2]\n", "ts=np.linspace(0,1000,50000)\n", "\n", "# Trade-off\n", "curve=-400\n", "grad=75\n", "\n", "def resdyn(t,x):\n", " \"\"\"Function to run population dynamics\"\"\"\n", " sdot=a-q*(np.exp(x[0])+np.exp(x[1]))-b-beta*np.exp(x[2])+gamma*np.exp(x[1]-x[0])\n", " idot=beta*np.exp(x[0]+x[2]-x[1])-(b+alpha+gamma)\n", " pdot=theta*np.exp(x[1]-x[2])-mu\n", " return [sdot,idot,pdot]\n", "\n", "solstore=[]\n", "\n", "for d in range(21):\n", " theta=4.0+d*0.2\n", " if d==0:\n", " print(\"Theta = 4.0\", end=\" \")\n", " else:\n", " print(\", %.1f\" %(theta), end=\" \")\n", " for res in range(strains):\n", " beta=0.05+0.0025*res\n", " a=10 -(grad)**2/(curve)*(1 - np.exp((beta - 0.1)*(curve)/grad))\n", " sols=solve_ivp(resdyn,[ts[0],ts[-1]],IC,t_eval=ts,method='Radau')\n", " peaks,_=find_peaks(np.exp(sols.y[1]),height=0.2)\n", " if len(peaks)<1 or ts[peaks[-1]]-ts[peaks[-2]]<1:\n", " period=1\n", " else:\n", " period=ts[peaks[-1]]-ts[peaks[-2]]\n", " solstore.append(period)\n", " \n", "for i in range(861):\n", " if solstore[i]>20:\n", " solstore[i]=1\n", "\n", "plt.contour(np.array(solstore).reshape(21,41).transpose())" ] }, { "cell_type": "markdown", "id": "957c31f1", "metadata": {}, "source": [ "Finally we put the last two results together, plotting the CSS points on top of the periods." ] }, { "cell_type": "code", "execution_count": null, "id": "c72bf6a1", "metadata": {}, "outputs": [], "source": [ "#plt.rcParams[\"pcolor.shading\"]='nearest'\n", "q=np.linspace(4,8,21)\n", "beta=np.linspace(0.05,0.15,41)\n", "plt.pcolormesh(q,beta,np.array(solstore).reshape(21,41).transpose(),shading='gouraud')\n", "plt.colorbar(label='Period')\n", "plt.scatter(q,cssstore,color='white')\n", "plt.xlabel('Production, $\\\\theta$')\n", "plt.ylabel('CSS Transmission, $\\\\beta$')\n", "plt.tight_layout()\n", "#plt.savefig('varytheta_sips.png')" ] }, { "cell_type": "code", "execution_count": null, "id": "d960a195", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.4" } }, "nbformat": 4, "nbformat_minor": 5 }