example-8-multi-model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "345e05bd",
   "metadata": {},
   "source": [
    "## Example 8: Multi-model simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5d7101a",
   "metadata": {},
   "source": [
    "In this example we show how to work with multiple models using the same code for each.\n",
    "\n",
    "To let the script know which variables have particular meanings (e.g. which is the membrane potential), we've added _labels_ (for outputs) and _bindings_ (a silly name for inputs) to each model.\n",
    "\n",
    "We place the following requirements on the models:\n",
    "\n",
    "- Each must have a meta data property `name` with a simple model name.\n",
    "- Each must have its membrane potential variable labelled as `membrane_potential`.\n",
    "- Each must have its L-type calcium current (all species and compartments combined) labelled as `ICaL`.\n",
    "- Each must contain a variable labelled `gCaL`, that can be scaled to scale `ICaL`.\n",
    "- If a model has modes/switches, it must be in baseline mode (no disease models etc.) and set to epicardial cells (if possible).\n",
    "\n",
    "We also assume that each `mmt` file contains"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4117fd75",
   "metadata": {},
   "source": [
    "Using the assumptions, in the example below we:\n",
    "\n",
    "- iterate over a list of model file names\n",
    "- load the model and protocol from each file\n",
    "- obtain the membrane potential and ICaL variable from each model\n",
    "- run a 1Hz simulation to obtain the baseline action potential duration (APD)\n",
    "- modify the ICaL in the model to be half its size, run another simulation and measure the APD again\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c1ec5943",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loaded priebe-1998\n",
      "  Pre-pacing...\n",
      "  Measuring baseline apd...\n",
      "  Measuring modified apd...\n",
      "  Difference: -57.03690247866763 ms\n",
      "Loaded iyer-2004\n",
      "  Pre-pacing...\n",
      "  Measuring baseline apd...\n",
      "  Measuring modified apd...\n",
      "  Difference: -90.15790640803968 ms\n",
      "Loaded tentusscher-2004\n",
      "  Pre-pacing...\n",
      "  Measuring baseline apd...\n",
      "  Measuring modified apd...\n",
      "  Difference: -55.342723760400816 ms\n",
      "Loaded tentusscher-2006\n",
      "  Pre-pacing...\n",
      "  Measuring baseline apd...\n",
      "  Measuring modified apd...\n",
      "  Difference: -50.4367058138055 ms\n",
      "Loaded priebe-1998\n",
      "  Pre-pacing...\n",
      "  Measuring baseline apd...\n",
      "  Measuring modified apd...\n",
      "  Difference: -57.03690247866763 ms\n",
      "Loaded grandi-2010\n",
      "  Pre-pacing...\n",
      "  Measuring baseline apd...\n",
      "  Measuring modified apd...\n",
      "  Difference: -58.90907637623076 ms\n",
      "Loaded ohara-2011\n",
      "  Pre-pacing...\n",
      "  Measuring baseline apd...\n",
      "  Measuring modified apd...\n",
      "  Difference: -27.072142655352792 ms\n",
      "Loaded li-2017\n",
      "  Pre-pacing...\n",
      "  Measuring baseline apd...\n",
      "  Measuring modified apd...\n",
      "  Difference: -36.58475121189454 ms\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "\n",
    "import myokit\n",
    "\n",
    "models = [\n",
    "    'priebe-1998.mmt',\n",
    "    'iyer-2004.mmt',\n",
    "    'tentusscher-2004.mmt',\n",
    "    'tentusscher-2006.mmt',\n",
    "    'priebe-1998.mmt',\n",
    "    'grandi-2010.mmt',\n",
    "    'ohara-2011.mmt',\n",
    "    'ohara-cipa-v1-2017.mmt',\n",
    "]\n",
    "\n",
    "# This will be used to store the results\n",
    "names = []\n",
    "delta_apds = []\n",
    "\n",
    "# Iterate over the list of models\n",
    "for filename in models:\n",
    "    path = os.path.join('models', 'c', filename)\n",
    "    model, protocol, _ = myokit.load(path)\n",
    "    \n",
    "    # Show the loaded model's name\n",
    "    name = model.meta['name']\n",
    "    names.append(name)\n",
    "    print('Loaded ' + name)\n",
    "\n",
    "    # Find the membrane potential, the l-type calcium current,\n",
    "    # and a variable that scales it.\n",
    "    # Instead of using .label, we use .labelx, which will raise\n",
    "    # an exception if model doesn't contain the label.\n",
    "    vm = model.labelx('membrane_potential')\n",
    "    ical = model.labelx('ICaL')\n",
    "    gcal = model.labelx('gCaL')\n",
    "\n",
    "    # Obtain the duration of a stimulus pulse from the embedded protocol\n",
    "    duration = protocol.head().duration()\n",
    "    # ...and use it to create a new protocol with cycle length 1000ms\n",
    "    protocol = myokit.pacing.blocktrain(1000, duration)\n",
    " \n",
    "    # Create a simulation\n",
    "    s = myokit.Simulation(model, protocol)\n",
    "\n",
    "    # Pre-pace the model for 50 beats\n",
    "    print('  Pre-pacing...')\n",
    "    #s.pre(50 * 1000)\n",
    "   \n",
    "    # Measure the baseline APD at V=-70mV\n",
    "    print('  Measuring baseline apd...')\n",
    "    x, apds = s.run(1000, log=myokit.LOG_NONE, apd_variable=vm, apd_threshold=-70)\n",
    "\n",
    "    # The returned object `apds` contains more info that we need, extract only\n",
    "    # the duration of the first AP\n",
    "    apd_baseline = apds['duration'][0]\n",
    "\n",
    "    # Update the model, reducing ICaL to 50%\n",
    "    ical.set_rhs(myokit.Multiply(ical.rhs(), myokit.Number(0.5)))\n",
    "    \n",
    "    # Scale the L-type calcium current\n",
    "    s.set_constant(gcal, 0.5 * gcal.rhs().value())\n",
    "    \n",
    "    # Measure the new apd\n",
    "    print('  Measuring modified apd...')\n",
    "    s.reset()\n",
    "    x, apds = s.run(\n",
    "        1000, log=myokit.LOG_NONE, apd_variable=vm, apd_threshold=-70)\n",
    "    apd_modified = apds['duration'][0]\n",
    "\n",
    "    # Store the difference between the APDs\n",
    "    delta = apd_modified - apd_baseline\n",
    "    delta_apds.append(delta)\n",
    "    print(f'  Difference: {delta} ms')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a61d6fb0",
   "metadata": {},
   "source": [
    "Next, we create a bar graph of the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c27deed1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhoAAAIECAYAAABfZYfNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABGK0lEQVR4nO3deZxkVX3//9d7hmZ0FBFQAq7jElGCEQOCqCijDphg/LokXxU3CNFEf3HFfWMUFcQlaIz6FTUoEZW4oTJhGAyoEdSwKURUQAZEFtkHGBjamc/vj3tLaopeqqHvdE/P6/l41KOr7j33nlOnu+u+69xzq1JVSJIkdWHeTDdAkiTNXQYNSZLUGYOGJEnqjEFDkiR1xqAhSZI6Y9CQJEmdMWhIkqTOGDQkSVJnNpvpBmjjkiTA/YAbZ7otkqQZtwVwWU3w6Z8GDU3V/YBLZ7oRkqRZ4wHA78ZbadDQVN0IcNFFF7H11lvPdFtmrdHRUU488UT23ntvRkZGZro5s5b9NBz7aXL20XCms59WrVrFAx/4QJhkhNugoTtliy224F73utdMN2PWGh0dZeHChdzrXvfyRW8C9tNw7KfJ2UfDmYl+cjKoJEnqjEFDkiR1xqAhSZI6Y9CQJEmdMWhIkqTOGDQkSVJnDBqSJKkzBg1JktQZg4YkSeqMQUOSJHXGoCFJkjpj0JAkSZ0xaEiSpM4YNCRJUmcMGpIkqTMGDUmS1BmDhiRJ6oxBQ5IkdcagIUmSOmPQkCRJnTFoSJKkzhg0JElSZwwakiSpMwYNSZLUGYOGJEnqjEFDkiR1xqAhSZI6s9lMN0CSNP0WvfX4mW7CBrVgfnH4brDT0uWsWZuZbs6s1eunDckRDUmS1BmDhiRJ6oxBQ5IkdcagIUmSOmPQkCRJnTFoSJKkzhg0JElSZwwakiSpMwYNSZLUGYOGJEnqjEFDkiR1xqAhSZI6Y9CQJEmdMWhIkqTOGDSmWZJTkhzR3l+Z5HUz2yJJkmaOQaNbjwM+M97KJIuSfC7JRUluSXJhkvck2Xyg3IOSfCfJzUmuTvLx/jJJ7pbkqCTnJPlDkm+NUddRSWqM2/9O5xOWJKnfZjPdgLmsqq6apMgjacLePwAXADsBRwL3AN4IkGQ+cDxwFfAkYBvgC0CAV7f7mQ/cAnwceN44db0WeGvf482AnwH/MZXnJEnSVBg0OpRkJXBEVR0x1vqqOgE4oW/Rb5LsALySNmgAewM7Ag+sqsva/R4EHJXkHVW1qqpubrchyROBe49R1w3ADX1tezawFfBvd/4ZSpI0MYPG7LMlcG3f4z2Ac3sho7UcWADsApx8J+s5EDipqi6eqFCSBW1dPVsAjI6OMjo6eiernvt6fWMfTcx+Gs6d6acF86ur5sxKC+bVej81tl7/TMf/3LD7MGjMIkkeRnM65KC+xdsBV/aXq6rrktzWrrsz9WwP/CWw3xDF3wYcPLjw5JNPZuHChXem+k3KihUrZroJGwX7aThT6afDd+uwIbPYIbuum+kmbBSm439u9erVQ5UzaGwgST4NvLj3uKruObD+fjSnUf6jqj47sPlYET3jLB/G/sD1wLeGKHso8NG+x1sAly5evJhtttnmTlY/942OjrJixQqWLFnCyMjITDdn1rKfhnNn+mmnpcs7btXssmBecciu63jX6fNYsy4z3ZxZq9dP0/E/t2rVqqHKGTQ2nHcDHx5rRRsyTgZOA14xsPoKYPeB8lsBIwyMdAwjSYC/A46uqtsmK19Va4A1fdsDMDIy4oFhCPbTcOyn4Uyln9as3TQPtmvWZZN97lMxHf9zw27v5a0bSFX9vqou6N16y5PcHzgFOBM4oKoGx/1OA3ZqT3f07E1z8D/jTjTlKcDDgc/diW0lSZoSRzRmUDuScQpwCc1VJvftjRhU1RVtsROBXwBHJ3kTsDXNyMiRVbWqb187Apu367dIsnO7n7MHqj0Q+ElVndvJk5IkqY9BY2btTTO68HDg0oF1AaiqtUn2BT4J/Ijm8zKO4fbLX3uWAQ/ue3xW/34AkmxJ8zkbr52m9kuSNCGDxjSrqr367i+apOxRwFFD7PMS4JmTlJmwrrbMDYCXikiSNhjnaEiSpM4YNCRJUmcMGpIkqTMGDUmS1BmDhiRJ6oxBQ5IkdcagIUmSOmPQkCRJnTFoSJKkzhg0JElSZwwakiSpMwYNSZLUGYOGJEnqjN/eKklz0MrD9p3pJmxQo6OjLFu2jHOX7sPIyMhMN2fW6vXThuSIhiRJ6oxBQ5IkdcagIUmSOmPQkCRJnTFoSJKkzhg0JElSZwwakiSpMwYNSZLUGYOGJEnqjJ8MKknq1KK3Ht95HQvmF4fvBjstXc6atem8vo1Vr582JEc0JElSZwwakiSpMwYNSZLUGYOGJEnqjEFDkiR1xqAhSZI6Y9CQJEmdMWhIkqTOGDQkSVJnDBqSJKkzBg1JktQZg4YkSeqMQUOSJHXGoCFJkjqzyQWNJHslqST3num2SJI0121yQWM2SvLyJD9Mcl17OynJbmOUe1WSi5LcmuSMJHsOrH9ukuVJrm7D1M5j7OMVSU5JssrAJUnqmkFjmiTZ/C5svhfwZWAxsAdwCXBikvv37f/5wBHA+4HHAj8E/jPJg/r2cw/gR8BbJ6hrIXAC8IG70F5JkoYyJ4NGkgVJPp7k9+27//9O8riBYrskOT3J6iSnJtmhb/uHJTkuyZVJbkryP0mePlDHyiTvTHJUkhuAI9vlH0zy63a/v0lySJKRidpbVS+qqk9W1dlV9Uvg5TS/m6f1FXsD8Lmq+mxVnVdVrwN+C7yybz9HV9V7gZMmqOuIqjoM+PFEbZIkaTpsNtMN6MjhwPOAlwEXA28Glid5eF+Z9wMHAVcBnwY+DzyxXXdPYBnwTuDWdj/fSbJDVV3St483AYcA7+tbdiOwP3AZ8GiaAHJj26ZhLQRGgGvhj6MluwCHDZQ7EXjCFPY7ZUkWAAv6Fm0BMDo6yujoaJdVb9R6fWMfTcx+Gs7G3k8L5lf3dcyr9X5qbL3+mY6/pWH3kaq59UtJcg/gOmD/qjqmXTYCrKQ59fA/wMnA06vqe+36vwKOB+5eVbeOs9//BT5VVZ9oH68Ezqqq50zSnjcBz6+qXafwHP4V2AfYqapuTXI/4HfAE6vq1L5ybwdeVlU7DGy/CLgIeGxVnT1OHXvR9MNWVXX9BG1ZChw8uPyYY45h4cKFwz4lSdIcs3r1avbbbz+ALatq1Xjl5uKIxsNoRgN+1FtQVaNJfgo8iiZoAPy8b5vL25/bApe0YeVg4JnA/Wj66e5A/3wIgNMHK0/yN8DrgIfTjIxsBqxq1z0I+EVf8Q9U1QcGtn8z8EJgrzFCz2AqzBjLptuhwEf7Hm8BXLp48WK22WabjqveeI2OjrJixQqWLFnCyMiEZ842afbTcDb2ftpp6fLO61gwrzhk13W86/R5rFmXzuvbWPX6aTr+llatGjdbrGcuBo3eX9hkB+X+MZ/e8t6clQ/RjCi8EbgAuAX4GjA44fPm9SpIHg98hSakLAduAF5Ac4oGmtMpO/dtcu3A9m8E3k4z2tIfhK4G1gLbDdS/LXAlHaqqNcCavjYCMDIyslG+4G1o9tNw7KfhbKz9tGbthjvwr1mXDVrfxmo6/paG3X4uTga9ALgNeFJvQXvqZFfgvCH3sSdwVFV9s6rOAa4AFg2x3ROBi6vq/VV1elWdDzy4t7Kq/lBVF/Td/hg02lMs7wKeUVXrjZRU1W3AGcCSgfqWAKciSdIsNedGNKrq5iSfAj6U5FqaS0XfTDPB8nPAY4bYzQXAc5N8h2a04xCGC2UXAA9K8gKaUzT7AhPO4YA/ni45BNgPWJmkN3JxU1Xd1N7/KHB0ktOB04BX0JzK+XTffrZul92vXbRDOwJxRVVd0ZbZjmZkpDcx9tFJbgQu6Q8+kiRNh7k4ogHN50h8HTgaOJPmoLpPVV035Pavp5lQeirwHZrTIGdOtlFVHQf8M/AJ4GyaK0IOGaK+V9GclvkazXyR3u2Nffv+Ks3cj3e3+34y8FdVdXHffp4FnEUzsRWa0zhnAf/YV+Yf22VHto9/0D5+1hDtlCRpSubciAZAO4nyNe1tcN0p3D6Po7fs7P5lVbUSeOrApv86sM2icep+M80ISr8jJmnvmPsao9wngU9OsP4o4KhJ9rEUWDpMfZIk3VVzdURDkiTNAgYNSZLUGYOGJEnqjEFDkiR1xqAhSZI6Y9CQJEmdMWhIkqTOGDQkSVJnDBqSJKkzBg1JktQZg4YkSeqMQUOSJHXGoCFJkjpj0JAkSZ2Zk18TL0maPVYetm/ndYyOjrJs2TLOXboPIyMjnde3ser104bkiIYkSeqMQUOSJHXGoCFJkjpj0JAkSZ0xaEiSpM4YNCRJUmcMGpIkqTMGDUmS1BmDhiRJ6oyfDCppVln01uNnugmzzoL5xeG7wU5Ll7NmbWa6ObOSfTScXj9tSI5oSJKkzhg0JElSZwwakiSpMwYNSZLUGYOGJEnqjEFDkiR1xqAhSZI6Y9CQJEmdMWhIkqTOGDQkSVJnDBqSJKkzBg1JktQZg4YkSeqMQUOSJHXGoDHNkpyS5Ii+xyuTvG7mWiRJ0swxaHTvccBnxluZZFGSzyW5KMktSS5M8p4kmw+Ue1CS7yS5OcnVST7eXybJ3ZIcleScJH9I8q1x6ntKkjOS3JrkN0n+cbqeqCRJgzab6QZsaEk2r6rbNlR9VXXVJEUeSRP4/gG4ANgJOBK4B/BGgCTzgeOBq4AnAdsAXwACvLrdz3zgFuDjwPPGqijJQ4Bl7f5fDDwR+GSSq6rq63fuGUqSNL6NfkQjyRZJvtS+0788yev7T1+0py7e2b7bv4HmIEuSDyb5dZLV7Tv7Q5KM9O13aZKzk7yk3ccNSb6SZIu+MvdI8sUkN7V1HzRG+yY8dVJVJ1TVAVV1YlX9pqq+DXwYeG5fsb2BHYEXV9VZVXUScBDw8iT3avdzc1W9sqqOBK4Yp7p/BC6pqtdV1XlV9Vng87SBRpKk6TYXRjQ+SvPO/FnAlcB7gb8Azu4r8ybgEOB9fctuBPYHLgMeTRNAbgQO7yvzMODZwDOBrYBjgbcC72jXfwhYDDyH5uD+AWCXgbrvjC2Ba/se7wGcW1WX9S1bDixo6zt5yP3uAZw4sGw5cGCSkaoaHdwgyYK2np4tAEZHRxkdvUNxtXp9Yx9NbKx+WjC/Zqo5s9aCebXeT92RfTScXv9Mx2vTsPvYqINGO7rwMmC/qvpeu+wAmvDQ77+q6sP9C6qqP3SsTPIR4PmsHzTmAftX1Y3tvo8Gnga8I8k9gQOBl1bVinb9y4BL7+JzehjN6ZD+0ZHtaEJUf/uvS3Jbu25Yd9hP+3gz4D7A5WNs8zbg4MGFJ598MgsXLpxC1ZumFStWzHQTNgr9/XT4bjPYkFnukF3XzXQTZj37aDjT8dq0evXqocpt1EEDeCgwAvy0t6Cqbkjyq4Fypw9umORvgNcBDwfuSdMXqwaKreyFjNblwLbt/YcBmwOn9dV97Rh199f5aZq5Eb3y9xxYfz/gBOA/2tMa/caK6Rln+UQGy2eC/QMcSjNq1LMFcOnixYvZZpttplj1pmN0dJQVK1awZMkSRkZGJt9gEzVWP+20dPkMt2r2WTCvOGTXdbzr9HmsWZfJN9gE2UfD6fXTdLw2rVo1eMgc28YeNMY7SA7+ld283srk8cBXaN6pLwduAF7A+qMIAIPjQsXt81ruzF/yu2nmX9xBGzJOpgkurxhYfQWw+0D5rWhC1uAIxUSu4I4jINsCfwCuGWuDqloDrOmrF4CRkREPoEOwn4bT309r1nqQGM+adbF/JmEfDWc6XpuG3X5jnwx6IU0Y+ONgazs58k8n2e6JwMVV9f6qOr2qzgcePMW6L2jrfnxf3VsBjxhvg6r6fVVd0Lv1bXd/4BTgTOCAqhoc+zsN2CnJ9n3L9qYJAGdMoc2nAUsGlu0NnD7W/AxJku6qjXpEo6puTPIF4ENJrgV+D7wHWMfEpxQuAB6U5AXA/wD70kzonErdNyX5XFv3NTQjC+9v6x5aO5JxCnAJzdUf9+2NGlRV7+qRE4FfAEcneROwNc3IyJFVtapvXzvSnM7ZGtgiyc7tfs5ui3wa+KckH6WZ/LoHzTyTF06lzZIkDWujDhqtN9AcQL9LM8ficOCBwK3jbVBVxyX5Z+ATNFdUHE9zVcrSKdb9Jpr5Hd+muWLlIzRXjEzF3jTzRB7OHSeSpm3v2iT7Ap8EfkTzeRnHcMfLUpex/sjMWQP7uSjJXwH/DPx/NJNmX+NnaEiSurLRB412suaLeo+T3INm7sVn2vWLxtnuzcCbBxYf0bd+KQPBo6qOGChzE/CS9tbzoYFtxqy/b/1RwFETlWnLXUJzme1EZSasqy3zfZrLfyVJ6txGHzSSPJbm0zV/SjOa8O521XEz1ihJkgTMgaDReiOwA3AbzeTIPavq6pltkiRJ2uiDRlWdRfPpmJIkaZbZ2C9vlSRJs5hBQ5IkdcagIUmSOmPQkCRJnTFoSJKkzhg0JElSZwwakiSpMwYNSZLUGYOGJEnqjEFDkiR1ZqP/CHJJc8vKw/ad6SbMOqOjoyxbtoxzl+7DyMjITDdnVrKPhtPrpw3JEQ1JktQZg4YkSeqMQUOSJHXGoCFJkjpj0JAkSZ0xaEiSpM4YNCRJUmcMGpIkqTMGDUmS1Bk/GVTSRmPRW4+f6SbMiAXzi8N3g52WLmfN2sx0c2Yl+2g4vX7akBzRkCRJnTFoSJKkzhg0JElSZwwakiSpMwYNSZLUGYOGJEnqjEFDkiR1xqAhSZI6Y9CQJEmdMWhIkqTOGDQkSVJnDBqSJKkzBg1JktQZg4YkSerMJhk0kixKUkl2nsI2pyQ5ortWSZI092ySQQP4LbA9cO5MNyTJ9kmOSfKrJOvGCjNJRpK8O8mFSW5N8rMkzxgos0WSI5JcnOSWJKcmedxAmXsm+USSS9sy5yV5ZcdPUZK0CdvkgkaSzatqbVVdUVV/mOn2AAuAq4D3Az8bp8z7gH8AXg3sCHwa+GaSx/aV+SywBHgJ8GjgROCkJPfvK/PPwDOAFwOPah//S5L/M23PRpKkPht90GhPaXyivV2f5Jok70uSdv3KJO9MclSSG4Ajxzp1kmTHJMuS3JTkyiRHJ7nPQHWbjVdPu4/Nkxye5HdJbk7ykyR7TdT+qlpZVa+tqi8CN4xT7CXAB6pqWVX9pqo+BSwHDmrrvTvwPODNVfWDqrqgqpYCFwH9IxZ7AF+oqlPaej9DE252nbCTJUm6kzab6QZMk5cBnwN2pzlofga4GDiyXf8m4BCakYE7SLI98P22/BuAuwMfBI4FnjqFev4NWAS8ALgMeA5wQpJHV9X5d+H5LQBuHVh2C/Ck9v5mwPxJygD8N/CsJJ9v27cX8AjgteNVnGRBW3/PFgCjo6OMjo5O6UlsSnp9Yx9NbKr9tGB+ddmcWWvBvFrvp+7IPhpOr3+m47Vp2H2kauP+pSQ5BdgW+LNqn0ySw4BnVdWOSVYCZ1XVc/q2WUTzbv+xVXV2kvcCu1fVPn1lHkAzl2OHqvr1EPU8DDgfeEBVXda3n5OAn1bV24d8LmdX1esGlh8DPAZ4NnAh8DTgOGB+VS1oy5wK3AbsB1wJvBD4InB+Ve3QltmcJhS9FPgDsA74+6o6eoI2LQUOHlx+zDHHsHDhwsmekiRpjlq9ejX77bcfwJZVtWq8cnNlROPHtX5iOg04KMn89vHpk2y/C7A4yU1jrHsY8Osh6vkLIMCv+86mQDMacA3AwP7/var+cZJ29byWJiD8EiiasPFvwAF9ZV4CfB74HbAWOBM4pm1Xz2uAxwPPohmJeTLwySSXV9VJ49R9KPDRvsdbAJcuXryYbbbZZsjmb3pGR0dZsWIFS5YsYWRkZKabM2tNtZ92Wrp8A7Rq9lkwrzhk13W86/R5rFmXyTfYBNlHw+n103S8Nq1aNW62WM9cCRqTuXmS9fOA7wBvGWPd5UPWMY/mAL9L+7NfL2Ds3LdsuN8QUFVXAc9OcjdgG5rTHofRjMr0ylwIPCXJPYB7VdXlSb7aK9PO4/gA8JyqOr7d7OftPJU3AmMGjapaA6zpPe6FqJGREQ+gQ7CfhjNsP61Zu2kfQNasyybfB5Oxj4YzHa9Nw24/V4LG48d4fH5VrR0YXRjPmTSTKVdOciXKRPWcRTNPYtuq+uFYG1fVBcM0ZjxVdSvwuyQjbXuPHaPMzcDNSbYC9gHe3K4aaW/rBjZZyxyYFCxJmp3mygHmgUk+mmSHJC+kuQz0Y1PY/l+BrYEvJ9ktyUOT7J3k832nXyasp6p+DXwJ+GKS5yZ5SJLHJXlLkr+aqPIkO7cjC/cE7ts+3rFv/e7tPh+aZE/gBJrf3eF9ZfZJ8oy23iXAycCvaE6x0J4/+z7woSR7teX2p5mv8c0p9JUkSUObKyMaX6S5UuSnNO/Q/4XmipChVNVlSZ5Ic6XJcpp5FRfTHND7RwAmq+cA4J3AR4D708zNOA1YNkkTzuq7vwvNhM6Laa5gAbgbzRUzD6U5DbMMeElVXd+33ZY08ykeAFwLfB14R1X1Twt+QVvmSzTB6mLgHTSfyyFJ0rSbK0FjtL1S4w6fcllVi8ZYtpJm4mb/svOB545XQVXt1fdwzE/TbA/qBzPGVRoTqaoJz+9U1fdpPqhrojLHMsaplIEyV7D+BFJJkjo1V06dSJKkWcigIUmSOrPRnzoZOKUhSZJmEUc0JElSZwwakiSpMwYNSZLUGYOGJEnqjEFDkiR1xqAhSZI6Y9CQJEmdMWhIkqTOGDQkSVJnDBqSJKkzG/1HkEvadKw8bN+ZbsKMGB0dZdmyZZy7dB9GRkZmujmzkn00nF4/bUiOaEiSpM4YNCRJUmcMGpIkqTMGDUmS1BmDhiRJ6oxBQ5IkdcagIUmSOmPQkCRJnTFoSJKkzvjJoJLUWvTW42e6CWNaML84fDfYaely1qzNTDdnVrKPhtPrpw3JEQ1JktQZg4YkSeqMQUOSJHXGoCFJkjpj0JAkSZ0xaEiSpM4YNCRJUmcMGpIkqTMGDUmS1BmDhiRJ6oxBQ5IkdcagIUmSOmPQkCRJnTFoSJKkzhg0xpBkZZLXzXQ7JEna2E0paCQ5JckR092IJJXk2dO939kuyUiSDyY5J8nNSS5L8sUk9xsotyDJvyS5ui337SQPGCizVZKjk9zQ3o5Ocu8x6tw/yc+T3JrkiiSf6PhpSpI2YY5obCBJNh9j8ULgL4BD2p/PBR4BfHug3BHAc4AXAE8C7gl8N8n8vjLHADsDz2hvOwNHD7ThDcD7gcOAPwOeBiy/009KkqRJDB00khwFPAV4bTsCUUkWJdkxybIkNyW5sn0nfZ++7U5J8vEkhye5tn0XvbRv/cr27jfbfa7s1ZfkWwNtOCLJKX2P/6YdDbglyTVJTkpyj3bdXkl+2o4AXJ/kR0ke3Lfts5Kc3r6zvzrJNwae8sIkn09yY5JLkrxioC33T/LVJNe1dR+XZFF/fyX5VpK3JbkM+PVgn1bVDVW1pKqOrapfVdWPgVcDuyR5ULufLYEDgYOq6qSqOgt4MfBo4OltmUfRhIu/r6rTquo04OXAM5Ps0JbZCngf8NKqOqaqLqyq/62q7wy2S5Kk6bLZFMq+lubd9rnAu9tl84HvA0cCbwDuDnwQOBZ4at+2LwM+CuwO7AEcleRHVbUCeBzwe+AA4ARg7TCNSbI98GXgzcA3gS2APZtV2Qz4VtuuFwKbA7sB1W67L/ANmnf3L2nX7ztQxUHAu4APAH8DfCrJD6rql0kWAicDPwSeDPwBeCdwQpI/r6rb2n08DVgFLAEyzPMCtmzbeX37eBdgBDixV6CqLktyLvAEmhGJPYAbquonfWV+nOSGtsyv2jbMA+6f5Ly2v06lCTC/Ha8xSRYAC/oWbQEwOjrK6OjokE9p09PrG/toYrOtnxbMr5luwpgWzKv1fuqO7KPh9PpnOv7nht3H0EGjqm5IchuwuqquAEjyXuDMqnp7r1ySvwN+m+QRVdV7F//zqnpPe//8JP9EcxBeUVVXJQG4vrffIW3ftv8bVXVxu+yctg1b0xywv1tVF7brzuvb9h3AV6rq4L5lPxvY/7Kq+mS7vw8Crwf2An5JcwpjHc0IQi+8HEATDvbi9lBwc1vmNoaQ5G40pzWOqapV7eLtgNuq6rqB4le263plfj/GLn/fV+ahNEHj7TSh8QaaEY4VA+Fo0NuAgwcXnnzyySxcuHCYp7VJW7FixUw3YaMwW/rp8N1mugUTO2TXdTPdhFnPPhrOdPzPrV69eqhyUxnRGMsuwOIkN42x7mHcfrrg5wPrLge2vYt1/wz4HnBOkuU0B/evVdV1VXVte6pneZIVwEnAsVV1ebvtzjSjHRP5Y5urqpJc0dfmXYCHAze2IannbjTPu+ec3gE8yYuA/9e37i+r6oe9B0lGgK/QhIFXTdI2aEZI+qP7WDG+v8w8mpGR11TViW2dLwSuABYz/lyNQ2lGo3q2AC5dvHgx22yzzRDN3DSNjo6yYsUKlixZwsjIyEw3Z9aabf2009LZOWVpwbzikF3X8a7T57Fm3bCDo5sW+2g4vX6ajv+5VatWTV6Iux405gHfAd4yxrrL++4Pjq8Uk88PWccdTzf8sVeqam2SJTSnBvammdvw/iS7V9VFVXVAko/TzF14PvC+JEvaeRC3TFL3ZG2eB5wBvGiM7a7qu39z3/1vAz/pe/y73p02ZBwLPAR4at9oBjRBYPMkWw2MamxLc+qjV+ZPxmjLfWlGPuD238cv/viEmtGkq4EHjbFtr8waYE1fWwEYGRmZFQeG2c5+Gs5s6ac1a2f3AWrNusz6Ns40+2g40/E/N+z2U73q5DaaeRk9Z9JcvbCyqi4YuN089i7GNDqwX2gO2NsPLNu5/0E1ftSeAnls277n9K0/q6oOraon0Mwt2a9d9XOaUzd31pnAnwK/H+N53zDWBlV140C5W2C9kPGnwNOr6pqBTc+g6Z8lvQXt/JSduD1onAZsmWS3vjK705w+6pX5Uftzh74yWwP3AXqnniRJmlZTDRorgd3bq03uA/wrsDXw5SS7JXlokr3bqzUGg8Nk+31aku3aqyMA/gvYNclLk/xpkvfQHFyB5kCa5O1Jdm2v0HguzTv485I8JMmhSfZI8uAke9NMZO3N03gP8MIk70nyqCSPTvLmKbT3S8DVwHFJ9mzre0qSj2Xg8y0m0k5a/RqwK83oyPy2D7ZLezlsG1w+B3wkydOSPBb4d5r5KCe1Zc6jmUh7ZJLHJ3k8zamh71bVr9oyvwaOAz6W5AlJdgK+QDPn5OQpPHdJkoY21aDxYZqrQn5BM+KwOfBEmtGI5TSjBh+jmWg4lRk5B9G8Y/8tcBZAVS2n+XyJw4H/oZkb8MW+bVbRXPGxjGYuyPtorqD4T2A18Ejg6+26zwCfoJ0jUVWnAH8LPAs4mybU7D5sY6tqdVv3JTRXr5wHfJ7mqpvhTlo1HtC24QFtOy7vuz2hr9zraa6iOZZmZGI18NdV1X+FzotowseJ7e3nNFfU9Hspzemb42muFhoFnlFVs2PKvyRpzpnSHI32XfEeY6x67gTb7DXGsmcPPP4OzVyPwXIHM8YVD+2682jmX4y17kr6TqGMU+YbNCFhrHWLxli288DjK2gu2x1v//tPVH9bZiVDXPZaVbfSzEF59QRlrqX5fI2J9rOK5jM5DpysTkmSpoOfDCpJkjpj0JAkSZ0xaEiSpM4YNCRJUmcMGpIkqTMGDUmS1BmDhiRJ6oxBQ5IkdcagIUmSOmPQkCRJnTFoSJKkzhg0JElSZwwakiSpMwYNSZLUmSl9TbwkzWUrD9t3ppswptHRUZYtW8a5S/dhZGRkppszK9lHw+n104bkiIYkSeqMQUOSJHXGoCFJkjpj0JAkSZ0xaEiSpM4YNCRJUmcMGpIkqTMGDUmS1BmDhiRJ6oyfDCppo7forcfPdBM6tWB+cfhusNPS5axZm5luzqxkHw2n108bkiMakiSpMwYNSZLUGYOGJEnqjEFDkiR1xqAhSZI6Y9CQJEmdMWhIkqTOGDQkSVJnDBqSJKkzBg1JktQZg4YkSeqMQUOSJHXGoCFJkjpj0BhDkpVJXjfT7ZAkaWM3paCR5JQkR0x3I5JUkmdP935nuyQjST6Y5JwkNye5LMkXk9xvoNyCJP+S5Oq23LeTPGCgzFZJjk5yQ3s7Osm9x6l3mySXtv0+ZhlJkqaDIxobSJLNx1i8EPgL4JD253OBRwDfHih3BPAc4AXAk4B7At9NMr+vzDHAzsAz2tvOwNHjNOdzwM+n/iwkSZqaoYNGkqOApwCvbd8JV5JFSXZMsizJTUmubN9J36dvu1OSfDzJ4UmuTXJFkqV961e2d7/Z7nNlr74k3xpowxFJTul7/DftaMAtSa5JclKSe7Tr9kry03YE4PokP0ry4L5tn5Xk9CS3tiMF3xh4yguTfD7JjUkuSfKKgbbcP8lXk1zX1n1ckkX9/ZXkW0neluQy4NeDfVpVN1TVkqo6tqp+VVU/Bl4N7JLkQe1+tgQOBA6qqpOq6izgxcCjgae3ZR5FEy7+vqpOq6rTgJcDz0yyw0C7XwncG/jwYHskSZpuUxnReC1wGnAksH17GwW+D5wN7EpzsPsT4NiBbV8G3AzsDrwZeHeSJe26x7U/D2j3+TiGkGR74MvA54FHAXsB32hWZTPgW23b/hzYA/gMUO22+7ZljwceCzwNOH2gioPaZY8FPgl8Kskj2+0XAicDNwFPphlluAk4YWDk4mlt25YAzxzmeQFbtu28vn28CzACnNgrUFWXAecCT2gX7QHcUFU/6SvzY+CGvjIk2RF4N/BSYN2Q7ZEk6U7bbNiCVXVDktuA1VV1BUCS9wJnVtXbe+WS/B3w2ySPqKreu/ifV9V72vvnJ/knmoPwiqq6KgnA9b39Dmn7tv3fqKqL22XntG3YmuaA/d2qurBdd17ftu8AvlJVB/ct+9nA/pdV1Sfb/X0QeD1NmPklzSmMdTQjCL3wcgBNONiL20PBzW2Z24Z5QknuBhwGHFNVq9rF2wG3VdV1A8WvbNf1yvx+jF3+vlcmyQKaYPamqrokyUOHbNMCYEHfoi0ARkdHGR0dHWYXm6Re39hHE5uuflowv6ajObPWgnm13k/dkX00nF7/TMdr07D7GDpojGMXYHGSm8ZY9zBuP10wOB/gcmDbu1j3z4DvAeckWU5zcP9aVV1XVde2p3qWJ1kBnAQcW1WXt9vuTDMyM5E/trmqKskVfW3eBXg4cGMbknruRvO8e87phYwkLwL+X9+6v6yqH/YeJBkBvkIzyvSqSdoGENoRml4zJylzKHBeVf37EPvu9zbg4MGFJ598MgsXLpzirjY9K1asmOkmbBTuaj8dvts0NWSWO2RXByInYx8NZzpem1avXj1UubsaNOYB3wHeMsa6y/vuD8aeYvLTNutoDpT9Rv64g6q17emXJwB708xteH+S3avqoqo6IMnHaU7nPB94X5Il7SmFWyape7I2zwPOAF40xnZX9d2/ue/+t4Gf9D3+Xe9OGzKOBR4CPLVvNAPgCmDzJFsNjGpsC5zaV+ZPxmjLfWlGPgCeCjw6yd/0qm1/Xp3k/QOjO/0OBT7a93gL4NLFixezzTbbjLOJRkdHWbFiBUuWLGFkZGTyDTZR09VPOy1dPo2tmn0WzCsO2XUd7zp9HmvWDb4sCuyjYfX6aTpem1atWjV5IaYeNG4D+q90OBN4HrCyqv4wxX31Gx3YLzQH7J0Glu1MXwBoT1v8CPhRexrnYpqrMz7arj8LOAs4NMlpwH7Aj2lGK54G/NudbO+ZNOHl9wOhYFxVdSNw4+DyvpDxp8DiqrpmoMgZNM95SVuuNz9lJ5r5LtDMndkyyW5V9dO2zO40p496YeR5wN379vs4mvktewIXMo6qWgOs6WsvACMjIx5Ah2A/Deeu9tOatZvGgWXNumwyz/XOso+GMx2vTcNuP9XLW1cCu7dXm9wH+Fdga+DLSXZL8tAke7dXawwGh8n2+7Qk2yXZql32X8CuSV6a5E+TvIe+4JFk9yRvT7Jre4XGc2newZ+X5CFJDk2yR5IHJ9mb5rLR3jyN9wAvTPKeJI9K8ugkvYP2ML4EXA0cl2TPtr6nJPlYBj7fYiLtpNWv0UykfREwv+2D7XqTSqvqBprLUT+S5GlJHgv8O818lJPaMucBJwBHJnl8ksfTnBr6blX9qi1zYVWd27sBF7XNOK+qxprfIUnSXTbVoPFhYC3wC5oRh82BJ9KMRiynuRLiYzRXO0zlRNlBNO/Yf0szAkFVLaf5fInDgf+hGbL/Yt82q2iu+FhGMxfkfTSXgP4nsBp4JPD1dt1ngE/QzpGoqlOAvwWeRXPFzH/RXBEzlKpa3dZ9Cc3VK+fRjA7cvW3XsB7QtuEBbTsu77s9oa/c62muojmWZgRnNfDXVbW2r8yLaMLHie3t58BLptAWSZKm3ZROnbRXkewxxqrnTrDNXmMse/bA4+/QzPUYLHcwY0xEbNedRzP/Yqx1V9KcQhlXVX2DJiSMtW7RGMt2Hnh8Bc1lu+Ptf/+J6m/LrOSO81DGKncrzRyUV09Q5lqaz9cYShu2HF+UJHXKTwaVJEmdMWhIkqTOGDQkSVJnDBqSJKkzBg1JktQZg4YkSeqMQUOSJHXGoCFJkjpj0JAkSZ0xaEiSpM4YNCRJUmcMGpIkqTMGDUmS1BmDhiRJ6syUviZekmajlYftO9NN6NTo6CjLli3j3KX7MDIyMtPNmZXso+H0+mlDckRDkiR1xqAhSZI6Y9CQJEmdMWhIkqTOGDQkSVJnDBqSJKkzBg1JktQZg4YkSeqMQUOSJHXGTwbVRm/RW4+f6SbcwYL5xeG7wU5Ll7NmbWa6ObOW/TQc+2ly9tFwev20ITmiIUmSOmPQkCRJnTFoSJKkzhg0JElSZwwakiSpMwYNSZLUGYOGJEnqjEFDkiR1xqAhSZI6Y9CQJEmdMWhIkqTOGDQkSVJnDBqSJKkzBg1JktQZg0afJKckOWKm2yFJ0lxh0Fjfc4F3bcgKk+yV5Lgklye5OcnZSV40RrmnJDkjya1JfpPkH8co87wkv0iypv35nAnqfVuSMlhJkrpk0OhTVddW1Y1d7T/JyBiLnwD8HHge8OfA54EvJvnrvu0eAiwDfgg8FvgA8PEkz+srswfwVeBo4DHtz2OT7D5GOx4HvKKtV5Kkzhg0+vROnSR5d5Jzxlh/RpL39j0+IMl57SjDL5O8qm/donbE4P+2+70VePHgPqvqA1X1rqo6taourKqPAycA/aMR/whcUlWvq6rzquqzNIHkjX1lXgesqKpDq+qXVXUo8L12ef9zuCfwJeDlwHVT7SNJkqZis5luwCz1eeDgJI+rqv8BSPLnNKMJf9s+fjnwHuCfgLPadUcmubmqvtC3rw8CBwEHAGuGrH9L4Ly+x3sAJw6UWQ4cmGSkqkbbMv88RpnXDSz7V+D4qjopyTsna0iSBcCCvkVbAIyOjjI6OjrZ5hvEgvk10024gwXzar2fGpv9NBz7aXL20XB6/TMdr9/D7sOgMYaqujTJcppw8D/t4gOA71fVb9rH7wIOqqpvtI8vSrIj8A9Af9A4oq/MpJL8DfC4dj892wFXDhS9kub3dx/g8gnKbNe37xcAf9Huf1hvAw4eXHjyySezcOHCKeymO4fvNtMtGN8hu66b6SZsFOyn4dhPk7OPhrNixYq7vI/Vq1cPVc6gMb4jgc8neQOwFngRzcgESe4LPBD4XJIj+7bZDLhhYD+n9+4k+V/gwe3DH1bVX/YXTLIXcBTw8qr634H9DMb0jLF8rDLV7vuBwMeAvavqVoZ3KPDRvsdbAJcuXryYbbbZZgq76c5OS5fPdBPuYMG84pBd1/Gu0+exZl0m32ATZT8Nx36anH00nF4/LVmyhJGRsaYNDm/VqlVDlTNojO87NKc6ntP+XAB8vV3Xm9vycuAnA9utHXh8c9/9vwJ6v9lb+gsleUpb5xuq6osD+7iCvpGJ1rbAH4BrJinTG+XYpX18RvLHf8L5wJOT/BOwoKoG205VraHvlE9v25GRkbv8Rzpd1qydvS8qa9ZlVrdvtrCfhmM/Tc4+Gs50vIYPu71BYxxV9YckX+D2uRVfqarV7bork/wOeGhVfWkK+7x4rOXtSMZ3gbdU1WfGKHIa8NcDy/YGTm/nZ/TKLGH9eRp7A6e2978HPHpgH/8G/BL44FghQ5Kku8qgMbHPcvukzCcOrFtKc4npKuA/aUY8dgW2qqqPMqQ2ZBxPc1rj60l6oxK3VdW17f1PA/+U5KM0p3T2AA4EXti3q48BP0jyFuA44P8ATweeBNBetnvuQN03A9dU1XrLJUmaLl7eOoGqOp9mROBXVfWTgXWfBf4e2B84B/h+e/+iKVazP7CQZtLl5X23P04graqLaE677AWcTTMR9TVV9fW+MqcCL6AZgfl5u9/nD7ZbkqQNyRGNPlW1V//jNBMS/gT4f+OUPwY4Zpx1K7l9wuZEde5PEwomK/d9mitGJirzNeBrk+2rr/xew5aVJOnOMGiMI8m2wEuA+9PMZZAkSVNk0BjflcDVwCuqyk/QlCTpTjBojKOqvD5KkqS7yMmgkiSpMwYNSZLUGYOGJEnqjEFDkiR1xqAhSZI6Y9CQJEmdMWhIkqTOGDQkSVJnDBqSJKkzBg1JktQZP4JcG72Vh+070024g9HRUZYtW8a5S/dhZGRkppsza9lPw7GfJmcfDafXTxuSIxqSJKkzBg1JktQZg4YkSeqMQUOSJHXGoCFJkjpj0JAkSZ0xaEiSpM4YNCRJUmcMGpIkqTN+Mqikjcaitx4/002YEQvmF4fvBjstXc6atZnp5sxK9tFwev20ITmiIUmSOmPQkCRJnTFoSJKkzhg0JElSZwwakiSpMwYNSZLUGYOGJEnqjEFDkiR1xqAhSZI6Y9CQJEmdMWhIkqTOGDQkSVJnDBqSJKkzBg1JktSZTTJoJFmUpJLsPIVtTklyRHetkiRp7tkkgwbwW2B74NyZbkiS7ZMck+RXSdaNFWaSjCR5d5ILk9ya5GdJnjFQZoskRyS5OMktSU5N8riBMvdM8okkl7Zlzkvyyo6foiRpE7bJBY0km1fV2qq6oqr+MNPtARYAVwHvB342Tpn3Af8AvBrYEfg08M0kj+0r81lgCfAS4NHAicBJSe7fV+afgWcALwYe1T7+lyT/Z9qejSRJfTb6oNGe0vhEe7s+yTVJ3pck7fqVSd6Z5KgkNwBHjnXqJMmOSZYluSnJlUmOTnKfgeo2G6+edh+bJzk8ye+S3JzkJ0n2mqj9VbWyql5bVV8Ebhin2EuAD1TVsqr6TVV9ClgOHNTWe3fgecCbq+oHVXVBVS0FLgL6Ryz2AL5QVae09X6GJtzsOmEnS5J0J2020w2YJi8DPgfsTnPQ/AxwMXBku/5NwCE0IwN3kGR74Ptt+TcAdwc+CBwLPHUK9fwbsAh4AXAZ8BzghCSPrqrz78LzWwDcOrDsFuBJ7f3NgPmTlAH4b+BZST7ftm8v4BHAa8erOMmCtv6eLQBGR0cZHR2d0pPYlPT6xj6a2FT7acH86rI5s9aCebXeT92RfTScXv9Mx2vTsPtI1cb9S0lyCrAt8GfVPpkkhwHPqqodk6wEzqqq5/Rts4jm3f5jq+rsJO8Fdq+qffrKPIBmLscOVfXrIep5GHA+8ICquqxvPycBP62qtw/5XM6uqtcNLD8GeAzwbOBC4GnAccD8qlrQljkVuA3YD7gSeCHwReD8qtqhLbM5TSh6KfAHYB3w91V19ARtWgocPLj8mGOOYeHChZM9JUnSHLV69Wr2228/gC2ratV45ebKiMaPa/3EdBpwUJL57ePTJ9l+F2BxkpvGWPcw4NdD1PMXQIBf951NgWY04BqAgf3/e1X94yTt6nktTUD4JVA0YePfgAP6yrwE+DzwO2AtcCZwTNuuntcAjweeRTMS82Tgk0kur6qTxqn7UOCjfY+3AC5dvHgx22yzzZDN3/SMjo6yYsUKlixZwsjIyEw3Z9aaaj/ttHT5BmjV7LNgXnHIrut41+nzWLMuk2+wCbKPhtPrp+l4bVq1atxssZ65EjQmc/Mk6+cB3wHeMsa6y4esYx7NAX6X9me/XsDYuW/ZcL8hoKquAp6d5G7ANjSnPQ6jGZXplbkQeEqSewD3qqrLk3y1V6adx/EB4DlVdXy72c/beSpvBMYMGlW1BljTe9wLUSMjIx5Ah2A/DWfYflqzdtM+gKxZl02+DyZjHw1nOl6bht1+rgSNx4/x+PyqWjswujCeM2kmU66c5EqUieo5i2aexLZV9cOxNq6qC4ZpzHiq6lbgd0lG2vYeO0aZm4Gbk2wF7AO8uV010t7WDWyyljkwKViSNDvNlQPMA5N8NMkOSV5Icxnox6aw/b8CWwNfTrJbkocm2TvJ5/tOv0xYT1X9GvgS8MUkz03ykCSPS/KWJH81UeVJdm5HFu4J3Ld9vGPf+t3bfT40yZ7ACTS/u8P7yuyT5BltvUuAk4Ff0ZxioT1/9n3gQ0n2asvtTzNf45tT6CtJkoY2V0Y0vkhzpchPad6h/wvNFSFDqarLkjyR5kqT5TTzKi6mOaD3jwBMVs8BwDuBjwD3p5mbcRqwbJImnNV3fxeaCZ0X01zBAnA3mitmHkpzGmYZ8JKqur5vuy1p5lM8ALgW+Drwjqrqnxb8grbMl2iC1cXAO2g+l0OSpGk3V4LGaHulxh0+5bKqFo2xbCXNxM3+ZecDzx2vgqraq+/hmJ+m2R7UD2aMqzQmUlUTnt+pqu/TfFDXRGWOZYxTKQNlrmD9CaSSJHVqrpw6kSRJs5BBQ5IkdWajP3UycEpDkiTNIo5oSJKkzhg0JElSZwwakiSpMwYNSZLUGYOGJEnqjEFDkiR1xqAhSZI6Y9CQJEmdMWhIkqTOGDQkSVJnDBqSJKkzG/13nUjadKw8bN+ZbsKMGB0dZdmyZZy7dB9GRkZmujmzkn00nF4/bUiOaEiSpM4YNCRJUmcMGpIkqTMGDUmS1BmDhiRJ6oxBQ5IkdcagIUmSOmPQkCRJnTFoSJKkzhg0JElSZwwakiSpMwYNSZLUGYOGJEnqjEFDkiR1xqAhSZI6Y9CQJEmdMWhIkqTOGDQkSVJnDBqSJKkzBg1JktQZg4YkSeqMQUOSJHXGoCFJkjpj0JAkSZ0xaEiSpM4YNCRJUmcMGpIkqTMGDUmS1JnNZroB2jjdeOONjIyMzHQzZq3R0VFWr17NqlWr7KcJ2E/DsZ8mZx8NZzr7adWqVUOVS1XdpYq0aUmyCLhoptshSZo1HlBVvxtvpSMamqpr258PAG6cyYbMclsAl2I/TcZ+Go79NDn7aDjT3U9bAJdNVMCgoTvrxqoabtxsE5Skd9d+moD9NBz7aXL20XA66KdJ9+FkUEmS1BmDhiRJ6oxBQ1O1BnhP+1Pjs5+GYz8Nx36anH00nA3eT151IkmSOuOIhiRJ6oxBQ5IkdcagIUmSOmPQkCRJnTFoaExJ3pHk1CSrk1w/TpkHJflOkpuTXJ3k40k2Hyjz6CTfT3JLkt8leXf6PjFmLkryiCTHtX2yKsmPkiweKDNp320Kkuyb5Cft38fVSb4xsN5+ApIsSHJ2kkqy88C6TbqPkixK8rkkF7V/Rxcmec8Yr0WbdD/1JHlV21e3JjkjyZ5d1+kng2o8mwP/AZwGHDi4Msl84HjgKuBJwDbAF4AAr27L3AtYAZwMPA54BHAUcDPwka6fwAw6Hvg18FTgFuB1wHeTPKyqrhim7zYFSZ4HHAm8Hfgvmuf/6L719tPtDqf5mOfH9C+0jwB4JM2b5n8ALgB2ovm7ugfwRrCfepI8HzgCeBXwI5o++88kO1bVJZ1VXFXevI17A/YHrh9j+V8Ca4H79S17AXArcK/28SuB64EFfWXeCvyO9tLquXYD7gMUsGffsi3aZU8btu/m+o3mTc6lwIETlNnk+6mvH84Ddmz/jna2jybtszcBv7Gf7tAvPwE+NbDsPODQLuv11InurD2Ac6uq/8t0lgMLgF36yny/qtYMlLkfsGhDNHIGXEPzj/vSJPdIshnNu4YrgTPaMsP03Vz3F8D9gXVJzkpyeZL/TPJnfWU2+X5K8ic0785fAqweo8gm30fj2JLbvwAS7Cfa00S7ACcOrDoReEKXdRs0dGdtR3Pw/KOqug64rV03Zpm+x9sxB1XzFmEJ8Fiab0a8FXg98Iyqur4tNkzfzXUPbX8uBd4HPBO4Dvh+kq3bdZt0P7VzmY4CPl1Vp49TbJPuo7EkeRjN6ZBP9y22n5rR1vmM/ZrcaR8YNDYhSZa2k8kmuu06hV2O9bGyGVg+WCbjLJ/Vhu279uDwSeD3wJ7AbsBxNHM0tu/b5TB9t9GZwt9Y77Xn/VX19ao6AziA5vn/bd8u51w/TaGPXg3cCzh0kl3OuT6CO/d6leR+wAnAf1TVZwd2OSf76U4Y6zW50z5wMuim5RPAVyYps3LIfV0B7N6/IMlWwAi3J+YruGNS3rb9OZiqZ7th++6pNO/Ot6rbv4L5VUmWAC8DDmO4vttYDdtPW7T3f9FbWFVrkvwGeFC7aK7207B99E7g8cCarH+h1ulJvlRVL2Pu9hFM8fWqDRkn00xgf8VAubncT8O6mmaeylivyZ32gUFjE1JVV9P8sU2H04B3JNm+qi5vl+1N80U9Z/SV+UCSzavqtr4ylzF8oJkVhu27JAvbu+sGVq3j9nfxw/TdRmkK/XQGzfPdAfjvdtkIzdydi9tic7KfptBHr6EJGz33o5lX8HyaSX0wR/sIpvZ6leT+NCHjDOCAqhr8/5uz/TSsqrqt/b9bAnyzb9USmlHXTiv35u0ON5p3lTsD76aZa7Bze7tnu34+cA5wEs18hKcBvwX+pW8fW9K8kziG5pKz5wA3AAfN9PPrsN/uQ/Pi+HWaSxEfAXyI5lzwY4btu03hRnOZ3aU0L/g7AJ+leWe1lf00Zn8t4o5XnWzyfUQTwM4HvkczwXi73s1+ukNfPb99Lfo74FHAPwM3AQ/utN6ZfuLeZueNZhJajXHbq6/Mg4Dv0syGvwb4F/ouZW3LPBr4Ac2kyMuBg5mjl7b2Peddad55XgOsonk39ZcDZSbtu7l+oxm2/nAbLlbRfObKn9lP4/bXHYKGffTHS/DHeq0q+2nM/noVzYhybzTnyV3X6dfES5KkznjViSRJ6oxBQ5IkdcagIUmSOmPQkCRJnTFoSJKkzhg0JElSZwwakiSpMwYNSbNWkkXtl2ftPNNtuTOS7NW2/94z3ZbpkOTAJINfMz6d+1+Q5JIkm8RXt28qDBrSHJXkCUnWJjmhg31v1AFgAzoV2J7mo/fvsiQntr/Tx4+x7qi+bzUdTfKbJB9Oco92/aKs/82nNyb53yT/muRPh6h7AfBe4JDpeC5jqao1NJ8W+8Gu6tCGZ9CQ5q6/o/mY5ScledBkhTX9quq2qrqipuEjmNvf4R4032p64DjFTqAJNg+l+UK2V9EcuPs9vS3zGODtNN958bMkT5ukCc8DbqqqH96pJzC8LwF7JnlUx/VoAzFoSHNQ+y72/wKfovl+h/3HKPOsJKcnuTXJ1Um+0beukjx7oPz1SXr7uaj9eVZb9pS+cgckOa/d7y+TvGqSts5L8pYkFyRZ0w6dv2Og2EOTnJxkdZKfJdmjb/ttknw5yaXt+nOSvHCgjlOSfDzJ4UmuTXJFkqUDZR6Z5L/bdv8iydMH+yHJ/ZN8Ncl1Sa5JclySRRM8t/VOnSTZv+3Hfdo+uinJCUm2n6iPWgfQ/C4/BTy/N1IxYE0bbH5bVcfQHLSfPVDmmrbMb6rqOJrg8RPgc0nmT1D/C4BvDzy/o5J8K8nbk1zZPreDk2yW5ENtX1+a5O/6ttk8ySeSXN729cokb+utr6praEaC1vsdauNl0JDmpucDv6qqXwH/DhyQJL2VSfYFvgEcz+3fZnn6FPa/W/uz9+74ue1+Xw68H3gHzTvltwOHJHnZBPs6FHgLzZD8jsB+NF+01u/9NO/MdwZ+DXw5yWbturvRfDnUM2m+JfgzwNFJdh/Yx8uAm4HdgTcD706ypG33POBbNF+4tTvwirbOP0qykOaryG8Cngw8qb1/QpLNJ3h+gxYCbwRe0u7nQdxx1GE97e/uAODfq+qXbR/83yHquoXmy+vGVc1Xqn8MeDAw0dyIPRn7b+SpNN+g+mTgDcBSmkB0HU1ffhr4dJIHtuVfAzyrbf8OwItpvuSr30/b+jQXzPQ3yXnz5m36b8CPgNe29zcDrgKe3rf+VJqD1njbF/DsgWXXA/u39xcx9jeJXgK8cGDZO4FTx6lnC5pv9v37cdb36jmwb9mO7bJHTtD+44EP9z0+BfjhQJmfAoe1958BjLL+V4s/vb8faE5F/ZK+bx8GNqcJJ3uP04692n3cu328f/v4YX1lXgVcMcnvcwnwe2Cz9vHrgP8eKHMU8K2+x7sBVwNfneh31q57ZLvu/45T/73b9XuOUedKYF7fsl8CP+h7PJ8mkL2gffxxmq90H/dbnGnCyEUz/X/kbXpujmhIc0ySHWgOMl8BqKo/AF+lOVD27EzzYj+d9d4XeCDNEPxNvRtN0HjYOJs9ClgwRFt+3nf/8vbntm2985O8I8nP29MZNwF704wUjLeP3n62be/vAPy2qq7oW//TgfK7AA8Hbux7btfSjKiM9/zGsrqqLhynHeM5kCYw/KF9/GVg9/Z33e+ZbdtuBU4DfgC8eog29Ua7xptLcvf2561jrPvfakZFeq4Ezuk9qKq1NF/L3nuOR9H8/f2qPZ219xj7vIVm5EdzwGaTF5G0kTmQ5n/7d/1nS4DRJFtV1XU0L+QTKW4/+PRMOATP7adiX05zzr/f2nG2mawdPaN993sHw159BwGvp3mXfw7N6ZEjaEYbxttHbz+9fYTxD7I982hO0bxojHVXTbLtZO0Y7Os/SrI1zTyLkSSv7Fs1nyY8vqVv2cnAK9s6LquqwbrG05t4edE4669p27nVGOvGej7j9nVVnZnkIcBf0owaHZvkpKr6m77yWzO1PtUs5oiGNIe08xZeSnPw3bnv9hjgYm4/SP6cZl7GeK6imXvR2++fsv47zNvan3+cPFhVVwK/Ax5aVRcM3MY7gJ1PEzYmu+JhInsCx1XVv1fVz4DfAJNerjngl8CDkvxJ37LHDZQ5s93v78d4ftNy+eo4XgRcSvM73Lnv9jrgZX1zVQBubttz8bAho52f8hqakHHWWGWq6jbgFzSnre6yqlpVVV+tqpfTzCd6XhuoenYary3a+DiiIc0tz6R51/m5wYNfkq/RjHZ8AngP8L0kF9KcYtkM+MuqOrwt/l/APyX5Mc0bkg+y/rvU39MEhGckuRS4ta1vKfDxJKuA/6Q5LbIrsFVVfXSwsVV1a5IPAocnuY1mbsl9gT+rqs8N+ZwvoDlQPYFmAuIbgO2A84bcHmAFcCHwhSRvppk70psM2hvp+BLwJuC4JO+mOfg/iGYi7Ieq6tIp1DcVBwJfq6pz+xcmuZjm97IvcNwU9rdNku1oguNONIFlN2Df9jTHeJbTTIA9Ygp13UGS19OcLjobWAf8LXAFzRygnj2Bd92VejR7OKIhzS0HAieN8w7768DOSf6iqk6heYF/Fs0L/n/RXCHQcxDwW5pz/MfQXBWxureynSvwGuAfgMtoD3RV9Vng72kmPZ4DfL+9P96IBjRXm3yE5sOgzqOZTzLZnIXB7c+kORCeQnPQ+tYUtu/NI3g2cE/gf4DPAu9rV9/alllNc2XFJTRX7JwHfJ5m/sKqqdQ3rDSfkPkYmt/dYJtvBE5k/M/UGM9JNAf6c4DDaJ7Hn1fVyZNsdyTwV0m2nGJ9g26iOd1zOk1fLwL+qjfPo710eUvga3exHs0SqbrLnyMjSXNOkicC/w08fGDy5iYrybHAWVV1aId1/Edbxwe6qkMblqdOJAlI8hyad9vn01xd8jHgR4aM9byJZhSsE2k+5vxnwD93VYc2PEc0JAlI8lKaeQEPpPn8iZOAg6r5pEpJd5JBQ5IkdcbJoJIkqTMGDUmS1BmDhiRJ6oxBQ5IkdcagIUmSOmPQkCRJnTFoSJKkzhg0JElSZwwakiSpM/8/iK5rYY6tMzQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 500x600 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig = plt.figure(figsize=(5, 6), dpi=100)\n",
    "ax = fig.add_subplot()\n",
    "ax.set_xlabel('Acute change in APD (ms)')\n",
    "ax.set_xlim(-100, 5)\n",
    "ax.set_ylim(-1, len(models))\n",
    "ax.grid(True)\n",
    "\n",
    "positions = list(range(len(models)))\n",
    "ax.set_yticks(positions, names)\n",
    "ax.barh(positions, delta_apds)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0b0dae7",
   "metadata": {},
   "source": [
    "Note that we could have gone a lot further here.\n",
    "\n",
    "For example, we made the implicit assumption that all models used milliseconds, but we could have checked this (e.g. with [model.time()](https://myokit.readthedocs.io/en/stable/api_core/Model.html#myokit.Model.time) and [variable.unit()](https://myokit.readthedocs.io/en/stable/api_core/Variable.html?highlight=convert_unit#myokit.Variable.unit)) and converted if necessary (with [variable.convert_unit()](https://myokit.readthedocs.io/en/stable/api_core/Variable.html?highlight=convert_unit#myokit.Variable.convert_unit)),\n",
    "\n",
    "Another tacit assumption was that `gCaL` pointed to a variable with a literal value, so that it could be changed with [Simulation.set_constant](https://myokit.readthedocs.io/en/stable/api_simulations/Simulation.html#myokit.Simulation.set_constant).\n",
    "Alternatively, we could have forced the variable to become a constant (e.g. `gCaL.set_rhs(gCaL.rhs().eval())`) or even modified the model and created a new simulation (e.g. `gCaL.set_rhs(myokit.Multiply(myokit.Number(0.5), gCaL.rhs()))`).\n",
    "\n",
    "A project that incorporates similar ideas in an organised way, is the [Cardiac Electrophysiology Web Lab](https://chaste.cs.ox.ac.uk/WebLab)."
   ]
  }
 ],
 "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.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}