{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# INCAWrapper validation case small toy model\n", "This notebook serves as a validation case the show that the INCAWrapper produce similar results as a model run through the INCA GUI. Here we compare the MCMC samples from the INCAWrapper to MCMC samples obtain from the GUI. \n", "\n", "This notebook is not meant as a tutorial and therefore code description is a bit more sparse. For a proper tutorial see the other examples at https://incawrapper.readthedocs.io/en/latest/examples/index.html.\n", "\n", "## Introduction\n", "In this validation case, we will use a simple toy model has been used in the literature to study 13C-MFA [1], [2]. The model contains 5 reactions and 6 metabolites. To ensure that we can obtain a good fits to the data, we employ a simulated dataset for this validation. The simulated dataset holds one experiment in which the network is fed a labelled version of the metabolite A ([2-C13]-A).\n", "\n", "## Method for INCA GUI based flux estimation and MCMC sampling\n", "The model, experiments and data was manually entered into the INCA GUI. This model was saved to a file (`docs/examples/Literature data/simple model/simple_model_mcmc_validation_gui_model.mat`). We then ran first the estimate, then clicked `Update model`, and the second `Monte carlo sampling` after sampling has converged we used `save fluxmap as` to save the results a file (`docs/examples/Literature data/simple model/simple_model_mcmc_validation_gui_fluxmap.mat`).\n", "\n", "### Extracting MCMC samples from INCA GUI\n", "As noted in the MCMC tutorial for the INCAWrapper, it is currently not capable of loading MCMC results if the MCMC has ben run through the GUI. Therefore we have exported the GUI samples through the Matlab command line, using the the following workflow:\n", "- Open INCA \n", "- In INCA open the fluxmap containing the MCMC samples\n", "- In Matlab command line: run `csvwrite('simple_model_gui_mcmc_results.csv', f.par.vals)`\n", "\n", "## Method for INCAWrapper flux estimation and MCMC sampling\n", "The MCMC was run through the incawrapper using the the Python script `docs/examples/Literature data/simple model/incawrapper_mcmc.py`.\n", "\n", "## Note about randomness in INCA\n", "When INCA estimates the flux distribution it deploys an optimisation algorithm, which searches for a local optimum in the parameter space. To increase the probability that the found flux distribution is a global optimum INCA can be configured to restart the optimisation algorithm at different point in the parameter space. Unfortunately, the INCA manual do not describe any method to set the random seed for random restarts, thus the best we can do is to a large number of restarts in both the INCA GUI and the INCAWrapper to improve the probability that the two executions finds the same optimum. Therefore, we used 1000 restarts during flux estimation in both the INCAWrapper and the INCA GUI. Once an optimum was found for the two models, we performed MCMC sampling both through the GUI and the wrapper. In this notebook we compare these samples to explore if they could be draw from the same population, i.e. the differences in flux solution between the GUI and the Wrapper are likely due to randomness.\n", "\n", "## Setting up the environment\n", "First, we will setup the coding environment, load packages, set pah to files and read-in the data." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "import pathlib\n", "import incawrapper\n", "import ast\n", "PROJECT_DIR = pathlib.Path().cwd().parents[1].resolve()\n", "data_directory = PROJECT_DIR / pathlib.Path(\"docs/examples/Literature data/simple model\")" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# File for which MCMC is started from\n", "incawrapper_fluxmap_file = data_directory / 'simple_model_incawrapper_fluxmap_mcmc.mat'\n", "# gui_fluxmap_file = data_directory / \n", "\n", "# MCMC samples\n", "incawrapper_mcmc_samples_file = data_directory / \"simulated_data\" / \"simple_model_incawrapper_mcmc_results.csv\"\n", "gui_mcmc_samples_file = data_directory / \"simulated_data\" / \"simple_model_gui_mcmc_results.csv\"\n", "\n", "# Hotelling T2 test results\n", "hotelling_result_file = data_directory / 'simulated_data' / 'hotelling_test_results.txt'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will load the MCMC results and save the samples to a csv file." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
samplereactionfluxsource
00R199.868174INCAwrapper
11R199.531294INCAwrapper
22R199.752950INCAwrapper
33R199.924800INCAwrapper
44R199.498449INCAwrapper
\n", "
" ], "text/plain": [ " sample reaction flux source\n", "0 0 R1 99.868174 INCAwrapper\n", "1 1 R1 99.531294 INCAwrapper\n", "2 2 R1 99.752950 INCAwrapper\n", "3 3 R1 99.924800 INCAwrapper\n", "4 4 R1 99.498449 INCAwrapper" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "incawrapper_flux_samples = pd.read_csv(incawrapper_mcmc_samples_file, index_col=None)\n", "\n", "# prepare long format\n", "incawrapper_flux_samples['sample'] = incawrapper_flux_samples.index.values\n", "incawrapper_flux_samples_long = pd.melt(incawrapper_flux_samples, id_vars='sample', var_name='reaction', value_name='flux')\n", "incawrapper_flux_samples_long['source'] = 'INCAwrapper'\n", "incawrapper_flux_samples_long.head()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
samplereactionfluxsource
00R199.584GUI
11R199.342GUI
22R199.934GUI
33R199.811GUI
44R199.748GUI
\n", "
" ], "text/plain": [ " sample reaction flux source\n", "0 0 R1 99.584 GUI\n", "1 1 R1 99.342 GUI\n", "2 2 R1 99.934 GUI\n", "3 3 R1 99.811 GUI\n", "4 4 R1 99.748 GUI" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Load GUI flux samples\n", "gui_flux_samples = pd.read_csv(gui_mcmc_samples_file, index_col=None)\n", "gui_flux_samples['sample'] = gui_flux_samples.index.values\n", "gui_flux_samples.columns = incawrapper_flux_samples.columns\n", "\n", "gui_flux_samples_long = gui_flux_samples.melt(id_vars='sample', var_name='reaction', value_name='flux')\n", "gui_flux_samples_long['source'] = 'GUI'\n", "gui_flux_samples_long.head()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGwCAYAAABPSaTdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABNrklEQVR4nO3deVhUZf8/8PfMsG+DoA6SILggbrhXSF9Ss1zS1EwTl9zKXcMllyxNElFze9zFvVwws9R6SjM3DMnUJ3BJxRRDRUBDBhDZZs7vD3+cGAEZYIYZDu/XdXE9nHPuOfOZm3k6b+9zn3NkgiAIICIiIpIouakLICIiIjImhh0iIiKSNIYdIiIikjSGHSIiIpI0hh0iIiKSNIYdIiIikjSGHSIiIpI0C1MXYA60Wi0SExPh6OgImUxm6nKIiIhID4IgICMjA+7u7pDLSx6/YdgBkJiYCA8PD1OXQUREROVw584d1K1bt8TtDDsAHB0dATztLCcnJxNXQ0RERPpIT0+Hh4eHeBwvCcMOIJ66cnJyYtghIiKqYkqbgsIJykRERCRpDDtEREQkaQw7REREJGmcs0NERCah0WiQl5dn6jLIjFlaWkKhUFR4Pww7RERUqQRBQFJSEtLS0kxdClUBzs7OcHNzq9B98Bh2iIioUhUEndq1a8POzo43c6ViCYKArKwspKSkAADq1KlT7n0x7BARUaXRaDRi0HF1dTV1OWTmbG1tAQApKSmoXbt2uU9pcYIyERFVmoI5OnZ2diauhKqKgu9KReZ3MewQEVGl46kr0pchvis8jWVkubm5OHjwIBITE+Hu7o7evXvDysqq0ut48uQJNm7ciLt376Ju3boYM2aMODxIREQkZQw7RrRhwwbs3bsXgiCI69atW4d3330XY8eOrbQ65syZg6ioKHH5/PnzOHDgAAICAhAaGlppdRAREZkCT2MZyYYNGxAREaETdICns8sjIiKwYcOGSqnj2aBTWFRUFObMmVMpdRAREZkKw44R5ObmIiIiQly2tLTEyJEjYWlpKa6LiIhAbm6uUet48uSJTtBp2rQpli1bhqZNm4rroqKi8OTJE6PWQUREZEo8jWUEW7ZsEX/fvXs33N3dAQDvvfceEhMTMWjQILHduHHjjFbHypUrxd9/+OEHODg4AADatm2LzMxM9OzZU2w3e/Zso9VBRETlp9FoIJPJIJdzfKK8TNpzkZGR6NWrF9zd3SGTyXDgwAGd7YIgYO7cuahTpw5sbW3RpUsX3LhxQ6dNamoqBg8eDCcnJzg7O2PUqFHIzMysxE8BZGdnIy4uTvzZu3cvAEAulyMzM1NnW2ZmpviF3bt3r862uLg4ZGdnG6SGuLg4nDhxAgDwwgsvIDExUWdbwYRpADhx4oRBaiAiqi6++eYbtGjRAra2tnB1dUWXLl3w+PFjaLVahISEoG7durC2tkarVq1w+PBh8XUnT56ETCbTuXt0TEwMZDIZbt++DQDYvn07nJ2dcejQITRt2hTW1tZISEhATk4OZs6cCQ8PD1hbW6Nhw4Y6/7i+fPkyunfvDgcHB6hUKgwdOhQPHz6srC4xayYd2Xn8+DFatmyJkSNH4u233y6yfcmSJVi1ahV27NgBb29vfPrpp+jatSv+/PNP2NjYAAAGDx6M+/fv4+jRo8jLy8OIESMwevRo7N69u9I+R0JCAkaPHl1kvVarLXZ9Yc9uDw8Ph4+Pj8FqAIB79+49t47c3Fyd7eWtgYioOrh//z6CgoKwZMkS9O3bFxkZGTh9+jQEQcB//vMfLFu2DBs3bkTr1q2xdetWvPXWW7hy5QoaNWqk93tkZWVh8eLF2Lx5M1xdXVG7dm289957iI6OxqpVq9CyZUvEx8eLYSYtLQ2dO3fG+++/jxUrVuDJkyeYOXMmBgwYgOPHjxurK6oMmfDsDFoTkclk+O6779CnTx8AT0d13N3dMW3aNEyfPh0AoFaroVKpsH37dgwcOBBXr15F06ZNce7cObRr1w4AcPjwYfTo0QN3794VRy6elZOTg5ycHHE5PT0dHh4eUKvVcHJyKvY1ycnJUKvVJe4vKSlJXC58hdOHH36I3bt348GDB6hVqxYGDRqE//znP+L2ZycIu7m5wdrautj3ycvL05n387waAGD9+vVITU0FAAwZMgQ7d+4UtxVednFx0Tmd9rwalEolVCpVsduIiEqTnZ2N+Ph4eHt7i/9orWr+97//oW3btrh9+zbq1auns+2FF17AhAkT8PHHH4vrXnzxRbRv3x5r167FyZMn0alTJzx69AjOzs4Ano7stG7dGvHx8fDy8sL27dsxYsQIxMTEoGXLlgCAuLg4NG7cGEePHkWXLl2K1LRgwQKcPn0aR44cEdfdvXsXHh4euH79epX+B+zzvjPp6elQKpXPPX4DZjxnJz4+HklJSTp/VKVSiZdeegnR0dEYOHAgoqOj4ezsLAYdAOjSpQvkcjnOnj2Lvn37FrvvsLAwzJ8/X+9akpOTMXjIUOTnlX1CceFg8+DBA51lAJV26XfhoPPscmpqqt51WFla4Kuduxh4iKjaatmyJV577TW0aNECXbt2xRtvvIF33nkHCoUCiYmJCAgI0GkfEBCA2NjYMr2HlZUV/Pz8xOWYmBgoFAq8+uqrxbaPjY3FiRMnxLmZhd28ebNKhx1DMNvZTgWjFM8eVFUqlbgtKSkJtWvX1tluYWEBFxeXIqMchc2ePRtqtVr8uXPnTqn1aDWasn4ESdJotKYugYjIpBQKBY4ePYqffvoJTZs2xerVq9G4cWPEx8eX+tqCOZuFT6oU9xgEW1tbnTsHl3YT2MzMTPTq1QsxMTE6Pzdu3EBgYKC+H02yzHZkx5isra1LPE1THJVKhXXr1uoVioCnQ27Lli0rtd20adPKNIyr0WjK9BC0lJQUbNq0qdR2H3zwQZHQWBIPDw+O6hBRtSeTyRAQEICAgADMnTsX9erVw7Fjx+Du7o6oqCidEZioqCi8+OKLAIBatWoBeDrvp0aNGgCejtqUpkWLFtBqtTh16lSxp7HatGmD/fv3w8vLCxYW1fLQ/lxm2yNubm4Anp5CKvxY9+TkZLRq1UpsU/Do9wL5+flITU0VX28ovr6+8PX1LXZbdnY2EhISxOWFCxeKv3/88ccICwuDIAiQyWSYPXu2uH3//v0653UBwNPTs1znsZ+toXAdcrkcS5cuxcKFC/Hw4UPUrFkTH3/8MaZPnw6tVovNmzdj48aNFa6BiKg6OHv2LI4dO4Y33ngDtWvXxtmzZ/HgwQM0adIEH330EebNm4cGDRqgVatW2LZtG2JiYrBr1y4AQMOGDeHh4YHPPvsMoaGhiIuL0+sfx15eXhg2bBhGjhwpTlD++++/kZKSggEDBmDChAnYtGkTgoKCMGPGDLi4uOCvv/5CREQENm/eXO6nhUuF2YYdb29vuLm54dixY2K4SU9Px9mzZ8XJtP7+/khLS8OFCxfQtm1bAMDx48eh1Wrx0ksvVVqtz7sSqnDwEQRBZ/n27duVcjWWVqvF1KlTxeWHDx/qLAuCwKuxiIj05OTkhMjISKxcuRLp6emoV68eli1bhu7du6Nr165Qq9WYNm0aUlJS0LRpUxw6dEi8EsvS0hJ79uzBuHHj4Ofnh/bt22PBggXo379/qe+7fv16fPzxxxg/fjz++ecfeHp6iv9gLhhRmjlzJt544w3k5OSgXr166NatG+/PAxNfjZWZmYm//voLANC6dWssX74cnTp1gouLCzw9PbF48WIsWrRI59Lzixcv6lx63r17dyQnJ2PDhg3ipeft2rUr06Xn+s7mLsmzoyqTJ09GdnY2GjRogJkzZxZpv2jRIty6dQs2NjZYtWqVzjZDjuyMGTMGgiBALpcX+3iKsWPHQqvVQiaTcWSHiCqFFK7GospV5a/GOn/+PDp16iQuF4w2DBs2DNu3b8eMGTPw+PFjjB49GmlpaXjllVdw+PBhnQ+7a9cuTJw4Ea+99hrkcjn69etXJEAYm42Njc5ISOfOnfHjjz/i5s2bqFu3Luzs7MRtWVlZuHXrltjOUCMoz9YAAJs2bcL7778PrVYLOzs71K1bV9x29+5daLVasV3Dhg0NUgcREZG5MZv77JhSRUd2nvXkyRN0795dXG7fvj2GDh2Kr776CufOnRPX//TTT6XOsK+ojh07ir8rFAr0798f+/btg6bQ1WUnT540ag1ERAU4skNlZYiRHZ7IMwJbW1ud+yycO3cOkydP1gk6AQEBRg86gG6Q0Wg0iIiIYNAhIqJqhWHHSEJDQ4vcWKpAQEBApd1MEHgaaDZv3izes0Emk2Hz5s0MOkREVC2Y7dVYUhAaGoonT55g48aNuHv3LurWrYsxY8ZUyojOsxo2bCg+GJSIiKg6YdgxMltbWwQHB5u6DCIiomqLp7GIiIhI0hh2iIiISNJ4GouIiEwuOTkZarW60t5PqVTyOX/VCMMOERGZVHJyMoYMfQ95uTmV9p6WVtbY+dWXZQo8w4cPR1paGg4cOIDhw4djx44dCAsLw6xZs8Q2Bw4cQN++fXWeai4IAjZt2oQtW7bgypUrsLCwQMOGDTFkyBCMHj1a58azd+/eRf369eHj44PLly8b5sMSww4REZmWWq1GXm4OntR/FVobpdHfT56tBm6dglqtrtDojo2NDRYvXowxY8aITzAvztChQ/Htt9/ik08+wZo1a1CrVi3ExsZi5cqV8PLyQp8+fcS227dvx4ABAxAZGYmzZ88a7DmPgiBAo9GY9RPRc3NzYWVlZZR9c84OERGZBa2NElr7msb/MVCg6tKlC9zc3BAWFlZim6+//hq7du3Cnj178PHHH6N9+/bw8vJC7969cfz4cZ1HJgmCgG3btmHo0KEYNGgQtmzZIm67fPky5HI5Hjx4AABITU2FXC7HwIEDxTYLFizAK6+8AuDp/dVkMhl++ukntG3bFtbW1vj1119x8+ZN9O7dGyqVCg4ODmjfvj1++eUXnZq9vLzw+eefIygoCPb29njhhRewdu1anTYymQzr169H9+7dYWtri/r16+Obb77RaXPnzh0MGDAAzs7OcHFxQe/evXH79m1x+/Dhw9GnTx+EhobC3d0djRs31rPny45hh4iIqBwUCgUWLlyI1atX4+7du8W22bVrFxo3bozevXsX2SaTyaBU/hu8Tpw4gaysLHTp0gVDhgxBREQEHj9+DABo1qwZXF1dcerUKQDA6dOndZYB4NSpUzqPCAKAWbNmYdGiRbh69Sr8/PyQmZmJHj164NixY/jjjz/QrVs39OrVq8iDpL/44gu0bNkSf/zxB2bNmoUPP/wQR48e1Wnz6aefol+/foiNjcXgwYMxcOBAXL16FQCQl5eHrl27wtHREadPn0ZUVBQcHBzQrVs35Obmivs4duwYrl+/jqNHj+KHH34orcvLjWGHiIionPr27YtWrVph3rx5xW6/ceOG3iMWW7ZswcCBA6FQKNC8eXPUr18f+/btA/A0GAUGBop3vj958iRGjBiBnJwcXLt2DXl5eThz5gxeffVVnX2GhITg9ddfR4MGDeDi4oKWLVtizJgxaN68ORo1aoTPP/8cDRo0wKFDh3ReFxAQgFmzZsHHxweTJk3CO++8gxUrVui06d+/P95//334+Pjg888/R7t27bB69WoAwN69e6HVarF582a0aNECTZo0wbZt25CQkKBz9357e3ts3rwZzZo1Q7NmzfTqp/Jg2CEiIqqAxYsXY8eOHeKoRmH6Pms7LS0N3377LYYMGSKuGzJkiM6prFdffVUMCqdOnULnzp3FAHTu3Dnk5eUVeUxRu3btdJYzMzMxffp0NGnSBM7OznBwcMDVq1eLjOz4+/sXWX728z2vTWxsLP766y84OjrCwcEBDg4OcHFxQXZ2Nm7evCm+pkWLFkabp1OY+c5UIiIiqgICAwPRtWtXzJ49G8OHD9fZ5uPjg2vXrpW6j927dyM7O1tnQrIgCNBqtYiLi4OPjw86duyI4OBg3LhxA3/++SdeeeUVXLt2DSdPnsSjR4/Qrl07nSu7gKcjJ4VNnz4dR48exdKlS9GwYUPY2trinXfe0Tm1ZAiZmZlo27Ytdu3aVWRbrVq1SqzPWDiyQ0REVEGLFi3C999/j+joaJ31gwYNQlxcHA4ePFjkNYIgiPcW2rJlC6ZNm4aYmBjxJzY2Fv/3f/+HrVu3Ang6ClKjRg0sWLAArVq1goODAzp27IhTp07h5MmTRebrFCcqKgrDhw9H37590aJFC7i5uelMGi7w22+/FVlu0qSJ3m3atGmDGzduoHbt2mjYsKHOT+F5SpWFYYeIiKiCWrRogcGDB2PVqlU66wcMGIB3330XQUFBWLhwIc6fP4+///4bP/zwA7p06YITJ04gJiYG//vf//D++++jefPmOj9BQUHYsWMH8vPzxXk7u3btEoONn58fcnJycOzYsSLzdYrTqFEjfPvtt2KYGjRoELRabZF2UVFRWLJkCeLi4rB27Vrs27cPH374oU6bffv2YevWrYiLi8O8efPw+++/Y+LEiQCAwYMHo2bNmujduzdOnz6N+Ph4nDx5EpMnTy5xMrcxMewQEZFZkGerIX/80Pg/2ca5U3NISEiR4CCTybB7924sX74cBw4cwKuvvgo/Pz989tln6N27N7p27YotW7agadOm8PX1LbLPvn37IiUlBT/++COAp/N2NBqNGHbkcjkCAwMhk8mKzNcpzvLly1GjRg106NABvXr1QteuXdGmTZsi7aZNm4bz58+jdevWWLBgAZYvX46uXbvqtJk/fz4iIiLg5+eHL7/8Env27EHTpk0BAHZ2doiMjISnpyfefvttNGnSBKNGjUJ2djacnJz06k9Dkgn6zp6SsPT0dCiVSqjVapP8EYiIqovs7GzEx8fD29sbNjY2AKrOHZSrCy8vLwQHByM4OLjENjKZDN99953ODRGNpbjvTAF9j9+coExERCalUqmw86sv+WwsMhqGHSIiMjmVSsXwQUbDsENERESi4q7OelZVmwHDCcpEREQkaQw7REREJGkMO0RERCRpDDtEREQkaQw7REREJGkMO0RERCRpvPSciIhMLjk5mTcVJKNh2CEiIpNKTk7Ge0OHICc3r9Le09rKEl9+tbPMgScpKQlhYWH473//i7t370KpVKJhw4YYMmQIhg0bBjs7uxIfpTB8+HCkpaXhwIEDAICOHTuiVatWWLlypWE+FJWIYYeIiExKrVYjJzcPY5tmwN1eY/T3S3yswIY/HaFWq8sUdm7duoWAgAA4Oztj4cKFaNGiBaytrXHp0iWEh4fjhRdewFtvvWXEyqm8GHaIiMgsuNtr4OVo/LBTXuPHj4eFhQXOnz8Pe3t7cX39+vXRu3fvKndX4eqEE5SJiIhK8c8//+Dnn3/GhAkTdIJOYTKZrJKrIn0x7BAREZXir7/+giAIaNy4sc76mjVrwsHBAQ4ODpg5c6aJqqPSMOwQERGV0++//46YmBg0a9YMOTk5pi6HSsA5O0RERKVo2LAhZDIZrl+/rrO+fv36AABbW1txnaOjY7GX0aelpUGpVBq3UCoWR3aIiIhK4erqitdffx1r1qzB48ePn9u2cePGuHDhgs46jUaD2NhY+Pj4GLNMKgHDDhERkR7WrVuH/Px8tGvXDnv37sXVq1dx/fp17Ny5E9euXYNCoQAATJ06FZs3b8a6detw48YNxMTEYPTo0Xj06BHef/99E3+K6omnsYiIyCwkPlaY9fs0aNAAf/zxBxYuXIjZs2fj7t27sLa2RtOmTTF9+nSMHz8eABAUFARBELB8+XLMmjULdnZ2aNu2LSIjI3nXZhORCbwxANLT06FUKqFWq+Hk5GTqcoiIJCs7Oxvx8fHw9vaGjY0NgKp1B2WqfMV9Zwroe/zmyA4REZmUSqXCl1/t5LOxyGgYdoiIyORUKhXDBxkNJygTERGRpDHsEBERkaQx7BARUaXjtTGkL0N8Vxh2iIio0lhaWgIAsrKyTFwJVRUF35WC7055cIIyERFVGoVCAWdnZ6SkpAAA7Ozs+LRwKpYgCMjKykJKSgqcnZ3FmzaWB8MOERFVKjc3NwAQAw/R8zg7O4vfmfJi2CEiokolk8lQp04d1K5dG3l5lXcjQap6LC0tKzSiU4Bhh4iITEKhUBjkQEZUGk5QJiIiIklj2CEiIiJJY9ghIiIiSWPYISIiIklj2CEiIiJJY9ghIiIiSWPYISIiIklj2CEiIiJJY9ghIiIiSWPYISIiIklj2CEiIiJJY9ghIiIiSeODQImIiCTkyZMn2LhxI+7evYu6detizJgxsLW1rdQakpKSMG7cOGRmZsLBwQHr16+Hm5tbpdZQmEwQBMFk724m0tPToVQqoVar4eTkZOpyiIiIymXOnDmIiooqsj4gIAChoaGVUkO3bt2QnZ1dZL2NjQ0OHz5s0PfS9/jN01hEREQSUFLQAYCoqCjMmTPH6DWUFHQAIDs7G926dTN6DcVh2CEiIqrinjx5ohN0mjZtimXLlqFp06biuqioKDx58sRoNSQlJekEHS8vLyxcuBBeXl7iuuzsbCQlJRmthpLwNBZ4GouIiKq2xYsX46effgIA/PDDD3BwcBC3ZWZmomfPngCA7t27Y+bMmUapoVevXsjIyAAAHDp0SOd4mp6ejrfeegsA4OjoiO+//94g76nv8ZsTlImIiKqY7OxsJCQkiMu//PILAKBu3bpITEws0t7d3R2JiYn45Zdf0LdvX3G9p6cnbGxsDFJDQdBRKpVISkoqMoLj6OiIjIwMZGRkIC4uziA16IsjO+DIDhERVS1xcXEYPXp0hfcTHh4OHx+fKluDvsdvhh0w7BARkXm5du0a7ty5U+L2vLw8PHz4UFzev38/1Gp1qftVKpXo16+fuFyzZk1YWloW21aj0UChUOhdw+7du0ucnFyYjY0NBg0apFcNHh4e8PX1LXFfkgg7Go0Gn332GXbu3ImkpCS4u7tj+PDh+OSTTyCTyQAAgiBg3rx52LRpE9LS0hAQEID169ejUaNGer8Pww4REZmL5ORkBAUNglarMXUpJqeQy7F7zx6oVKpit0vi0vPFixdj/fr1WLNmDa5evYrFixdjyZIlWL16tdhmyZIlWLVqFTZs2ICzZ8/C3t4eXbt21StdEhERmSP5c0ZUqhOFwjAxxaxHdnr27AmVSoUtW7aI6/r16wdbW1vs3LkTgiDA3d0d06ZNw/Tp0wEAarUaKpUK27dvx8CBA/V6H47sEBGROUlOTn7uaamcnBydCcALFy6EPodzmUyGjz/+WFx2c3ODtbV1sW3z8vJKPL1UXA2LFy9Gfn5+qTVYWFjoXBH2vBqUSmWJozqARK7G6tChA8LDwxEXFwcfHx/Exsbi119/xfLlywEA8fHxSEpKQpcuXcTXKJVKvPTSS4iOji4x7OTk5CAnJ0dcTk9PN+4HISIiKgOVSvXcg3xcXFy57ogsCILO6yo6Qbk8NeTn5xusBn2ZddiZNWsW0tPT4evrC4VCAY1Gg9DQUAwePBgAxET57BdCpVI996ZFYWFhmD9/vvEKJyIiMiJPT0+Eh4eLy2PHjoVWq4VCocD69euLtB83bhw0Gg3kcjk2bNigsx9D1TB+/Hjk5+fDzs4OK1euLNI+ODgYWVlZsLCwwLp16wxSg77MOux8/fXX2LVrF3bv3o1mzZohJiYGwcHBcHd3x7Bhw8q939mzZ2Pq1Knicnp6Ojw8PAxRMhERkdHZ2NjojIZs3rwZI0eOhEajgYODA9zd3cVtiYmJ0Gg0Yrv69esbpYaPPvoIYWFhyMrKgouLC2rWrClue/jwIbKyssR2xh7JeZZZh52PPvoIs2bNEk9HtWjRAn///TfCwsIwbNgw8QmqycnJqFOnjvi65ORktGrVqsT9Wltbl3h+kIiIqKopHGAGDRoES0tLDBw4EBEREcjLyyu2naF16dIFYWFhAIB33nkHjo6OGDFiBLZt2ybecLCgXWUz66uxsrKyIJfrlqhQKKDVagEA3t7ecHNzw7Fjx8Tt6enpOHv2LPz9/Su1ViIiIlM6efKk+HteXh6++uornaBTeLsxKBQKhISEiMsZGRlYtWqVTtAJCQl57r17jMWsw06vXr0QGhqK//73v7h9+za+++47LF++XLzVtUwmQ3BwMBYsWIBDhw7h0qVLeO+99+Du7o4+ffqYtngiIqJKdvLkSWzdulUcKJDL5di6davRg06BwMBAhISEwNXVVWd9zZo1ERISgsDAwEqp41lmfel5RkYGPv30U3z33XdISUmBu7s7goKCMHfuXFhZWQH496aC4eHhSEtLwyuvvIJ169aV6XwgLz0nIiIyHI1Gg4sXLyI1NRUuLi7w8/MzyoiOJO6gXFkYdoiIiKoeSdxBmYiIiKiiGHaIiIhI0hh2iIiISNIYdoiIiEjSGHaIiIhI0hh2iIiISNIYdoiIiEjSGHaIiIhI0hh2iIiISNIYdoiIiEjSGHaIiIhI0hh2iIiISNIYdoiIiEjSGHaIiIhI0hh2iIiISNIYdoiIiEjSGHaIiIhI0hh2iIiISNIYdoiIiEjSGHaIiIhI0hh2iIiISNIYdoiIiEjSGHaIiIhI0hh2iIiISNIYdoiIiEjSGHaIiIhI0hh2iIiISNIYdoiIiEjSGHaIiIhI0hh2iIiISNIYdoiIiEjSGHaIiIhI0hh2iIiISNIYdoiIiEjSGHaIiIhI0hh2iIiISNIYdoiIiEjSGHaIiIhI0hh2iIiISNIYdoiIiEjSGHaIiIhI0hh2iIiISNIYdoiIiEjSGHaIiIhI0hh2iIiISNIYdoiIiEjSGHaIiIhI0hh2iIiISNIYdoiIiEjSGHaIiIhI0hh2iIiISNIYdoiIiEjSGHaIiIhI0ixMXQBRVfLkyRNs3LgRd+/eRd26dTFmzBjY2tqauiwiInoOmSAIgqmLMLX09HQolUqo1Wo4OTmZuhwqQW5uLg4ePIjExES4u7ujd+/esLKyqrT3nzNnDqKiooqsDwgIQGhoaKXVQURET+l7/GbYAcNOVbBhwwZEREQUWT9w4ECMHTvW6O9fEHQsLCygUCiQl5cHS0tLaDQa5OfnM/AQEZmAvsdvztkhs1dS0AGAiIgIbNiwwajv/+TJE3FEJz8/Hzk5OdBqtcjJyUF+fj4AICoqCk+ePDFqHUREVD4MO2TWcnNzdYKOpaUlRo4cCUtLS3FdREQEcnNzjVbDxo0bdZbd3Nwwb948uLm5PbcdERGZB05QJrNWOEDs3r0b7u7uAID33nsPiYmJGDRokNhu0qRJRqnh1q1b4u8HDhyAs7MzAKBTp05IS0tDnz59irQjIiLzwbBDZiU7OxsJCQni8v79+wEACoUCmZmZiIuL02mvUCig0Wiwf/9+dO3aVVzv6ekJGxsbg9R09epVAICdnZ0YdAo4OzvDzs4OWVlZYjsiIjIvDDtkVhISEjB69Ogi6zUaTbHrCyu8PTw8HD4+PuWq4dnApdFoxPWXL1/WuQIsNzcX2dnZYrvCYcyQgYuIiMqPV2OBV2OZk2eDxpgxY1DwFV23bh3u3buH0NBQzJkzBy+88ALGjx8PAJDJZDqnvCoSNOLi4koNVvqoSOAiIqLS8dLzMmDYMV+3bt3CyJEjAQBWVlYYMGAAdu7ciSFDhuDrr78WJyZv3boV9evXN8h7Phu40tLSMGPGjFJft2TJEp3TXBzZISIyLn2P3zyNRZUuOTkZarW6zK/Lzc3Fzp07AUD83wL5+flF5vM8T8F9cvTh7OwMKysrMVhZWFggPz9f/F/gaRB7dj5P4cBUHKVSCZVKpXfNRERUPhzZAUd2KlNycjKGDH0Pebk5Jq1DBsDUX3xrK0t8+dVOBh4ionLiyA6ZJbVajbzcHGS/0AaClUPZXpyVBavk82JQyVW1A+zsylyDPDMF1g+uoZ/3Y9Sy1ZbptRm5wP5bdsjVymAlF9CvfhYcy/HEigdP5Ngfbw+1Ws2wQ0RkZAw7ZBI29/5X4X3YJp+v0Ov3x9tX6PXZWmDXX44V2gcRERlfme+gXHCZbXHu379foWJI+pRKJSwsK+/hnebMytICSqXS1GUQEUlemUd22rRpg927d6NVq1Y66/fv34+xY8fiwYMHhqqNJEilUmHXzq/KNEF5+vTpSE9PL3G7k5MTli5dWqY6yjJBGXj6INDnfbdr1apV5geBcoIyEVHlKHPY6dixI15++WXMnz8fM2fOxOPHjzFhwgR8/fXXfOoz6UWlUul9kE9NTRWDzosvvojOnTtj0aJFmDVrFo4fP47ff/8d6enpqFmzJlxcXIxSr1qtLjXEP3jwACqViiM1RERmqMynsdatW4f9+/dj5cqV+L//+z+0bNkSMTEx+P333zFlyhSDF3jv3j0MGTIErq6usLW1RYsWLXD+/L9zNQRBwNy5c1GnTh3Y2tqiS5cuuHHjhsHroMqRnZ2NuLg48WfChAkAnj6qYeTIkVAoFACePiZi5MiRsLW1BQBMmDBB53XPO91a1hoKP3PL3t4ePXr0AAD06NED9vb/zvuZNGmSwWogIiLDKdel51qtFpMmTcL69ethYWGB77//Xue5RIby6NEjtG7dGp06dcK4ceNQq1Yt3LhxAw0aNECDBg0AAIsXL0ZYWBh27NgBb29vfPrpp7h06RL+/PNPvW/oxkvPzYc53L3YHGogIqLSGe0Oyjdv3sSgQYOQlJSEzZs349SpU/jiiy/w4YcfIjQ0tEzzIEoza9YsREVF4fTp08VuFwQB7u7umDZtGqZPnw4A4qW827dvx8CBA/V6H4Yd8/Hs3YtnzJiBtLQ01KtXD3PmzCnSfsGCBUhISICzszOWLFkirq/I3YtLemSFq6srwsLCirSfOXMmHj16ZNBHVhARUemMdp+dVq1a4c0338SRI0fg7OyM119/HT169MB7772Ho0eP4o8//qhQ4YUdOnQIXbt2Rf/+/XHq1CnxWUgffPABACA+Ph5JSUno0qWL+BqlUomXXnoJ0dHRJYadnJwc5OT8e1O7501+pcplY2OjMxoydepUzJ07F3///Tfc3d3h4PDvvXkyMzPFUDJ16lSDjaI8W4NSqURaWhr++ecf1K1bF3aF7u2TlZWFR48eie04kkNEZH7KNWcnIiJC59b4HTp0wB9//IE2bdoYsjbcunUL69evR6NGjXDkyBGMGzcOkydPxo4dOwAASUlJAFBksqtKpRK3FScsLAxKpVL88fDwMGjdZDgBAQGQy59+TXv27Ilx48bh999/x7hx49CzZ08AgFwuR0BAgNFqePvtt8Xfe/TogY8++ggXL17ERx99JM7febYdERGZD7N+XISVlRXatWuHM2fOiOsmT56Mc+fOITo6GmfOnEFAQAASExNRp04dsc2AAQMgk8mwd+/eYvdb3MiOh4cHT2OZqcjISMydO7fE7SEhIQgMDDTa++fm5uKNN94otd3PP/8MKyveQ4iIqLIY7TTWl19+WeI2mUyGoUOHlnWXJapTpw6aNm2qs65JkybYv38/AMDNzQ3A0+ctFQ47ycnJRe4DVJi1tTWsra0NVicZV2BgIEJCQrBmzRqkpKSI61UqFSZMmGDUoAM8Dd0DBw5EREREiW0GDhzIoENEZKbKHHY+/PBDneW8vDxkZWXBysoKdnZ2Bg07AQEBuH79us66uLg41KtXDwDg7e0NNzc3HDt2TAw36enpOHv2LMaNG2ewOsj0AgMDERAQgIsXLyI1NRUuLi7w8/MTL0U3trFjxwJAkcAjk8nw7rvvituJiMj8lDnsFEzGLOzGjRsYN24cPvroI4MUVWDKlCno0KEDFi5ciAEDBuD3339HeHg4wsPDATw90AQHB2PBggVo1KiReOm5u7s7+vTpY9BayPQUCgVat25tsvcfO3YsRo4ciYMHDyIxMRHu7u7o3bs3R3SIiMycwebsnD9/HkOGDMG1a9cMsTvRDz/8gNmzZ+PGjRvw9vbG1KlTxauxgKeXn8+bNw/h4eFIS0vDK6+8gnXr1pXpqhheek5ERFT1GO0+OyWJiYlBYGBglbyMm2GHiIio6jHaBOVDhw7pLAuCgPv372PNmjVGvfyXiIiIqDzKHHaenQsjk8lQq1YtdO7cGcuWLTNUXUREREQGUeawo9VqjVEHERERkVGU+Q7KRERERFWJXiM7U6dO1XuHy5cvL3cxRERERIamV9jZtm0bmjdvDgsLC8hkMpR0AZdMJjNocUREREQVpVfYUavV2L9/P2rXro369evj3LlzcHV1NXZtRERERBWm15ydGjVqID4+HgBw+/ZtTlImIiKiKkOvkZ1+/fohMDAQ7u7ukMlkaNeuXYnPJLp165ZBCyQiIiKqCL3CTnh4ON5++2389ddfmDx5Mj744AM4OjoauzYiIiKiCtP7PjvdunUDAFy4cAEffvghww4RERFVCWW+qeC2bduMUQcRERGRUfCmgkRERCRpDDtEREQkaQw7REREJGkMO0RERCRpDDtEREQkaQw7REREJGkMO0RERCRpDDtEREQkaQw7REREJGkMO0RERCRpDDtEREQkaQw7REREJGkMO0RERCRpDDtEREQkaQw7REREJGkMO0RERCRpDDtEREQkaQw7REREJGkMO0RERCRpDDtEREQkaQw7REREJGkMO0RERCRpDDtEREQkaQw7REREJGkMO0RERCRpDDtEREQkaQw7REREJGkMO0RERCRpDDtEREQkaQw7REREJGkMO0RERCRpDDtEREQkaQw7REREJGkMO0RERCRpDDtEREQkaQw7REREJGkWpi6AiCqfRqPBxYsXkZqaChcXF/j5+UGhUJi6LCIio2DYIapmIiMjsW7dOiQlJYnr3NzcMH78eAQGBpqwMiIi4+BpLKJqJDIyEvPmzUP9+vWxdu1a/Pjjj1i7di3q16+PefPmITIy0tQlEhEZHMMOUTWh0Wiwbt06+Pv7Y/78+cjNzUV0dDRyc3Mxf/58+Pv7Y/369dBoNKYulYjIoHgai6iauHjxIpKSktCrVy8MHTq0yGmsXr164cyZM7h48SJat25twkqJiAyLYYeomkhNTQUAbN68Gf7+/vj000/h7e2N+Ph47Nq1C5s3b9ZpR0QkFTyNRVRNODs7AwCaN2+OBQsWoFmzZrCzs0OzZs2wYMECNG/eXKcdEZFUcGSHqBrSaDSIjY0VLz0vCDqVXQMvfyeiysCwQ1RNpKWlAQAuXbqEN998E7m5ueI2KysrcbmgnTHx8nciqkw8jUVUTbi4uIi/Fw46zy4XbmcMBZe/P3r0SGf9o0ePePk7ERkFR3aIqolmzZpBJpNBEIQS28hkMjRr1sxoNWg0GixfvhyCIKBVq1aoW7cucnNzYWVlhbt37+Ls2bNYsWIFAgICeEqLiAyGYYeomoiNjX1u0AEAQRAQGxuL9u3bG6WGmJgYpKWloWbNmjh37hzOnj0rbpPL5ahZsyYePnyImJgYtG3b1ig1EFH1w9NYRNXEoUOHDNquPGJiYgAADx8+LHZ7wfqCdkREhsCRHaJq4ty5cwZtVx6F78784osvYujQoeK9fr766iv89ttvRdoREVUUww5RNZGdnS3+3r59ewwfPlwMGtu3bxdDTuF2hpaRkQEAsLa2xoIFC2Bh8fQ/QQX3+nnzzTeRk5MjtiMiMgSexiKqhuRy+XOXjaXg7sw5OTn49NNPceXKFWRlZeHKlSv49NNPkZOTo9OOiMgQOLJDJFHZ2dlISEgQlxUKhXh66OzZszqTgwtTKBSIi4sTlz09PWFjY2OQGvLz88Xfz58/j+joaHHZyspKp52haiAikgmlXZ5RDaSnp0OpVEKtVsPJycnU5RAZRFxcHEaPHl3h/YSHh8PHx6fK1kBE0qXv8ZthBww7JE3PjqrExMRg3bp1pb5u/PjxaNWqlbhsyJEdrVaL4OBgZGdnw9HREU2bNsXZs2fx0ksv4c8//0RGRgZsbGywcuVKnVNrHNkhouIw7JQBww5VBxqNBq+//jq0Wm2JbeRyOY4ePWrUG/pFRkZi7ty5JW4PCQnhIyOISC/6Hr85QZmomlAoFPjss8+e2+azzz4z+p2LAwMDERISApVKpbPezc2NQYeIjIIjO+DIDlUvkZGRWL16NR48eCCuq127NiZOnFipQUOj0eDHH3/EsmXLMG3aNPTo0YOPiCCiMpHkyM6iRYsgk8kQHBwsrsvOzsaECRPg6uoKBwcH9OvXD8nJyaYrksjMBQYGIiIiAtOmTQMATJs2DXv27Kn0ERWFQoHGjRsDABo3bsygQ0RGU2XCzrlz57Bx40b4+fnprJ8yZQq+//577Nu3D6dOnUJiYiLefvttE1VJVDUwaBBRdVIlwk5mZiYGDx6MTZs2oUaNGuJ6tVqNLVu2YPny5ejcuTPatm2Lbdu24cyZM+Jt54mIiKh6qxJhZ8KECXjzzTfRpUsXnfUXLlxAXl6eznpfX194enrq3KzsWTk5OUhPT9f5ISIiImky+zsoR0RE4H//+1+xDydMSkqClZUVnJ2dddarVCokJSWVuM+wsDDMnz/f0KUSERGRGTLrkZ07d+7gww8/xK5duwx6Q7HZs2dDrVaLP3fu3DHYvomIiMi8mHXYuXDhAlJSUtCmTRtYWFjAwsICp06dwqpVq2BhYQGVSoXc3FykpaXpvC45ORlubm4l7tfa2hpOTk46P0RERCRNZn0a67XXXsOlS5d01o0YMQK+vr6YOXMmPDw8YGlpiWPHjqFfv34AgOvXryMhIQH+/v6mKJmIiIjMjFmHHUdHRzRv3lxnnb29PVxdXcX1o0aNwtSpU+Hi4gInJydMmjQJ/v7+ePnll01RMhEREZkZsw47+lixYgXkcjn69euHnJwcdO3aVa+HHRIREVH1UOXCzsmTJ3WWbWxssHbtWqxdu9Y0BRGZSHJyMtRqdblf//fff+v8b3nk5eXB0tLSpDUolcoiz9kiIiqMz8YCn41FVU9ycjKGDH0Pebk5Jq1DBsDU/wGxtrLEl1/tZOAhqob0PX5XuZEdInp69/C83Bw8qf8qtDZKk9SgUN+Fzb3/YWzTDLjba0xSQ+JjBTb86Qi1Ws2wQ0QlYtghqsK0Nkpo7Wua5L3lT9IAAO72Gng5mibsEBHpw6zvs0NERERUURzZIarCCkZXTEGWkwHg6akkUzHlexNR1cGwQ1SF2cZHmroEbPjT0dQlEBE9F8MOURX2xDsQWltnk7y3Iu0ObBL/MIsJykREz8OwQ1SFaW2dOUGZiKgUnKBMREREksaRHaIqTJ5d/jsoV5QsNxMAJygTkflj2CGqgpRKJSytrIFbp0xahwymn6BsbWUJpdI0N1YkoqqBYYeoClKpVNj51ZcVfjZWaGgo5syZg3r16pVrH4Z4NlZFa+CzsYioNAw7RFWUSqUyyEG+Xr168PHxMUBFVbsGIpIuTlAmIiIiSWPYISIiIklj2CEiIiJJY9ghIiIiSWPYISIiIklj2CEiIiJJY9ghIiIiSWPYISKT0Gg0uH79OgDg+vXr0Gj4MFEiMg6GHSKqdJGRkQgKCsKyZcsAAMuWLUNQUBAiIyNNXBkRSRHDDhFVqsjISMydOxcpKSk661NSUjB37lwGHiIyOIYdIqo0Go0GixYtAgDUqFEDQ4cOBQAMHToUNWrUAAAsWrSIp7SIyKD4bCwiicrOzkZCQkKJ2//++2+d/y2Jp6cnbGxsDFLDlStXkJWVBTs7O4SGhuLevXvie4SGhmL69OnIysrCoUOH0KxZM4PUQEQkEwRBMHURppaeng6lUgm1Wg0nJydTl0NkEHFxcRg9enSF9xMeHl7uh3SaQw1EJF36Hr8ZdsCwQ9L0vJEdrVaLGzduQK1WQ6lUolGjRpDLiz+rbciRnXXr1iEmJgbdu3fH77//jn/++Ufc5urqinbt2uHIkSNo1aoVxo8fb5AaiEi69D1+8zQWkUTZ2NgUOxoSGRmJdevWISkpSVzn5uaG8ePHIzAw0Kg1+Pv7IyYmBj/99BP8/f0xZMgQeHt7Iz4+Hjt37sSRI0fEdhzJISJD4QRlomokMjIS8+bNQ/369bF27Vr8+OOPWLt2LerXr4958+YZ/Uqot956S/xdq9XqbCu8XLgdEVFF8TQWeBqLqgeNRoPBgwejfv36WLBggc5pK61Wi08++UQcYVEoFEap4Y8//sCUKVMAAHK5XCfgFF5esWIFWrdubZQaiEg69D1+c2SHqJq4ePEikpKSMHjw4CLzc+RyOQYPHoz79+/j4sWLRqshNTUVANCvXz/IZDKdbTKZDP369dNpR0RkCJyzQ1RNFAQIb2/vYrcXrDdm0HBxcQEAdO7cGWPGjMHBgweRmJgId3d39O7dGzdu3MD+/fvFdkREhsCwQ1RNFASI+Ph4nXvYFIiPj9dpZwx+fn5wc3PDrl27sGDBAvTv31/cptVqsWvXLtSpUwd+fn5Gq4GIqh+exiKqJgoHjeImB1dG0FAoFBg/fjyio6PxySefiDcZvHLlCj755BNER0dj3LhxRpszRETVEycogxOUqfoouBrL398fgwcPFi/73rVrF6KjozF//nyDX35eUh3PXv5ep04djBs3rlLen4ikgTcVLAOGHapOzCVoaDQaXLx4EampqXBxcYGfnx9HdIioTBh2yoBhh6obBg0ikgLeQZmISqRQKHgfGyKqNjhBmYiIiCSNYYeIiIgkjWGHiIiIJI1hh4iIiCSNYYeIiIgkjWGHiIiIJI1hh4iIiCSNYYeIiIgkjWGHiIiIJI1hh4iIiCSNYYeIiIgkjWGHiIiIJI1hh4iIiCSNYYeIiIgkjWGHiIiIJI1hh4iIiCSNYYeIiIgkjWGHiIiIJI1hh4iIiCSNYYeIiIgkjWGHiIiIJI1hh4iIiCSNYYeIiIgkjWGHiIiIJI1hh4iIiCSNYYeIiIgkjWGHiIiIJI1hh4iIiCSNYYeIiIgkjWGHiIiIJI1hh4iIiCTNwtQFEBFVVbm5uTh48CASExPh7u6O3r17w8rKqtLrePDgASZOnAi1Wg2lUok1a9agVq1alV5HRZhDX0qhH6l4MkEQBFMXUZKwsDB8++23uHbtGmxtbdGhQwcsXrwYjRs3FttkZ2dj2rRpiIiIQE5ODrp27Yp169ZBpVLp/T7p6elQKpVQq9VwcnIyxkchIonZsGEDvv76a2i1WnGdXC7HgAEDMHbs2Eqro2fPnsjMzCyy3sHBAT/88EOl1VER5tCXUujH6kjf47dZn8Y6deoUJkyYgN9++w1Hjx5FXl4e3njjDTx+/FhsM2XKFHz//ffYt28fTp06hcTERLz99tsmrJqIpG7Dhg2IiIjQOTgDgFarRUREBDZs2FApdZR0gAaAzMxM9OzZs1LqqAhz6Esp9CM9n1mHncOHD2P48OFo1qwZWrZsie3btyMhIQEXLlwAAKjVamzZsgXLly9H586d0bZtW2zbtg1nzpzBb7/9ZuLqiUiKcnNzERERIS43bdoUy5YtQ9OmTcV1ERERyM3NNWodDx480DlAe3l5YeHChfDy8hLXZWZm4sGDB0atoyLMoS+l0I9Uuio1Z0etVgMAXFxcAAAXLlxAXl4eunTpIrbx9fWFp6cnoqOj8fLLLxe7n5ycHOTk5IjL6enpRqyaiKRkz5494u8//PADHBwcAABt27bVGQXYs2cPhg0bZrQ6xo0bJ/5+6NAhcQi/Q4cOSE9Px1tvvSW2++abb4xWR0WYQ19KoR+pdFUm7Gi1WgQHByMgIADNmzcHACQlJcHKygrOzs46bVUqFZKSkkrcV1hYGObPn2/McolIIrKzs5GQkCAu7927FwDg6emJxMTEIu09PDxw584d7N27F/7+/jrbPD09YWNjU+EaAODhw4cAAGdnZyQlJRX5b17BPIaHDx8iLi6uwjUYgqH6siKf4dkaqmI/UtmZ9QTlwsaNG4effvoJv/76K+rWrQsA2L17N0aMGKEzSgMAL774Ijp16oTFixcXu6/iRnY8PDw4QZmIioiLi8Po0aMNsq/w8HD4+PhUyRoMwVCfoyKfwRxqIMPRd4JylQg7EydOxMGDBxEZGQlvb29x/fHjx/Haa6/h0aNHOqM79erVQ3BwMKZMmaLX/nk1FlH1lpycLJ4mf1ZOTo7Ov/bXr1+P1NRUyOVyzJw5E8nJydi6dStGjhwJlUqFxYsXQ6vVwsXFRecUCQC4ubnB2tq62PfJy8uDpaWlXjUAwIoVK5CVlQUAmD59OlJTU8U6XFxcsHTpUgCAnZ2dzn8Ln1eDUqks05WsxamMvnzeZwDK1pfG6EfAMH1JpZNE2BEEAZMmTcJ3332HkydPolGjRjrb1Wo1atWqhT179qBfv34AgOvXr8PX1/e5c3aexbBDVH0lJydj8JChyM8z7oTiqsDK0gJf7dxV7oM0+/JfFe1L0o++x2+znrMzYcIE7N69GwcPHoSjo6OYxpVKJWxtbaFUKjFq1ChMnToVLi4ucHJywqRJk+Dv76930CGi6k2tVvPg/P/l5uVDrVaX+wDNvvxXRfuSDMusw8769esBAB07dtRZv23bNgwfPhzA0yFIuVyOfv366dxUkIhIH0qlEhaWVjxI4+lohFKpLPfr2Zf/qmhfkmGZ9WmsysLTWETV2/PmmTwrLS0NM2bMAAB4e3sjNTVVfLyAi4sL4uPjAQBLliwpcqXo8zxvnklxUlNTMWvWLHFZoVBAo9GI/1tg0aJF4u06SmPsOTvPMoe+NEY/ApyzU1kkMWensjDsEFFJnr1Uedq0acjIyIClpSVWr16NO3fuIDQ0FHPmzIGHhwcmTpyI/Px8ODo6YtmyZTr7MuSl55MmTRKvKl21ahWSk5PFOlQqFSZPngwAsLa2xurVqytcgyEYqi8Neel5VexH+hfDThkw7BBRSczhsm9zqMEQzOGyb3OogQyHYacMGHaIqCTPjgRMnjwZ2dnZcHFxwaJFi4q0nzlzJh49egQbGxusWrVKZ5shR3YKDtj29vZYsWJFkdcEBweLl1SHh4dXuAZDMFRfGnJkpyr2I/1LEldjERGZmo2Njc6/4GfOnIn58+cjNTUVtWvX1plLkpaWhkePHontDPUv/2drAJ4+UuHChQt4/PgxatasqTOfJDU1VTxAt23b1mxGIMyxL6tiP1LZcWQHHNkhIv1pNBq89tpr4rKbmxtGjRqFLVu26Nys7tixY1AoFEar49mncbu4uGDkyJHYunUrUlNTxfWFnzllbsyhL6XQj9WZvsdvs37qORGRuVEoFAgJCRGXk5KSEBoaqnNwDgkJMWrQAQAHBwf4+vqKy6mpqVi6dKnOAdrX19esD9Dm0JdS6EcqHcMOEVEZBQYGIiQkBDVq1NBZ7+LigpCQEAQGBlZKHRs2bNA5UBfm6+uLDRs2VEodFWEOfSmFfqTn42ks8DQWEZWPRqPBxYsXkZqaChcXF/j5+Rl9RKc4mZmZCAsLQ2JiItzd3TF79uwqNxJhDn0phX6sbng1Vhkw7BAREVU9nLNDREREBIYdIiIikjiGHSIiIpI0hh0iIiKSNIYdIiIikjSGHSIiIpI0hh0iIiKSNIYdIiIikjSGHSIiIpI0C1MXYA4KbiKdnp5u4kqIiIhIXwXH7dIeBsGwAyAjIwMA4OHhYeJKiIiIqKwyMjKgVCpL3M5nYwHQarVITEyEo6MjZDKZqcspVnp6Ojw8PHDnzh0+v6uC2JeGwX40HPal4bAvDaOq9KMgCMjIyIC7uzvk8pJn5nBkB4BcLkfdunVNXYZenJyczPqLV5WwLw2D/Wg47EvDYV8aRlXox+eN6BTgBGUiIiKSNIYdIiIikjSGnSrC2toa8+bNg7W1talLqfLYl4bBfjQc9qXhsC8NQ2r9yAnKREREJGkc2SEiIiJJY9ghIiIiSWPYISIiIklj2CEis3H79m3IZDLExMSYuhQikhCGHTMyfPhwyGQyyGQyWFpawtvbGzNmzEB2drbYJjQ0FB06dICdnR2cnZ1NV2wl0adPbt++jVGjRsHb2xu2trZo0KAB5s2bh9zcXKPX17FjRwQHBxv9fcrK3PutOtLnb/LWW2/B09MTNjY2qFOnDoYOHYrExEQTVm2e9OnLAjk5OWjVqhVDdAn06UsvLy+xTcHPokWLTFh12fEOymamW7du2LZtG/Ly8nDhwgUMGzYMMpkMixcvBgDk5uaif//+8Pf3x5YtW0xcbeUorU+uXbsGrVaLjRs3omHDhrh8+TI++OADPH78GEuXLjVx9abDfjM/pf1NOnXqhI8//hh16tTBvXv3MH36dLzzzjs4c+aMiSs3P6X1ZYEZM2bA3d0dsbGxJqrU/OnTlyEhIfjggw/EZUdHR1OUWn4CmY1hw4YJvXv31ln39ttvC61bty7Sdtu2bYJSqaycwkyoLH1S2JIlSwRvb+/ntgEgbNq0SejTp49ga2srNGzYUDh48KBOm0uXLgndunUT7O3thdq1awtDhgwRHjx4INYGQOcnPj6+zJ/RGIzZb48ePRJGjRol1KxZU3B0dBQ6deokxMTECIIgCCkpKYJKpRJCQ0PF9lFRUYKlpaXwyy+/CIIgCBqNRli8eLHQoEEDwcrKSvDw8BAWLFggCIIgxMfHCwCE/fv3Cx07dhRsbW0FPz8/4cyZM2XtArNTnr/JwYMHBZlMJuTm5hq5uqpF37788ccfBV9fX+HKlSsCAOGPP/6ovCKrCH36sl69esKKFSsqtzAD42ksM3b58mWcOXMGVlZWpi7FbOjbJ2q1Gi4uLqXub/78+RgwYAAuXryIHj16YPDgwUhNTQUApKWloXPnzmjdujXOnz+Pw4cPIzk5GQMGDAAA/Oc//4G/vz8++OAD3L9/H/fv34eHh0fFP6QRGLLf+vfvj5SUFPz000+4cOEC2rRpg9deew2pqamoVasWtm7dis8++wznz59HRkYGhg4diokTJ+K1114DAMyePRuLFi3Cp59+ij///BO7d++GSqXSeY85c+Zg+vTpiImJgY+PD4KCgpCfn1+xTjAzpf1NUlNTsWvXLnTo0AGWlpaVXF3VUlxfJicn44MPPsBXX30FOzs7E1ZXtZT0vVy0aBFcXV3RunVrfPHFF1Xv/4+mTlv0r2HDhgkKhUKwt7cXrK2tBQCCXC4XvvnmmyJtq9PIjr59UuDGjRuCk5OTEB4e/tx9AxA++eQTcTkzM1MAIPz000+CIAjC559/Lrzxxhs6r7lz544AQLh+/bogCILw6quvCh9++GE5P53xGKvfTp8+LTg5OQnZ2dk66xs0aCBs3LhRXB4/frzg4+MjDBo0SGjRooXYPj09XbC2thY2bdpU7P4LRnY2b94sriv4V/nVq1f1+uzmSt+/yYwZMwQ7OzsBgPDyyy8LDx8+NFHF5qu0vtRqtUK3bt2Ezz//XBCEf79XHNkpSp/v5bJly4QTJ04IsbGxwvr16wVnZ2dhypQpJqy67Dhnx8x06tQJ69evx+PHj7FixQpYWFigX79+pi7LpMrSJ/fu3UO3bt3Qv39/nfPLJfHz8xN/t7e3h5OTE1JSUgAAsbGxOHHiBBwcHIq87ubNm/Dx8SnnJ6ocxui32NhYZGZmwtXVVWf9kydPcPPmTXF56dKlaN68Ofbt24cLFy6It5y/evUqcnJyxFGekhT+u9SpUwcAkJKSAl9f3+d/aDOnz9/ko48+wqhRo/D3339j/vz5eO+99/DDDz9AJpOZqGrz9Ly+XL16NTIyMjB79mwTV1k1lPa9nDp1qvi7n58frKysMGbMGISFhVWZx0kw7JgZe3t7NGzYEACwdetWtGzZElu2bMGoUaNMXJnp6NsniYmJ6NSpEzp06IDw8HC99v3s6QGZTAatVgsAyMzMRK9evYpMeAT+PQCbM2P0W2ZmJurUqYOTJ08W2Vb46sCbN28iMTERWq0Wt2/fRosWLQAAtra2etVe+O9ScJAv+LtUZfr8TWrWrImaNWvCx8cHTZo0gYeHB3777Tf4+/ubqmyz9Ly+PH78OKKjo4sciNu1a4fBgwdjx44dpijZbJX1uPPSSy8hPz8ft2/fRuPGjSuz1HLjnB0zJpfL8fHHH+OTTz7BkydPTF2OWSipT+7du4eOHTuibdu22LZtG+Tyin+127RpgytXrsDLywsNGzbU+bG3twcAWFlZQaPRVPi9jM1Q/damTRskJSXBwsKiSJ/UrFkTwNMrBocMGYJ3330Xn3/+Od5//31xtKxRo0awtbXFsWPHjPdhqwh9/v9dEPBycnIqs7Qq59m+XLVqFWJjYxETE4OYmBj8+OOPAIC9e/ciNDTUxNWaN32+lzExMZDL5ahdu3YlV1d+DDtmrn///lAoFFi7di0AICEhATExMUhISIBGoxH/z5yZmWniSivPs31ScMD29PTE0qVL8eDBAyQlJSEpKalC7zNhwgSkpqYiKCgI586dw82bN3HkyBGMGDFCDDheXl44e/Ysbt++jYcPH5r16IMh+q1Lly7w9/dHnz598PPPP+P27ds4c+YM5syZg/PnzwN4OrlYrVZj1apVmDlzJnx8fDBy5EgAgI2NDWbOnIkZM2bgyy+/xM2bN/Hbb79Vm9soPKvw3+Ts2bNYs2YNYmJi8Pfff+P48eMICgpCgwYNOKqjh8J96enpiebNm4s/BaecGzRogLp165q4UvNXuC+jo6OxcuVKxMbG4tatW9i1axemTJmCIUOGoEaNGqYuVX+mnjRE/yruEkBBEISwsDChVq1aQmZmZrGXOwMQTpw4Uen1VgZ9+mTbtm3F9klpX28AwnfffaezTqlUCtu2bROX4+LihL59+wrOzs6Cra2t4OvrKwQHBwtarVYQBEG4fv268PLLLwu2trZmf+m5IBim39LT04VJkyYJ7u7ugqWlpeDh4SEMHjxYSEhIEE6cOCFYWFgIp0+fFtvHx8cLTk5Owrp16wRBeHrp+YIFC4R69eoJlpaWgqenp7Bw4UKxLZ6ZSPro0SNJfMdL+5vExMQInTp1ElxcXARra2vBy8tLGDt2rHD37t3KL9bM6fP9LowTlEtWWl/++uuvwksvvSQolUrBxsZGaNKkibBw4cIiFymYO5kgCEJlhisiIiKiysTTWERERCRpDDtEREQkaQw7REREJGkMO0RERCRpDDtEREQkaQw7REREJGkMO0RERCRpDDtEREQkaQw7RFStyWQyHDhwwNRlEJERMewQUbXw2WefoVWrVkXW379/H927d6/8goio0liYugAiotzcXFhZWZnkvd3c3EzyvkRUeTiyQ0SVrmPHjpg4cSKCg4NRs2ZNdO3aFZcvX0b37t3h4OAAlUqFoUOH4uHDh+JrDh8+jFdeeQXOzs5wdXVFz549cfPmTZ393r17F0FBQXBxcYG9vT3atWuHs2fPYvv27Zg/fz5iY2Mhk8kgk8mwfft2AEVPY126dAmdO3eGra0tXF1dMXr0aGRmZorbhw8fjj59+mDp0qWoU6cOXF1dMWHCBOTl5Rm1z4io/Bh2iMgkduzYASsrK0RFRWHRokXo3LkzWrdujfPnz+Pw4cNITk7GgAEDxPaPHz/G1KlTcf78eRw7dgxyuRx9+/aFVqsFAGRmZuLVV1/FvXv3cOjQIcTGxmLGjBnQarV49913MW3aNDRr1gz379/H/fv38e677xap6fHjx+jatStq1KiBc+fOYd++ffjll18wceJEnXYnTpzAzZs3ceLECezYsQPbt28XwxMRmSFTP3adiKqfV199VWjdurW4/PnnnwtvvPGGTps7d+4IAITr168Xu48HDx4IAIRLly4JgiAIGzduFBwdHYV//vmn2Pbz5s0TWrZsWWQ9AOG7774TBEEQwsPDhRo1agiZmZni9v/+97+CXC4XkpKSBEEQhGHDhgn16tUT8vPzxTb9+/cX3n333dI/OBGZBEd2iMgk2rZtK/4eGxuLEydOwMHBQfzx9fUFAPFU1Y0bNxAUFIT69evDyckJXl5eAICEhAQAQExMDFq3bg0XF5dy13T16lW0bNkS9vb24rqAgABotVpcv35dXNesWTMoFApxuU6dOkhJSSn3+xKRcXGCMhGZROFAkZmZiV69emHx4sVF2tWpUwcA0KtXL9SrVw+bNm2Cu7s7tFotmjdvjtzcXACAra1t5RQOwNLSUmdZJpOJp9OIyPww7BCRybVp0wb79++Hl5cXLCyK/mfpn3/+wfXr17Fp0yb83//9HwDg119/1Wnj5+eHzZs3IzU1tdjRHSsrK2g0mufW0aRJE2zfvh2PHz8Ww1hUVBTkcjkaN25c3o9HRCbG01hEZHITJkxAamoqgoKCcO7cOdy8eRNHjhzBiBEjoNFoUKNGDbi6uiI8PBx//fUXjh8/jqlTp+rsIygoCG5ubujTpw+ioqJw69Yt7N+/H9HR0QAALy8vxMfHIyYmBg8fPkROTk6ROgYPHgwbGxsMGzYMly9fxokTJzBp0iQMHToUKpWqUvqCiAyPYYeITM7d3R1RUVHQaDR444030KJFCwQHB8PZ2RlyuRxyuRwRERG4cOECmjdvjilTpuCLL77Q2YeVlRV+/vln1K5dGz169ECLFi2waNEicW5Nv3790K1bN3Tq1Am1atXCnj17itRhZ2eHI0eOIDU1Fe3bt8c777yD1157DWvWrKmUfiAi45AJgiCYuggiIiIiY+HIDhEREUkaww4RERFJGsMOERERSRrDDhEREUkaww4RERFJGsMOERERSRrDDhEREUkaww4RERFJGsMOERERSRrDDhEREUkaww4RERFJ2v8DZZes2kRHzC4AAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# setup a plotting data frame\n", "plot_df = pd.concat([incawrapper_flux_samples_long, gui_flux_samples_long])\n", "sns.boxplot(data=plot_df, x='reaction', y='flux', hue='source')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The plot shows that the MCMC samples appear very comparable. However, we now wish to conduct a statistical comparison between the samples of flux distributions obtained from the INCAWrapper and the GUI. To do this we will use a Hotelling's $t^2$ test, which is capable of comparing multivariate means. The most appropriate implementation that we could find was in R. Therefore the test is conducted the R script, `compare_mcmc_samples.R`. Now we fetch the results of the test." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
value
p_value0.364359
degree_of_freedom_13.000000
degree_of_freedom_25696.000000
n_permutations1000000.000000
\n", "
" ], "text/plain": [ " value\n", "p_value 0.364359\n", "degree_of_freedom_1 3.000000\n", "degree_of_freedom_2 5696.000000\n", "n_permutations 1000000.000000" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pd.read_csv(hotelling_result_file, sep=' ', index_col=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that the hotelling's $t^2$ test find that there is no evidence that samples from the GUI and INCAWrapper come from different distributions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Conclusion\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As the MCMC samples from the INCAWrapper and the GUI was not significantly different (alpha = 0.05), we conclude that the two methods explore the flux space and thus the differences in the results can be attributed to randomness." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\n", "[1] M. R. Antoniewicz, J. K. Kelleher, and G. Stephanopoulos, “Determination of confidence intervals of metabolic fluxes estimated from stable isotope measurements,” Metabolic Engineering, vol. 8, no. 4, pp. 324–337, Jul. 2006, doi: 10.1016/j.ymben.2006.01.004.\n", "\n", "[2] M. R. Antoniewicz, J. K. Kelleher, and G. Stephanopoulos, “Elementary metabolite units (EMU): A novel framework for modeling isotopic distributions,” Metabolic Engineering, vol. 9, no. 1, pp. 68–86, Jan. 2007, doi: 10.1016/j.ymben.2006.09.001.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [] } ], "metadata": { "kernelspec": { "display_name": "bfair-testing", "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.10.8" }, "orig_nbformat": 4, "vscode": { "interpreter": { "hash": "820c70ec08a0eb018d8ec3c5d089748cbb2d1e243fbedf0a7cbeb7c8948c3b84" } } }, "nbformat": 4, "nbformat_minor": 2 }