diff --git a/.github/workflows/ci_push.yml b/.github/workflows/ci_push.yml index 933c74f..ad37df8 100644 --- a/.github/workflows/ci_push.yml +++ b/.github/workflows/ci_push.yml @@ -7,7 +7,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: [3.7, 3.8] + python-version: [3.6, 3.7, 3.8, 3.9] steps: - name: Check out repository diff --git a/.gitignore b/.gitignore index 28df307..395603f 100644 --- a/.gitignore +++ b/.gitignore @@ -8,24 +8,39 @@ yaml2sbml/tests/test_yaml2sbml/sbml_test.xml ### Files generated when running the notebooks (AMICI + PEtab files) ## Lotka Volterra Python -./doc/examples/Lotka_Volterra_python/amici_models -./doc/examples/Lotka_Volterra_python/amici_models/* +./doc/examples/Lotka_Volterra/Lotka_Volterra_python/amici_models +./doc/examples/Lotka_Volterra/Lotka_Volterra_python/amici_models/* -./doc/examples/Lotka_Volterra_python/Lotka_Volterra_AMICI -./doc/examples/Lotka_Volterra_python/Lotka_Volterra_AMICI/* +./doc/examples/Lotka_Volterra/Lotka_Volterra_python/Lotka_Volterra_AMICI +./doc/examples/Lotka_Volterra/Lotka_Volterra_python/Lotka_Volterra_AMICI/* -./doc/examples/Lotka_Volterra_python/Lotka_Volterra_basic.xml +./doc/examples/Lotka_Volterra/Lotka_Volterra_python/Lotka_Volterra_basic.xml -./doc/examples/Lotka_Volterra_python/Lotka_Volterra_PEtab/* -!./doc/examples/Lotka_Volterra_python/Lotka_Volterra_PEtab/measurement_table.tsv +./doc/examples/Lotka_Volterra/Lotka_Volterra_python/Lotka_Volterra_PEtab/* +!./doc/examples/Lotka_Volterra/Lotka_Volterra_python/Lotka_Volterra_PEtab/measurement_table.tsv ## Lotka Volterra CLI -./doc/examplesLotka_Volterra_CLI/Lotka_Volterra_PEtab/* -./doc/examplesLotka_Volterra_CLI/Lotka_Volterra_PEtab.yml -./doc/examplesLotka_Volterra_CLI/Lotka_Volterra_PEtab/*.xml +./doc/examples/Lotka_Volterra/Lotka_Volterra_CLI/Lotka_Volterra_PEtab/* +./doc/examples/Lotka_Volterra/Lotka_Volterra_CLI/Lotka_Volterra_PEtab.yml +./doc/examples/Lotka_Volterra/Lotka_Volterra_CLI/Lotka_Volterra_PEtab/*.xml +## Lotka Volterra Model Editor +./doc/examples/Lotka_Volterra/Lotka_Volterra_Model_Editor/*.yaml +./doc/examples/Lotka_Volterra/Lotka_Volterra_Model_Editor/*.xml +./doc/examples/Lotka_Volterra/Lotka_Volterra_Model_Editor/Lotka_Volterra_PEtab + + +## Sorensen example +./doc/examples/Sorensen/*.xml +./doc/examples/Sorensen/amici_models + + +## FSP examples +./doc/examples/Finite_State_Projection/amici_models +./doc/examples/Finite_State_Projection/*.xml + ## Format Features ./doc/examples/Format_Features/__pycache__ diff --git a/doc/examples/Finite_State_Projection/Finite_State_Projection.ipynb b/doc/examples/Finite_State_Projection/Finite_State_Projection.ipynb new file mode 100644 index 0000000..1c569fe --- /dev/null +++ b/doc/examples/Finite_State_Projection/Finite_State_Projection.ipynb @@ -0,0 +1,160 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Finite State Projection of a Gene Transmission Model\n", + "\n", + "## Model\n", + "\n", + "In this application example, we consider the two-stage model of gene transmission, discussed e.g. in. [[1]](#References). The two-stage model describes the stochastic transcription and translation of a gene by the reactions\n", + "\n", + "\n", + "\\begin{align*}\n", + "R_1 &: \\emptyset \\to mRNA \\\\\n", + "R_2 &: mRNA \\to \\emptyset\\\\\n", + "R_3 &: mRNA \\to mRNA + P \\\\\n", + "R_4 &: P \\to \\emptyset.\n", + "\\end{align*}\n", + "\n", + "The solution to the Chemical Master Equation (CME) [[2]](#References) gives the probabilities of the states at a given time point as the solution of an infinite dimensional ODE\n", + "\n", + "\\begin{align*}\n", + "\\frac{d}{dt}x_{r, p} =& - \\big(k_1 + (k_2+k_3) \\cdot r + k_4 \\cdot p \\big) \\cdot x_{r, p} \\\\\n", + "&+ k_1 \\cdot x_{r-1, p} \\\\\n", + "&+ k_2 \\cdot (r+1) \\cdot x_{r+1, p} \\\\\n", + "&+ k_3 \\cdot r \\cdot x_{r, p-1} \\\\\n", + "&+ k_4 \\cdot (p+1) \\cdot x_{r, p+1}\\\\ \n", + "\\end{align*}\n", + "\n", + "We assume the initial probabilites to be independent Poisson distributions for mRNA and Protein abundances.\n", + "\n", + "## Finite State Projection\n", + "\n", + "The Finite State Projection [[3]](#References) approximates this CME by restricting the states to a finite domain and approximating the probability of the left out states by zero. In our case, we ill restrict the states to $x_{r, p}$ for $(r, p) \\in [0, r_{max}-1] \\times [[0, p_{max}-1]]$.\n", + "\n", + "## Model Construction:\n", + "\n", + "The model is constructed using the `yaml2sbml` Model Editor and exported to SBML. This takes only a few lines of code: " + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from yaml2sbml import YamlModel\n", + "from itertools import product\n", + "from scipy.stats import poisson\n", + "\n", + "model = YamlModel()\n", + "\n", + "r_max = 20\n", + "p_max = 50\n", + "\n", + "lambda_r = 5\n", + "lambda_p = 5\n", + "\n", + "# add parameters\n", + "model.add_parameter(parameter_id='k_1', nominal_value=2)\n", + "model.add_parameter(parameter_id='k_2', nominal_value=1)\n", + "model.add_parameter(parameter_id='k_3', nominal_value=10)\n", + "model.add_parameter(parameter_id='k_4', nominal_value=3)\n", + "\n", + "# add ODEs & construct the rhs\n", + "for r, p in product(range(r_max), range(p_max)):\n", + " \n", + " rhs = f'-(k_1 + (k_2 + k_3)*{r} + k_4*{p}) * x_{r}_{p} '\n", + " \n", + " if r>0:\n", + " rhs += f'+ k_1 * x_{r-1}_{p}'\n", + " if r+1 < r_max:\n", + " rhs += f'+ k_2 * {r+1} * x_{r+1}_{p}'\n", + " if p>0:\n", + " rhs += f'+ k_3 * {r} * x_{r}_{p-1}'\n", + " if p+1 < p_max:\n", + " rhs += f'+ k_4 * {p+1} * x_{r}_{p+1}'\n", + " \n", + " model.add_ode(state_id = f\"x_{r}_{p}\", \n", + " right_hand_side=rhs,\n", + " initial_value=poisson.pmf(r, lambda_r)*poisson.pmf(p, lambda_p))\n", + "\n", + "model.write_to_sbml('gene_expression.xml', overwrite=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Model Simulation\n", + "\n", + "We now use AMICI to simulate the generated model and plot the marginal distributions of the mRNA and protein abundance." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWsAAAEGCAYAAACjLLT8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd3xU15nw8d9zy6hLoC4kUKF3DNjY4IJbEtuJHSdOcbIpXu9uNsmmvFuzm7x5nbaOnTibOInT7MTeOC5x7xUbN8AgUUQVRTQJECBAqJeZ8/4xIzwWaDRCmiY93w/3MzNXd859ZubMw51zzz1HjDEopZSKb1asA1BKKTUwTdZKKZUANFkrpVQC0GStlFIJQJO1UkolACcShebm5pqysrJIFK1U2Kqqqo4aY/LC2VbrrIoHoepsRJJ1WVkZlZWVkShaqbCJyN5wt9U6q+JBqDqrzSBKKZUANFkrpVQC0GStlFIJICJt1kqNdMYYdh1ppbPHC4AgTC5Ix7X1+EdFhiZrpc7Ca9sOc/N97z8h+Q8XV/BfV0+PUURqpNPDAKXOwkNr9pNqJXGpu4BL3QWUWHncv2ofJzu6Yx2aGqE0WSs1SI0tnby27TDlUkypXUipXcg8ZyptXT08vHp/rMNTI5Qma6UG6ekNB/D6DJPsklPrcq0sCq1s7nl7Nz1eXwyjUyOVJmulBumRyjpyrSzGWhnvWz/TruDQyQ5e2HQoRpGpkUyTtVKDsO3QSbYcPMlEq/i0v5VY+WRZafzhzVp0Ug813DRZKzUIj1XVYSOU26cnaxFhulVOdX0TlXuPxyA6NZJpslYqTD1eH4+vPUCxlU+yeM64zSS7hGRxueet3VGOTo10mqyVCtObO47Q2NrJxKATi305YlNhlbBs62FaOnuiGJ0a6TRZKxWm5zceItlyKbHyQ25XahfS7fPxRs2RKEWmRgNN1kqFafXuY+STjS2hvzZ5MpYUy8PLW7RXiBo+mqyVCsOR5k72HWsj3xo74LaWCCWSz7Kth+nq0T7XanhoslYqDFV7jwGQb2WHtf0Eq5CWzh7e3d0YybDUKKLJWqkwrNlzHEcsciQrrO2LrFxcsXl5c0OEI1OjRb+j7omIA9wMXA+MC6yuB54C7jHGxPWINYeaOlhecxhfiGsTROCiybmUjE2NXmAqIa3Zc5xcGTNge3UvR2zGSR4vbW7ge9fOxLIkwhGqkS7UEKl/Bk4AtwB1gXUlwBeA+4FPRTSyIXhqfT3feWITzWF0nfLYFl+6pIIvL51IqkdHjFWna+vqYfOBJmZIxaCeN8Eu4K3mQ2ysb2Lu+DERik6NFqGy0wJjzJQ+6+qAVSKyPYIxnbUTbV1858lNPFt9kHxrDJd6Zvd78QJAl+mh2ruDX762k4fX1PHta6Zx7dxxiOhRkHrP+v0n8PoMBW547dW9Sqx8LISXtxzSZK2GLNRvumMi8gmR9373iYglIp8C4u5a2rd2HOEDP3uT56sPMd+ZyofcC8i2MkmV5H6XMVY6F7vncJXnAnytSXzjofV8/Dcrqa47EeuXo+JI1R5/dc8LoydIsCTxUGBl89ImbbdWQxcqWX8auAFoEJHtgaPpQ8DHAn+LC+1dXm55ejOfu2c13W0u13iWMMeZhBVm2yJAgZXNNe4Sljhz2FbXynW/eod/e2QDh5s7Ihi5ShRr9hwn28ogSdxBP3eCVcDOIy3sOtISgcjUaNJvM4gxZg+BdmkRyQmsi6t+SNV1J/jmQ+upPdrKdLuMBc40HLHPqiwRYbIznlJTyIaenTxetZvnNh7in6+czN8uKdcTRKOU12eo2nucEhk38MZnMMEu5N2eLTxffZCvXT55mKNTo0m4h595xphGEZkW0WjCZIzh16/v5Ppfr6DhmJcPuItY5M4860QdzCMu57rTudZzMTk92fzwua185S9raevScR5Go5pDzbR29YR1McyZpEkKBdZYnq0+OMyRqdEm3GT9QJ/bmPrVazv5yUs1jLcKuda9mHF27rDvI8tK5zJ3Iec603lp8yE+dtdK6k+0D/t+VHw7dTGMDO7kYrBSq4iahmZtClFDMtiLYmLeFvDQ6n3c8cp2JlrFXOKcc1btiOESEWY6FVzunkttQxsfufPtU19eNTqs2XOcNCuJdEk56zJK7UIAntejazUECXUF46tbGvivJzZSbOWxxJ0TtS52JXY+V7uL8XY4fPp37/JIpU6KOhoYY3i39hh5ZA+prmlTiBoOCZOsq/Ye56sPrCVbsljqzh9Ub4/hMMbK4Gp3CXmM5d8ereZHz23BG+rySJXw9ja20dDcQaGVM+SytClEDdVgM15MstPOw8387Z/WkOxL5nL3XFyJzZWGSeLhCvc8ptml/OGt3dx83xpOdsT1VfdqCFbs8nd+KhqOZK1NIWqIwk3W0uc2ag41dfC5u9fQ3WVxhbOIFEmKdgjvY4nF+e4sLnBm8WbNUa7/1Qr2HG2NaUwqMlbsOkqalUSmpA25LG0KUUM1YLIWkduAiwIPLwpaF3FN7d18/o+raWzu5nLnXDKs+BlwaapTypXuIuobO/nwL9/mlS16ldpIYozhnZ2NFJA7bOdGtClEDUU4R9ZXGmNaAHpvgasiF5JfR7eXv7+vkl0NLSx1FpBjhTc0ZTQV2Tlc415Icncaf/+/ldz+4jZ6vDrY/EiwvaGF421dw9IE0qvULkKAR6vqBtxWqb76TdYi8mUR2QhMFZHqoGU3UB3JoLw+wzcfXs/qPcdY4s6LSD/q4ZJhpfIh9wKm2BO4a/kuPnfPao62dMY6LDVEK3YdBfz/IQ+XNElmglXI/Sv30qqT6apBCnVk/QDwEeDpwG3vssAY8zeRCsjnM3z7iY28uOkQ5zozqLDP7jLfaHLEZrE7myXOHNbsPs7Vv3ibtfvibqwrNQgrdjWSaaWSLsPb9DbTqaC5s0e7f6pB6zdZG2OajDF7jDE3GmP2Bi0RuyrE5zN8+8mNPLRmP3PsScx0yiO1q4iY7IznKncxHa0Wn/ztSu5bsQdjtHtfovH6DCt3NVIow3dU3SvfGku+NYZ73t6jXT/VoMRNP+veRP3gan+iPsfpO5R2YsixsrjGvZAi8vh/T2/mmw+t13FFEszmA020dPZQaEWm+W2GXcH+4228vFlnP1fhi4tkfaZEncgTACSJy2XuQuY7U3l6wwE+cuc7p9pAVfx7Z+fw9a8+kwlWIZlWKr9/szYi5auRKebJeqQl6l4iwhxnEle6izh8zMtn/vAuX/pzJfsa22IdmhrAil1HGWulR6xPvyXCdKucdftPULVXz22o8MQ0WfsT9aYRl6iDjbNzuc69hHOcKby25SiX3/EGt724jRbtDRCXunp8rN59jEKJbA+kSXYJyeLy81e363kNFZaYJev3EvW+EZuoezliM9eZzPWepUygiN8s38XS25fz18r9+PQkU1x5bVsDnT0+iiLUXt3LFYe59hTe2nGU+1bsiei+1MgQk2Q9mhJ1sFRJ5iLPPK7xLMFqT+HfH63m2l+9w5o9OuxqPDDG8MtlO8myUimx8iK+v2l2KeOtfP77+W1sPXgy4vtTiS3qybq1s4dvPV496hJ1sDxrDFe5i7nYnceeQ5184rcrufH3q/hr5X6adWComFm+/QibD55kljW4OTzPloiwxJ2DY1y+9sA62ru8Ed+nSlxRG76uq8fHg6v38YtXd3CsrWvUJupeIkKFXcx4q4Ct3j1s2bOflbXVfOeJTVw5s4Dr5xVz8ZQ8PE7MzwGPGne+uoMMK4UKuzhq+0yWJJbYc3nlyGq+/+wW/vv6WaP2O6FCi3iy9vkMT284wE9fqqHuRDtFVjbXeKaRd5Zz2o00rjjMcSYx20zkqDnBLm89r206yHPVB8lKcbl27jg+es445k8Yq1/iCGrp7GHd/hOc78zEjvJY6cV2HrN8FTy4upam9m5u+/hsMpIjNwOSSkwRS9bGGJbXHOHHL2yjpqGZHCuTK93ZjLOGbxSzkUREyJOx5FljOc/MoN53hNrOAzz47n7+vGovhZnJzC8dw6ziLGaNy2JWcRbZaZ5Yhz1iHD7ZSaGVxCR7fEz2v8CZRpJ4eHFjDZvqm/jt3yxgxrjMmMSi4pNEotvQjDnzzMyv/IY1e46RaaUyz55CuTVOk/RZ6DLd7PM1sN/bwHFOctL3Xj/twsxk5pT4E/es4kxKc9LIy0giI8nR9xoQkSpjzMJwtk0qmmyuvfklZjoVkQ4rpAbfMd7sWUu3dPPRc8bxsfklnFeWjWXp5zkahKqzEUnWSUWTzcSbfsMcaxKT7QlR/1k5knWabo75mmg0Jznma+IYTTT5Wt83hU+SY5GbnkR+ZhIFGcnkZSSRl5FETrqHNI9DsmuT6vEvvfdTPDaproPHsbAtwbUl4RP+YJJ1StEU88//sCVmsxAFazedrO2pYa/vAF3GS3FWChdPzWVSfgaT8tMpzU4lK8UlM8XF1iQ+ooSqsxGpmRlOMp/JXBoXFX+kScFlDLlU8F4/4C7TQ6P3JM3edtpMJ22+TtraOjjS0sk+WmjzNdLuG3wvE0vAtiwcS3AswbYE27KwxH8VngRuLQFO3fev73VqiiHpfSwE/r1P3/8XJMqTEhWkZJCZGh/1NYUkrmAO3WYmu7sPsaOtnicrD9HuO32kvlSPTZJj49qCx7awLDn1mYj4m9f6fgYQ/fdXDV1EjqxFpBmoGfaCIyMXSISBOxIlToifWEuNMWF1mI6TOhsP75vGENsY+q2zkTqUqAn352esiUhlIsSaKHFCYsUaJOZ1Nh7eN40hfmLoSxuTlVIqAWiyVkqpBBCpZP37CJUbCYkSa6LECYkVa694iFlj8NMYziAiJxhzc3NNWVnZsJer1GBUVVUdDfcEo9ZZFQ9C1dmInGAsKyujsrIyEkUrFTYR2RvutlpnVTwIVWe1zVoppRKAJmullEoAMblka+WuRnYfbaWls5uWjh5aOr3++53++80d3RRlpnDnjfNwbP3/RCmlop6s9zW28Zk/rHrfWBYecfCIgysOjnEQI6wzB1m8Joe/Ob802iEqpVTciXqyfmxtHQDXei4iQ1JxsE8bMMgYw0vdq/jpS9v5yNxxZKXo2L5KqdEtqm0MPp/h0co6iqxcsq1MXDnzUJ4iwrnODJrau/jVazuiGaJSSsWlqCbrd3cfo76pnUl2yYDb5lhZTLLH86d39rD7aGsUolNKqfgV1WT92No6POIwwSoMa/v5zlQsY/Oj57ZGODKllIpvUUvWrZ09PFd9kFKrCEfssJ6TIknMtifx6tYG3t4R6xETlVIqdqKWrF/cdIj2bm9YTSDBZthlZFqpfP+ZLfR4fRGKTiml4lvUkvUjVXVkWqnky+BmNbfFZoE9je2Hm3lozekzZSil1GgQlWRdd7yNVbWNTLRKzmpevwlWIYVWNj99aTtN7YOfnkoppRJdVJL142vrAZhoF5/V84O78v1ymXblU0qNPhFP1sYYHqmso8jKIV1Sz7qcU135Vuyh9kjLMEaolFLxL+LJunLvcfYfb2PiIE8snsl8Zyq2sfnv57cNQ2RKKZU4Ip6sH62swyM2pWH2rQ5Fu/IppUariCbr9i4vz1QfZIJVhCvDMwyJduVTSo1G/SZrEZkTdN8Vke+IyNMi8t8i4TU+v7T5EG1dPYPuWx1KcFe+B7Urn1JqlAh1ZH1v0P0fA5OAO4AU4LfhFP5IVR0ZVgoFkn3WAZ5Jb1e+n728neYO7cqnlBr5QiXr4A7RlwN/b4x5A/hnYN5ABR840c6KnUepkLPrWx2KiLDQmcHxti7uWr5rWMtWSql4FCpZZ4nI9SLycSDJGNMNYPzToQ84JfoT6+oxMKxNIMFyrSwmWsXc89Zu6o63RWQfSikVL0Il6zeAa4EPA6tEpABARAqBAbtiPFJZR6GVTYZ19n2rBzLfnYrPB7e/WBOxfSilVDzot4uGMeamftYfwt8s0q+2Li9HG1tZ4kwcYnihpUkKM6wKnt6wk5uWlHHOhMGNO6KUUokirK57IjIt+HYgx9u6cMWmzC4aSmxhme1MJNVK4gfPbsXfQqOUUiNPuP2sH+hzG9KJtm4mSOGw9a0OxRWHedYU1u47zgubDkV8f0opFQuDvSgmrG4dPmOG5fLycE2yx5NtZXDr89vo7PFGbb9q9DLGcOBEOweb2jnW2kVbV0+sQ1IjXEQOfW0siqycSBR9RpYIC+zpvHJ8Nf+7Yi9/f3FF1PatRg9jDJvqT/LcxoM8V32Q/X16IZXlpPGBmQVcOaOA+RPGYlvD22VVjW4RSda5nvRh71s9kGI7j2JvHncu28ENC0oYm+aJ6v7VyLb/WBv/+sgG3t19DAuhyMphkVOGjY0XL914aTjeyN1v7ub3b9ZSlJnCVy6byCcXlpDkhDeNnVKhDDZZh3UGT8JrLRl2C53pPN35Jr9YtoNbrp0ZkxjUyGKM4cHV+/nhs1vo6RHOc2ZQYReTLGc6GJhEl+mmzneEbS27+b9PbuKXy3bylaUVfPb8Ulw7qvNTqxEm3NojfW7j0lgrg8n2BP68cq+Oea2GrLGlky/+aQ3/9cRGxnjHcq3nYmY45f0kaj+PuFTY47jKXcwH3EVYranc8swWPvQ/b7Fip44Uqc7egMlaRG4DLgo8vChoXVw6x5mCjc2tL+iY1+rsvVvbyFU/f4u3dzSyyJnJle55pEtK2M8XEcbZuXzIcwGXuws5cszLZ+5+l6/+ZS0NJzsiGLkaqcJpBrnSGPMfAMaY3sPVq4D/iFhUQ5AiScyyJ/LKlhpW7mrkgonRO9GpEtfOwy1886F1lOWm0d7t5Q9v1pIhaVzjnke2lTmkssfbBRRZuWzy7uLFTbt4Y/sR/u+Hp/PJheOjfm5HJa5+k7WIfBn4ClAhItVBf8oA3ol0YEMxwy5nh28f339mC098dTHJrp7gUaH5vMKy6mM0+w4AUG6PY7Eze9iuFXDEZp4zhXJrHCu7N/Ifj23k6fUH+fHHZzM+O3JDMqiRI1RNfAB4AbgV+FbQ+mZjzLGIRjVEjtgstGfw+qEqvv7gOu767HwcPbmjQsjxpPFxz+X0GC+ddJFKckSOerOsdD7onk+Ndx9rardx5c/e5FtXTeXzF5RhaVc/FUK/GcwY02SM2WOMudEYszdoietE3avULuQ8ZwYvb2ngv57YqJeiq7A4YpMmKRFtnhARpjmlXOe5mFxfNrc8s4VP/HYlu/SkuAphRB9uznDKmWtP5q+Vdfz4RT3hqOJLmqRwuXsuF7pz2bS/mat+/ha/XLaDrh6drk6dbkQna4B5zmSm2aX87o1afveGTlSg4ouIMMku4TrPJYyjgDte2c7Vv3iLNXsS4gesiqIRn6xFhPOcmZTbRdz6wjb+qvM2qjiUKsksdedzubuQhkYvn/jtSv7lrxu0m586JfLD4sUBS4QLnXl0mW6+9Xg1WakuH5xZGOuwlDrNeLuAQiuHDT07eXLdbp7feJB/umwSN19Yrr2aRrkRf2TdyxaLpe4CcmUMX3tgHSt3NcY6JKXOyBWHhe40rnMvJs+by09equGS25dz/6q92p49io2aZA3+L8Fl7rmkmVRuvq+STfVNsQ5JqX5lWmlc5lnIB91FmNZkvvPkJi79yXIeWr2Pjm4dCni0GVXJGiBZPFzhLsLucfncPatZv/9ErENSKqQiO5er3MVc4Z5LZ7OHbz2+kcW3vsb/vLKdI82dsQ5PRcmoS9YAaZLMFc55dHYIH/31O9z0pzWs23c81mEp1S8RocTO5xp3CR9wF5HWMYZfLNvB4ltf4yt/qWLZ1ga6vdpEMpKNihOMZ5JlpfNR9xK2evewcsduXq9ZwYWTcvn65ZM5rzw71uEpdUa9A0SNs3Np8rWwzbuX1zcf4PmNh8hO9XD1nEKunFHIBRU5eJxReSw2Yo3aZA3+4SznOpOZYcqp8e5lbW0tn9y5kkXl2Xz98sksnpijA+2ouJVlpbPImslCM51632F2ddTz8Lv13L9qH6keh4sn53LBxBwWVWQzJT9DL2dPcKM6WfdyxWGWM5Fppozt3n1s2ruLz979LueMH8OXLqlg3vixFGQmaeJWcckWiwl2IRPsQnqMl4O+o+zzNvDO1qO8uNk/iXRWssvskiz/UpzFlIJ0JmSn6dF3AtFkHcQRmxlOOVPMBHZ669hUv5N/vH8tAJlJLtOKMphWlMHUwgymFWYwpSCDjGQ3xlEr9R5HbMbbBYy3CwBo9rXRYI7R0H2MrbVNrNhZiy8w4ZMlwoSxqZTnpVI8NoXiMamMG5NMXkYSeelJ5KYnkZXi6hF5nNBkfQaO2ExzSplsxnPEHOe4r5njPc3s2dfM+r31dJn3ZrIel5VCQVYSGcku6Uk26UkO6Uku6clO4LFLqsfGtgTHEqzArR28iJz6QgjgP4AXRIKm6JHTJ0vre6Afq+nUVPzKsFLJIJVJdgkAXuPluGmhyTTT5Gul6UQLG0+0sZITdPi6T3u+ABnJLmNSXDJTXNKT/XU6LckmxbVJDiwexyIpsLi2f3FswbUFSwTHsrAt/38QduB7YIlgiX+dv66/V+dFgu/7Iwn+PpyKL+hLIO9b3/97kqjfk4gk69RUuOCCSJQcbRaQE1j8jDE0drRT39pMfUszB1qbaWnu4sSJbjp97XR4e2jv7qGjpye8CStVXEhJGSl1diA2kBVY3q+jp4fjne2c7OrkZFcXTV2dtHV30drTTWt3N+0d3TS2eTngbafT10O310un10un10ePT3uiRJpEYuhQEWkGaoa94MjIBRJhcrxEiRPiJ9ZSY0xeOBvGSZ2Nh/dNY4htDP3W2Ug1g9QYYxZGqOxhJSKViRBrosQJiRVrkJjX2Xh43zSG+ImhLz0VrJRSCUCTtVJKJYBIJevfR6jcSEiUWBMlTkisWHvFQ8wag5/GcAYROcGYm5trysrKhr1cpQajqqrqaLgnGLXOqngQqs5G5ARjWVkZlZWVkShaqbCJyN5wt9U6q+JBqDqrbdZKKZUAYnIF48Gmdo61doXcZlxWCmPTPFGKSCml4lvUk/XRlk6W/mQ5nQNMT5SR5LDsXy4hPzM5SpEppVT8inqyfmr9ATp7fCx2ZpMkZz5y7sHLO50buOPl7dx2w5woR6iUUvEn6sn6kco68qwspjgTQm7X6Gvir5W7+eKSMqYXZUYpOqWUik9RPcG4+UAT2w6dpMIqGXDbuc5kksTlh89uJRLdC5VSKpFENVk/VlWPLRbl9rgBt00Slzn2ZN7ZdZTlNUeiEJ1SSsWvqCXrbq+PJ9bVUyL5JPfTVt3XVLuULCuVHz67lR6dDFQpNYpFLVkvrznC8bauU4Ogh8MWi/n2dHYdbeGhNfsjGJ1SSsW3qCXrR6v2k2olUWyFdfXvKROsAgqtbO54eTvNHafPZKGUUqNBVJL1sdYulm09TLmMw5LB7VJEWOjM4HhbF3ct3xWhCJVSKr5FJVk/vb6eHp9h4iCaQILlWllMtIq5563d1B1vG+bolFIq/kUlWT9aVU+OlUm2dfb9pee7U/H54PYXYz3zklJKRd+AyVpECkRkfmApGOwOth06yaYDTUwMo291KGmSwgyrgqc3HGD9/hNDKksppRJNv8laROaJyCpgOXB7YHlDRFaJyPxwd/BYVR0WQkUYfasHMtuZSKqVxA+e2aIXyiilRpVQl5vfC3zJGPNu8EoROR/4EzB3oMJ7vD4eX3uAEiufZEkaUqAArjjMs6awYt9GXth0iKtnFw25TKWUSgShmkHS+iZqAGPMKiAtnMLf3HGExtbOsz6xeCaT7PFkWxnc+vw2Onu8w1auUkrFs1DJ+gUReU5EPiUiiwPLp0TkOeDFcAp/rKqeFMtDiZU/PNEClggL7OnsP97G/64IeyIQpZRKaP02gxhjvi4iVwHXAcWB1fXAr40xzw9U8Im2Ll7e3MAkmYA9yL7VAym28yj25nHnsh3csKBEJylQSo14IYdINca8ALxwNgU/U32Qbp+PSZ7hawIJttCZztOdb/KLZTu45dqZEdmHUkrFi7M65BWRAadpf6Syjmwrg2yJzFjUY60MJtsT+PPKvdQeaYnIPpRSKl6E6rqX3c+SA1wdqtDOHh/VdSeYaJUgIsMedK9znCnYWNz6wraI7UMppeJBqGaQI8BeIDjbmsDjkGcMj7d24UGosItDbTZkKZLELHsSr2ypYVVtI+dX5ER0f0opFSuhmkFqgaXGmPKgpcIYUw40hCr0eFsXxVYeKcPQt3ogM+xyMqwUfvDMFnw+vVBGKTUyhUrWPwfG9vO320MVOpRBmwbLEZtz7KlsPniSJ9bVR2WfSikVbaG67v06xN9+GapQQRg/jH2rB1JujWOrtZvbXqzh6tlFpHjsqO1bjS7GGHYcbuGVLQ1U7T1OkmORkeyQmewyv3Qsl07N1/qnIiKs2c1FZJoxZlvv7UDb57jp2BK9Cts75vULzSu5+61avnb55KjtW40OXp/hvhV7+NM7e9gfGKZ3rJWOCHSZHjpMF3e/vZtkx+byGfncsKCEpVPyInqCXY0uYSVr4AFgftBt6EKH+SKYcBRY2ZRahdy1fBefOnc8+ZnJUY9BjUy7jrTwb49Us3bfcQqtbC5wJjLezidV3qtjPmNo8DWyx3eQ1zYd4rnqg8wtGcO/fXAqSyblaNJWQzbYrBrXNW6BM42ubh8/e2V7rENRI4DPZ7j7rVqu+vlbbNnfwkXuXD7ons9UZ8L7EjX4h0EosnO5wJ3NJ9zLWezMZteBDv7mnnf59O9Xse3QyRi9CjVSRP8QOIIyrTSm2mX8tXI/Ww/ql0OdvcaWTm66dw0/fG4rhSaPaz0XM9EO77oBSyymOBO43l3KImcGG/Y2c/Uv3ubW57fS1tUThejVSDSikjXAXGcyHnH50XNbdcxrFbbtDS18+4mNPLGujle2NHDVz9/i7R2NnO/M4lJ3wWlH0uGwxWa6U85H3aVMtEr43Zu1XP7TN3m95nAEXoEa6cJts+4V99kvSVzmWJN5e+cWlm8/wqVTo9crRSUwn/DI6gP85d19AIyx0rjGPW9IU9H1ShYPS9w5TLJLWNWykZv+tIZPLCjhOx+eQVaKO+Ty1egQbrKWPrdxbapdSo1vDz94Zgvnllf9HoUAABwLSURBVGWTnjTY/5PUaJPjSePTng9wwjTTZFoosfJxZXjrTYGVzYfdC9nQs4PHqmp5o+YIP75hNpdNG/RseWoUCmcOxtuAiwIPLwpaF7dssVhkz2bP0Tb+7r5KOrp1kgI1MEuEbCuTcnvcsCfqXrbYzHencbVnCT1tHv723kr+5a8baGrrjsj+1MgRTpv1lcaYFoDeW+CqyIU0PMbZuSxx57KqtpGvPbiOHq8v1iEpdUqulcU17hLm2JN4Ym09V/zsDV7dEnIUBzXKhRp178sishGYKiLVQctuoDp6IZ69iXYxi5yZvLKlgX9/tFrHDlFxxX+UPZVrPEvwtnn4u/+t5Kt/Wcvh5o5Yh6biUKjfeg/gn3jgVuBbQeubjTHHIhrVMJrulNFFN4+v205misv/+8gMvUBBxZUcK4tr3AvZ5N3Fi5t28sb2I/zn1dO48dwJWJbWVeUXamyQJqAJuDF64UTGHHsSXaabe1fsJjPF5Z+vnBLrkJR6H1ss5jqTKbOKWNW9kW8/sYmHVu/nlmtnsKA0O9bhqTgw4vpZn4l/7JDpTLZLuHPZDu55e3esQ1LqjLKsdD7gns9F7lxqD3bw8d+s5OsPrqP+RHusQ1MxNmr6tIkIFzhz6DI9/ODZLWQkO3xy4fhYh6XUaUSEiXYJE6xCNvXs4vnqWp7feJAbz5vAVy6dSFFWSqxDVDEwKo6se1kiXOzOY5yVy7ceq+bFTQdjHZJS/XLF4Rx3Ktd7llIh4/nLqn1ccvtyvvvUJvYcbY11eCrKRlWyBv8Z+EvdBeTKGL72wHruX7WXzh7th63iV5qksNidzfWepZRSzP0r93HpT5dz831reGfnUe3lNEqMmmaQYK44XO6ex2s9a/jOk5u4c9lOvnRJBTeeN55Uz6h8S1QCyLBSWWLN4RxnCjU9e1lRs49lW9+lOCuFjy8o5vr5JZTnpsU6TBUhEonBjkpLF5qbbqoc9nKHmzGGA76jbPTu5JDvGGNSXG6+sJzPLy7TMRtGABGpMsYsDGfb8vKF5gtfiP86G6zHeNnrO0itt54DvqMYYHphJlfOyOfy6QXMLs7Srn8JJlSdHdXJOliD7xibenax33eYVI/D5y8o5eYLy8nLiPykvyoyRnqyDtZqOtjtPcB+XwOHfccwwJgUD+dXZHNueTbnlWUztTADjzPqWj4TSqg6q7/5AwqsbAo82TT6mtjk3cXv3tjFH9/ezXXzxjF3/BimFGQwtSCDrFQ94lbxJ02SmeVUMIsKOkwX9b7DHOxqZMXWRl7cfAgAxxKmFGQwqziTyfkZVOSlUZ6bxvjsVFxbk3i802TdR46VxSXWfM6xW9no3cVTaw/ySFXdqb/npScxtTDDvxRkMLkgnbyMJNKTHNKSHK30KuaSxcNEu4SJdgkAraadw77jNPpOcqyhiWcaDtPue69OC5CbnkTx2BTGjUkmNz2J3PQkctI9jEnxMCbVJSvFPVXH05Mckl1LrwSOMk3W/ci00lhizWGxmU0rHZzwNXPCNHO8rYVttc2s2rWXHnP64FAe2yLV45CWZJOe5JCe7JDqsXEswbYsbAvs3vvCqXVWoOL7b967L6fW+WeND9b3u6JfHXUmaZJCuZ1CuT3u1LpO00WTaaXJ10KLaae1vZ2Gtnb21DfTQSMdvtCjAAqQ5NgkORbJro3HsfA4FkmOhWtbuI7gWhaOLf76LoIVuLUt/31L/PXeX8/lffX91OMzfB/eF4cExyRnXD/Q60gUEWmzFpFmoGbYC46MXOBorIMIQ6LECfETa6kxJi+cDeOkzsbD+6YxxDaGfutspI6sa8I9sRNrIlKZCLEmSpyQWLEGiXmdjYf3TWOInxj60gZWpZRKAJqslVIqAUQqWf8+QuVGQqLEmihxQmLF2iseYtYY/DSGM4jICcbc3FxTVlY27OUqNRhVVVVHwz3BqHVWxYNQdTYiJxjLysqorEzcq8HUyCAie8PdVuusigeh6qy2WSulVALQZK2UUgkg6lcwGmO4+63dA05TNKMok0+eqzO5KKUUxCBZV9c18aPnt+IRu9+xBYwxdBkvxWNTWDIpN8oRKqVU/Il6sn60qg5HLG7wXI5HzjyCXY/x8lT3G/zg2a089/ULsXVMXqXUKBfVNuvOHi9PrT/AeCnsN1EDOGIz357GtkMneWxtXb/bKaXUaBHVZL1s62FOdnQzKTB0YyhlVhH51hh+8mINbV09UYhOKaXiV1ST9aNVdaRZyRRZA7dDiwgLnekcaenkd2/URiE6pZSKX1FL1keaO1lec4QKKT41dvNA8q1syuwifvdGLQ0nOyIcoVJKxa+oJeun1tfjM+bU7BXhWmBPo9tr+OlLsR5qWCmlYicqydoYw1/X1JFvjWGMlT6o52ZYqUyzyni0qo7NB5oiFKFSSsW3qCTrzQdOsv1wMxXW4I6qe81xJpFkufzw2a1EYuAppZSKd1FJ1o9W1WGL9b454AYjSVzmWlNYWdvIa9sOD3N0SikV/yKerLt6fDy5rp7xUkBSiL7VA5lqT2CMlcYPn91Kt/f0iWqVUmokG1SyFpHswe7g9ZrDnGgPr291KJZYzLens7uxlQdX7xtSWUoplWj6TdYi8p2g+zNEZDtQJSJ7RGRRuDt4tKqOVCuJcWH0rR7IeCufIiubn728g5Md3UMuTymlEkWoI+uPBd3/CfANY0w58Engf8IpvLGlk9e2HaZcirFk6C0u/gtlZtDU3sWvX9855PKUUipRhJtBxxljXgAwxqwGUsJ50lPrD+D1mSE3gQTLsbKosEr441t72H+sbdjKVUqpeBYqWVeIyNMi8gxQIiKpQX8L60zhI5V15FpZjLUyhhRkX/PdqRgDt724bVjLVUqpeBVqiNTr+jy2AESkAPjNQAVvOXCSrYdOssiZOYTwzixNkplpTeTZ6h3ctOQ4C0rHDvs+lFIqnvR7ZG2MeaPP0hJY32CM+fVABT+2tg4LOeu+1QOZ5VSQaiXxg2e36IUySqkR76zO+onI70P93QBPrK1nvFVAsnjOKrCBuOIwz5rC+v0neG7jwYjsQyml4kWornvZ/Sw5wNWhCm3u6OZYW9egB20arEn2eLKtDH78/DY6ur0R3ZdSSsVSqCPrI0AlUBW0VAaW/FCFHm/rJsXyUGLlDVecZ2SJsNCeQd2Jdu5bsSei+1JKqVgKlaxrgaXGmPKgpSLQ17ohVKHN7d3D1rd6IOPsXMZb+dy5bCeNLZ0R359SSsVCqGz6c6C/bha3hyrUwLD2rR7IAmca7V1efrFsR9T2qUY3YwyNLZ20dPboCW4VFf123QvV48MY88tQhSaJQ7aVOZS4BmWMlcFkezz3r9rH5y8oY1L+4MbMViocXp+hcs8xXt3awMubD7P3WCsAlkB6ksuC0jFcM2ccV84oICvl7ActU+pMQvWzPkVEphljtvXeDrT9WDdt6JEN0jxnCru7DnDr81u554vnRn3/amTbevAk//rXDWw+eBJbLAolh4XOeAC6TA8d3V1U7jjC6zUbcCzhAzML+eYVk5lSMLwXhKnRK6xkDTwAzA+6jTspksRsexLLtm1jxc6jLJ409IGjlOrq8XHX8p388rWdJOFyoTuXUqsQV07/6hhjOGqa2O09wKub9vPCxoNcN28c37hiCuW50T+AUSPLYM8AhjfTbYxMt8vIsFL4wbNb8fq0HVENTe2RFj7663f4+as7KJUirnUvYZJdcsZEDf6BxvKsMZznzuBjnkuZaU/k2Q2HuOKON/jxC9to6+qJ8itQI0nUJsyNBkdszrGnsfXQSR5fWxfrcFQCae/y0hM0qcVT6+v58J1vU3uonUvdBVzsnjOoC7ySxcNCdxof81xKuVXMb9/YxRV3vMlr20J2pFKqX+E2gySMcquIbdZubn+xhmvmFJHqGXEvUUXAziMtzL7lZeaNH0NGssPLWxoosMZysecc0iSsQSbPKFWSudCdyyS7hHebN/G391by4TlFfO/ameSkJw3jK1Aj3WCPrOO+bcE/5vV0jrR08vs3a2MdjkoQY91Uynzj2b6nm9e3HmWOPZEPuucPKVEHK7Ry+LB7Eec4U3i++hBX/OxNntdhEtQghHvYKX1u41q+lU2ZXcRdr+/iwkm5LCwb9GxkapRJsV0Wuf4RIo0xiAx/VbfFYq4zmQlWAe90VPOVv6zlqlmF/OCjs8jVo2w1gAGPrEXkNuCiwMOLgtbFtUXOTFJMCl/80xo2H2iKdTgqgUQiUQcba2VytbuY+c5UXt58mCvueJNnqw/oxTUqpHCaQa4MGh61JbDuqsiFNDxSJIkr3UXQ7fC5u1dTe6Rl4CcpFSWWWMxxJvFh90LczhT+6YF1fPn+tRw+2RHr0FScCjXq3pdFZCMwVUSqg5bdQHX0Qjx76ZLClc4i2jvgs39YzYET7bEOSan3GWtlcJW7mAXONF7ZcpjL7niD+1ftxaddT1UfoY6sHwA+AjwduO1dFhhj/iYKsQ2LLCudK5zzaGzu5rN/eFcHe1JxxxKL2c5ErnMvJqM7i+88uYkbfrtSm+/U+4SaKabJGLPHGHOjMWZv0HIsmgEOhxwri8vchew/1s7n71nNyY7uWIek1GkyrTQ+4C7iQncuW+pa+PCdb/Ofj2/UAwwFjLCLYkIptHJY6ixg68Fm/u7eSp2sQMUlEWGSXcL17qVMt8t5ePV+Lr59OXct30lrp14BOZqNmmQNUGLnc6E7jzV7jvHl+9fSHXTFmlLxJElcznNncK3nYrJ7srn9xRouvO11/vBmLe1deqAxGo2qZA1QYY/jfGcWr9cc5l//ukFP5Ki4NsZK53LPuVztWUxKRyY/en4ri3/8Gj97uUZ7jowyo/Ja7KlOKZ1089SGGvYda+Prl09m6dS8iPevVeps5Vtj+YBnEQ2+Y2zq2MUvX9vJXct38ZG54/jkwvEsKs/GsrT+jmSjMlkDzLYnkoyH6vod3HTvGqYXZvLVyyZy1awibK30Kk4VWNkUeLI56Wtlq3cPz62v44l19RRlpvCxBeO4ZvY4phdl6IHHCCSRuGqqtHShuemmymEvNxK8xkett57Nvl2c8LVSmp3GVy+dyEfPKcbjjLpWohFFRKqMMQvD2ba8fKH5whcSo84G6zFe9vkOsctbzwHfEQxQmJnMlTMKWDo1j4Vl2TprTQIJVWdHfbLu5TOGfb5DbPTupNF3koKMZP5xaQWfPncCKR471uGpszAaknWwdtNJnfcw+30NHDRH6TZeBJhamMmi8rHMKs5iVnEWk/LTcW09EIlHoersqG0G6csSocwuotQqpN53hE2tu/jeM1u4/cUaphRkMLUwnSkFGYH7GeRnJOlPTRVXUiSJyc54JjOeHuPlqDnBId8xGg438kBDHd1mLwAe26I0J41J+WmU56ZRmpPKuDEpjBuTQlFWsg4rHKf0U+lDRCix8ymx82nwHWOP9yCHDzSz48Bh2nzvTWiQkeQwtTCDKYUZ5KUnkZ7kkJ7skJbkkJ5kk+bpve+Q6rGxLcGxLCwLbEv8i/hvNemr4eaITaHkUGjlAJPxGcNJ08ox08Qx30majraw8mgzL/ka8PUZ+TjFtclJSyIvw8OYVJesFP/yXv12SPU4JLsWKa5NsmvjcSw8toXHsXBtf113bMG1LSwRHEuwguq9ZYEgWOI/UBKJ/ABaiS4iydp1oagoEiVHVxHZzOO94VXbvJ0c6WrhSHczR7uaaWhoZkvdQdq8Q7sisrfCAv5K2zsSrfjHpJVTD99fmfvWba3qZ89xRkad7Z9QTDqQDhSfWus1Ppp7OjjZ085JbzvNPR20erto7e7k+NEuDvk66TQttHu76fD2RHxAezlV5+VU3e/7fei7/an7QX89m7wf79+fiLRZi0gzUDPsBUdGLnA01kGEIVHihPiJtdQYkxfOhnFSZ+PhfdMYYhtDv3U2Us0gNeGe2Ik1EalMhFgTJU5IrFiDxLzOxsP7pjHETwx96SlhpZRKAJqslVIqAUQqWf8+QuVGQqLEmihxQmLF2iseYtYY/DSGM4jICUallFLDS5tBlFIqAWiyVkqpBDDsyVpEPiQiNSKyU0S+NdzlDxcR2SMiG0VkvYjE1aAQIvJHETksIpuC1mWLyCsisiNwOzaWMfbqJ9ZbRKQ+8N6uF5GrYxnjQGJRZ+PhMxaR8SLyuohsEZHNIvKNaMchIskislpENgRi+F5gfbmIvBv4TB4WEU+kYgjszxaRdSLybCz2H45hTdYiYgO/Bq4CZgA3isiM4dzHMLvUGDMv3vpTAvcCH+qz7lvAMmPMZGBZ4HE8uJfTYwX4n8B7O88Y83yUYwpbDOvsvcT+M+4B/sUYMwM4H/hq4LVHM45O4DJjzFxgHvAhETkfuA1/HZoEHAdujmAMAN8AtgY9jvb+BzTcR9bnATuNMbXGmC7gIeC6Yd7HiGeMeRPoOzHxdcB9gfv3AR+NalD96CfWRBKTOhsPn7Ex5qAxZm3gfjP+ZFUczTiMX0vgoRtYDHAZ8Gg0YhCREuAa4O7AY4nm/sM13Mm6GNgf9LiO4IEI4osBXhaRKhH5h1gHE4YCY8zBwP1DQEEsgwnDP4lIdeDnflw02fQjnupszD5jESkDzgHejXYcgSaI9cBh4BVgF3DCGNM7Q3CkP5OfA/8O9E7KmhPl/YdlNJ9gvNAYMx//z9+visjFsQ4oXMbf3zKe+1z+BpiI/2ftQeCO2IaTeKL5GYtIOvAY8E1jzMlox2GM8Rpj5gEl+H/pTIvk/oKJyIeBw8aYqmjt82wNd7KuB8YHPS4JrIs7xpj6wO1h4An8lSSeNYhIEUDg9nCM4+mXMaYh8AX0AX8gvt/beKqzUf+MRcTFn6j/Yox5PFZxABhjTgCvAxcAY0Skd+yiSH4mS4BrRWQP/iawy4BfRHH/YRvuZL0GmBw4k+oBPg08Pcz7GDIRSRORjN77wAeATaGfFXNPA18I3P8C8FQMYwmp94secD3x/d7GU52N6mccaJu9B9hqjPlZLOIQkTwRGRO4nwJcib/t/HXghkjHYIz5T2NMiTGmDP9n/5ox5rPR2v+gGGOGdQGuBrbjb3f69nCXP0wxVgAbAsvmeIsTeBB/80E3/vaym/G3oy0DdgCvAtmxjjNErH8GNgLV+L/4RbGOc4DXEPU6Gw+fMXAh/iaOamB9YLk6mnEAc4B1gRg2Ad8NrK8AVgM7gUeApCh8JkuBZ2O1/4EWvdxcKaUSwGg+waiUUglDk7VSSiUATdZKKZUANFkrpVQC0GStlFIJQJN1EBH5oogcCYwUt01E/k/Q324RkTYRyQ9a19Ln+R8VESMig74Cq29ZwyXwmn4VibJVfAp85uPC2O77InLFMOxvae9odcNNRO4VkRsG3nLk02R9uoeN/9LXJcC3RST46rajwL+EeO6NwNuBW6UiJjBaYH++CAyYrI0x3zXGvDpsQamIGhXJWkTKAkfK94rIdhH5i4hcISLvBMbsPe1yaGNMI/4O8cFX4/0R+JSIZJ9hH+n4LzK4Gf+VUP3F8mRg8KjNfQeQEpH/CaxfJiJ5gXXLRWRh4H5u4LLY3qOnx0XkxcBruD2onJsCr3M1/v90etd/JDBG7zoReVVECgLrbwkMuLRcRGpF5OtBz/l8YECmDSLy58C6PBF5TETWBJZT+1BDE1RX/yIiW0XkURFJDfxtj4jcJiJrgU+IyDwRWRX4fJ4QkbGBo9CFwF8CvxBTRGSBiLwRqHcvBV1KfuqoNVD290RkrfjHeT/t12EgtrcC26wVkcVBf84UkefEPy74b0XECjynJej5N4jIvUH7vlNEVgTqXG8cIiK/CpTzKhD8S/a7gfq2SUR+LyISWL888L6sDtT7iwLrbRH5aWD7ahH5WmD9Gd+PuBfrq3KisQBl+MfunY3/P6gq/IlX8A8H+WRguy8Cvwrcn4D/iq7kwONbgH8Fvgt8L7CuJWgfnwXuCdxfASzoJ5bswG0K/iu2cgKPDfDZwP3vBsWxHFgYuJ8L7AmKtRbIApKBvfjHuCgC9gF5gAd4J6issbw37+bfAXcEvbYVQFJgH434h6qcif/Kvtw+sT+AfyCs3vdpa6w/45GyBOqqAZYEHv8R+NfA/T3AvwdtWw1cErj/feDnZ6gzbuCzzQs8/hTwx8D9e4Ebgsr+WuD+V4C7zxBbatD3YTJQGbi/FOjAf9WfjX/kvN5yg78jNwD3Bu37Efzfxxn4h6kF+Fjg+Tb+XwcngsrKDirrz8BHgl5vb12+Gng1cP/L+Ic5dXqfH+r9iPeld6CS0WC3MWYjgIhsxj+4uhGRjfi/IL0+Jf4R+KYB/2SM6ehTzp3AehH5aZ/1N+IfAAb8A8LciP8/hb6+LiLXB+6Px1/pG/EPz/hwYP39wONneG5fy4wxTYHXtAUoxZ9slxtjjgTWPwxMCWxfAjwcOJLwALuDynrOGNMJdIrIYfzDYl4GPGKMOQpgjOkdf/kKYEbgwAb8R1Xp5r1xidXQ7DfGvBO4fz/wdaC3vj0MICJZwBhjzBuB9ffhT359TQVmAa8EPi8b/2XuZ9Jb56rwJ82+XOBXIjIP8PJevQJYbYypDcT2IP5fmY+eXsT7PGn8g31t6f2VB1wMPGiM8QIHROS1oO0vFZF/x/+fRjb+oSKeOUPsZYH7VwC/NYGhTo0xx0RkFuG/H3FlNCXrzqD7vqDHPt7/PjxsjPkn8Tc9vCwiTxtjDvX+0RhzQkQeAL7au078zSKXAbNFxOCvAEZE/s0E/vsObLcUfwW6wBjTJiLL8R8Vn0nv83p4r7mq77bBr8nLwJ/nL4GfGWOeDsRyy1mWZQHnn+E/MjU8+o4BEfy4dZBlCbDZGHNBGNv21oH+Pv//AzQAc/HXgeDPv7+Yg9eHqr9CCCKSDNyF/xfDfhG5pU95A8UevJ9w34+4MirarM+GMaYS/0+tb5zhzz8DvsR7leIG4M/GmFJjTJkxZjz+o9aL+jwvCzgeSNTT8E+l1MvivVG+PoP/RCX4f54uCNrPQN4FLhGRHPEPf/mJPvvvHerxC6c983Sv4W8bzYFT/ykBvAx8rXejwJGWGj4TRKQ3mQTXhVMCv6iO97bPAp8Deo+ym4GMwP0aIK+3PBFxRWTmWcaVBRwMHA1/Dv9BSa/zxD9yoYW/aaE35gYRmR5Yfz0DexP/r1s78Avw0sD63sR8VPznh8L5LrwCfEkCQ50G6u9wvh9Rpck6tNuAmyQwnGqvQLPAE/jbeMHf5PFEn+c+xum9Ql4EHBHZCvwYWBX0t1b8FX4T/qP07wfW/xT4soisw9/EEZLxz/BxC7ASf3t18LxytwCPiEgV/p4tA5W1GfgR8IaIbMD/nxT4f5YvDJy02QL840BlqUGpwT8hxlb85xl+0892XwB+IiLV+Cd66K0z9wK/Ff/sKzb+xHZb4DNcDyw+Q1nhuAv4QqCcabz/KH8N8Cv89W03730fvgU8i7+dOJzmhifwj/a3Bfhf/PUY4x/r+g/4z/O8FNjfQO7Gf/6mOhDzZ4x/6rbhej+iSkfdUyqOiH96rWeNMbNiHIqKM3pkrZRSCUCPrJVSKgHokbVSSiUATdZKKZUANFkrpVQC0GStlFIJQJO1UkolgP8Pm5Xhp4+V1EAAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "%matplotlib inline\n", + "\n", + "import numpy as np\n", + "from utils import plot_AMICI\n", + "\n", + "# sim, model = compile_and_simulate('gene_expression.xml', 0.1 * np.arange(10), r_max, p_max)\n", + "\n", + "plot_AMICI('gene_expression.xml', [0, 5, 10, 15], r_max, p_max)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# References\n", + "\n", + "[1] Shahrezaei, V., and Swain, P.S. (2008). \"Analytical distributions for stochastic gene expression\" PNAS 17256–17261 https://www.pnas.org/content/105/45/17256\n", + "\n", + "[2] Gillespie, D. T. (1992). \"A rigorous derivation of the chemical master equation\" Physica A: Statistical Mechanics and its Applications https://www.sciencedirect.com/science/article/abs/pii/037843719290283V\n", + "\n", + "[3] Munsky, B. and Khammash, M. (2006) \"The finite state projection algorithm for the solution of the chemical master equation\" The Journal of chemical physics https://aip.scitation.org/doi/full/10.1063/1.2145882" + ] + } + ], + "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.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/doc/examples/Finite_State_Projection/utils.py b/doc/examples/Finite_State_Projection/utils.py new file mode 100644 index 0000000..256c7ae --- /dev/null +++ b/doc/examples/Finite_State_Projection/utils.py @@ -0,0 +1,121 @@ +import amici +import amici.plotting + +import importlib +import matplotlib.pyplot as plt +import numpy as np +import os +import sys + + +def plot_AMICI(sbml_dir: str, + t: np.ndarray, + r_max: int, + p_max: int): + """ + Compiles, simulates and plots the AMICI model for the FSP example. + """ + n_t = len(t) + # compile and Simulate model + simulation, model = compile_and_simulate(sbml_dir, + t, + r_max, + p_max) + + def get_obs_idx_by_id(obs_id: str): + try: + return model.getObservableIds().index(obs_id) + except ValueError: + raise IndexError(f'No observable with id {obs_id}') + + for i in range(n_t): + + r_marginal = [simulation.y[i, get_obs_idx_by_id(f'x_r{r}')] for r in range(r_max)] + p_marginal = [simulation.y[i, get_obs_idx_by_id(f'x_p{p}')] for p in range(p_max)] + # rna + plt.subplot(n_t, 2, 2*i+1) + plt.fill_between(np.arange(r_max), 0, r_marginal, facecolor='blue', alpha=0.5) + plt.plot(r_marginal) + plt.ylabel(f't={int(t[i])}') + plt.yticks([]) + plt.xlim(0, r_max-1) + plt.ylim(0, 0.3) + + if i < n_t-1: + plt.tick_params(labelbottom=False) + else: + plt.xlabel('mRNA abundance') + + # protein + plt.subplot(n_t, 2, 2*(i+1)) + plt.fill_between(np.arange(p_max), 0, p_marginal, facecolor='blue', alpha=0.5) + plt.plot(p_marginal) + plt.yticks([ ]) + plt.xlim(0, p_max-1) + plt.ylim(0, 0.2) + if i < n_t-1: + plt.tick_params(labelbottom=False) + else: + plt.xlabel('protein abundance') + + plt.subplots_adjust(wspace=0.07, hspace=0.1) + plt.show() + + +def compile_and_simulate(sbml_dir: str, + t: np.ndarray, + r_max: int, + p_max: int): + """ + Utility function, that compiles and simulates the FSP model. + """ + model_name = 'Gene_regulation_FSP' + + # define marginal rna concentrations + observables_rna = {f'x_r{r}': get_marginal_rna_obs(r, p_max) + for r in range(r_max)} + + # define marginal protein concentrations + observables_protein = {f'x_p{p}': get_marginal_protein_obs(p, r_max) + for p in range(p_max)} + + observables = {**observables_rna, **observables_protein} + + sbml_importer = amici.sbml_import.SbmlImporter(sbml_dir) + sbml_importer.sbml2amici( + model_name=model_name, + output_dir=os.path.join('amici_models', model_name), + observables=observables) + + # Import the compiled model. + sys.path.insert(0, + os.path.abspath(os.path.join('amici_models', model_name))) + model_module = importlib.import_module(model_name) + model = model_module.getModel() + + # Simulate. + model.setTimepoints(t) + solver = model.getSolver() + simulation = amici.runAmiciSimulation(model, solver) + + return simulation, model + + +def get_marginal_rna_obs(r: int, p_max: int): + """Get the observable for a marginalized RNA abundance.""" + marginal = '' + for p in range(p_max-1): + marginal += f'x_{r}_{p} + ' + marginal += f'x_{r}_{p_max-1}' + + return {'name': f'x_r{r}', 'formula': marginal} + + +def get_marginal_protein_obs(p: int, r_max: int): + """Get the observable for a marginalized protein abundance.""" + marginal = '' + for r in range(r_max-1): + marginal += f'x_{r}_{p} + ' + marginal += f'x_{r_max-1}_{p}' + + return {'name': f'x_p{p}', 'formula': marginal} diff --git a/doc/examples/Lotka_Volterra_CLI/Lotka_Volterra_CLI.ipynb b/doc/examples/Lotka_Volterra/Lotka_Volterra_CLI/Lotka_Volterra_CLI.ipynb similarity index 100% rename from doc/examples/Lotka_Volterra_CLI/Lotka_Volterra_CLI.ipynb rename to doc/examples/Lotka_Volterra/Lotka_Volterra_CLI/Lotka_Volterra_CLI.ipynb diff --git a/doc/examples/Lotka_Volterra/Lotka_Volterra_Model_Editor/Lotka_Volterra_Model_Editor.ipynb b/doc/examples/Lotka_Volterra/Lotka_Volterra_Model_Editor/Lotka_Volterra_Model_Editor.ipynb new file mode 100644 index 0000000..787a336 --- /dev/null +++ b/doc/examples/Lotka_Volterra/Lotka_Volterra_Model_Editor/Lotka_Volterra_Model_Editor.ipynb @@ -0,0 +1,205 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Model Editor\n", + "\n", + "This notebook demonstrates functionality of the `yaml2sbml` model editor using the Lotka Volterra equations as an example. The \"Lotka-Volterra\" equations are given by \n", + "\n", + "\\begin{align*}\n", + "\\frac{d}{dt} x_1 &= \\alpha x_1 - \\beta x_1x_2, \\\\\n", + "\\frac{d}{dt} x_2 &= \\delta x_1x_2 - \\gamma x_2.\n", + "\\end{align*}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## ODE model\n", + "\n", + "We first generate a basic model, that only contains all necessary information for generating an SBML file for model simulation." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from yaml2sbml import YamlModel\n", + "\n", + "# generate model\n", + "model = YamlModel()\n", + "\n", + "# add ODEs\n", + "model.add_ode(state_id='x_1', \n", + " right_hand_side='alpha*x_1 - beta*x_1*x_2', \n", + " initial_value=2)\n", + "model.add_ode(state_id='x_2', \n", + " right_hand_side='delta*x_1*x_2 - gamma*x_2', \n", + " initial_value=2)\n", + "\n", + "# add parameters\n", + "model.add_parameter(parameter_id='alpha', nominal_value=2)\n", + "model.add_parameter(parameter_id='beta', nominal_value=4)\n", + "model.add_parameter(parameter_id='gamma', nominal_value=3)\n", + "model.add_parameter(parameter_id='delta', nominal_value=3)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`yaml2sbml` can export the `model` object either to YAML or to SBML directly, via " + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# write to YAML\n", + "model.write_to_yaml('Lotka_Volterra_basic.yaml', overwrite=True)\n", + "\n", + "# write to SBML\n", + "model.write_to_sbml('Lotka_Volterra_basic.xml', overwrite=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There are further functions to: \n", + "* get all `parameter_id`s via `model.get_parameter_ids()` \n", + "* get a parameter by its id (`model.get_parameter_by_id('alpha')`) \n", + "* delete a parameter by its id (`model.delete_parameter('alpha')`)\n", + "\n", + "Similar functions also exist for the other model components." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Parameter Estimation Problem\n", + "\n", + "Now we want to extend the current `model` to include all the necessary information for parameter estimation in PEtab. Therefore we load the model from the `.yaml` file and modify the parameters, such that it also contains all information that is going to be written into the PEtab parameter table." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "model = YamlModel.load_from_yaml('Lotka_Volterra_basic.yaml')\n", + "\n", + "# extend parameters \n", + "model.add_parameter(parameter_id='alpha',\n", + " nominal_value=2,\n", + " parameter_scale='log10',\n", + " lower_bound=0.1,\n", + " upper_bound=10,\n", + " estimate=1, \n", + " overwrite=True)\n", + "\n", + "model.add_parameter(parameter_id='beta',\n", + " nominal_value=4,\n", + " parameter_scale='log10',\n", + " lower_bound=0.1,\n", + " upper_bound=10,\n", + " estimate=1, \n", + " overwrite=True)\n", + "\n", + "model.add_parameter(parameter_id='gamma',\n", + " nominal_value=3,\n", + " parameter_scale='log10',\n", + " lower_bound=0.1,\n", + " upper_bound=10,\n", + " estimate=1,\n", + " overwrite=True)\n", + "\n", + "model.add_parameter(parameter_id='delta',\n", + " nominal_value=3,\n", + " parameter_scale='log10',\n", + " lower_bound=0.1,\n", + " upper_bound=10,\n", + " estimate=1,\n", + " overwrite=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A parameter fitting problem in PEtab allows for the specification of observables and experimental conditions:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# specify an observable:\n", + "model.add_observable(observable_id='prey_measured', \n", + " observable_formula='log10(x_1)',\n", + " noise_formula='noiseParameter1_prey_measured', \n", + " noise_distribution='normal',\n", + " observable_transformation='lin')\n", + "\n", + "# specify trivial condition\n", + "model.add_condition(condition_id='condition1', \n", + " condition_dict={})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The modified model can either be exported to YAML or PEtab via" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "# write to YAML\n", + "model.write_to_yaml('Lotka_Volterra_PEtab.yaml', overwrite=True)\n", + "\n", + "# write to PEtab\n", + "model.write_to_petab(output_dir='./Lotka_Volterra_PEtab',\n", + " model_name='Lotka_Volterra',\n", + " petab_yaml_name='Lotka_Volterra_problem', \n", + " measurement_table_name='measurement.tsv')" + ] + } + ], + "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.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/doc/examples/Lotka_Volterra_python/Lotka_Volterra.ipynb b/doc/examples/Lotka_Volterra/Lotka_Volterra_python/Lotka_Volterra.ipynb similarity index 100% rename from doc/examples/Lotka_Volterra_python/Lotka_Volterra.ipynb rename to doc/examples/Lotka_Volterra/Lotka_Volterra_python/Lotka_Volterra.ipynb diff --git a/doc/examples/Lotka_Volterra_python/Lotka_Volterra_PEtab.yml b/doc/examples/Lotka_Volterra/Lotka_Volterra_python/Lotka_Volterra_PEtab.yml similarity index 100% rename from doc/examples/Lotka_Volterra_python/Lotka_Volterra_PEtab.yml rename to doc/examples/Lotka_Volterra/Lotka_Volterra_python/Lotka_Volterra_PEtab.yml diff --git a/doc/examples/Lotka_Volterra_python/Lotka_Volterra_PEtab/measurement_table.tsv b/doc/examples/Lotka_Volterra/Lotka_Volterra_python/Lotka_Volterra_PEtab/measurement_table.tsv similarity index 100% rename from doc/examples/Lotka_Volterra_python/Lotka_Volterra_PEtab/measurement_table.tsv rename to doc/examples/Lotka_Volterra/Lotka_Volterra_python/Lotka_Volterra_PEtab/measurement_table.tsv diff --git a/doc/examples/Lotka_Volterra_python/Lotka_Volterra_basic.yml b/doc/examples/Lotka_Volterra/Lotka_Volterra_python/Lotka_Volterra_basic.yml similarity index 100% rename from doc/examples/Lotka_Volterra_python/Lotka_Volterra_basic.yml rename to doc/examples/Lotka_Volterra/Lotka_Volterra_python/Lotka_Volterra_basic.yml diff --git a/doc/examples/README.md b/doc/examples/README.md index 5cac4f6..e94b30f 100644 --- a/doc/examples/README.md +++ b/doc/examples/README.md @@ -4,17 +4,23 @@ The [`yaml2sbml`](https://github.com/yaml2sbml-dev/yaml2sbml) package translates # Scope of the Notebooks -* [Lotka Volterra](./Lotka_Volterra_python/Lotka_Volterra.ipynb) (Python) - * Introduces the input format, syntax & capabilities of `yaml2sbml` - * Showcases ODE simulation & parameter fitting using [AMICI](https://github.com/AMICI-dev/AMICI) & [pyPESTO](https://github.com/ICB-DCM/pyPESTO) -* [Lotka Volterra CLI](./Lotka_Volterra_CLI/Lotka_Volterra_CLI.ipynb) - * Shows the command line interface of `yaml2sbml` -* Sorensen (planned) - * Application example, model of glucose and insulin metabolism -* Pom1p (planned, work in progress) - * Application example, discretization of a PDE model of Pom1 gradient formation -* [Format Features](./Format_Features/Format_Features.ipynb) - * several didactic examples, that show individual features of `yaml2sbml`: - * Time-dependent right hand sides - * Step functions in the right hand side - * Function definitions +* _Lotka Volterra Notebooks_ + * [Lotka Volterra (Python)](./Lotka_Volterra/Lotka_Volterra_python/Lotka_Volterra.ipynb) + * Introduces the input format, syntax & capabilities of `yaml2sbml` + * Showcases ODE simulation & parameter fitting using [AMICI](https://github.com/AMICI-dev/AMICI) & [pyPESTO](https://github.com/ICB-DCM/pyPESTO) + * [Lotka Volterra CLI](./Lotka_Volterra/Lotka_Volterra_CLI/Lotka_Volterra_CLI.ipynb) + * Shows the command line interface of `yaml2sbml` + * [Lotka Volterra Model Editor](./Lotka_Volterra/Lotka_Volterra_Model_Editor/Lotka_Volterra_Model_Editor.ipynb) + * Shows model construction using the model editor of `yaml2sbml` +* [Sorensen](./Sorensen/yaml2sbml_Sorensen.ipynb) + * Application example, model of glucose and insulin metabolism + * Shows the extension of an existing YAML model by the model editor +* [Finite State Projection](./Finite_State_Projection/Finite_State_Projection.ipynb) + * Application example of a stochastic gene transcription model having hundreds of states. + * Generates a complex and realistic model within a few lines of Python. + +* [Format Features](./Format_Features/Format_Features.ipynb) + * Several didactic examples, that show individual features of `yaml2sbml`: + * Time-dependent right hand sides + * Step functions in the right hand side + * Function definitions diff --git a/doc/examples/Sorensen/Sorensen1985.yaml b/doc/examples/Sorensen/Sorensen1985.yaml new file mode 100644 index 0000000..7c1f6a7 --- /dev/null +++ b/doc/examples/Sorensen/Sorensen1985.yaml @@ -0,0 +1,590 @@ +time: + variable: t + +odes: + + - initialValue: GlucBV0 + rightHandSide: (GlucH - GlucBV) * QfloGB / VolGBV - VolBI / (TdifB * VolGBV) + * (GlucBV - GlucBI) + stateId: GlucBV + + - initialValue: GlucBI0 + rightHandSide: 1 / TdifB * (GlucBV - GlucBI) - GammaBGU / VolBI + stateId: GlucBI + + - initialValue: GlucH0 + rightHandSide: (QfloGB * GlucBV + QfloGL * GlucL + QfloGK * GlucK + QfloGP + * GlucPV - QfloGH * GlucH - GammaRBCU) / VolGH + stateId: GlucH + + - initialValue: GlucJ0 + rightHandSide: (GlucH - GlucJ) * QfloGJ / VolGJ - GammaJGU / VolGJ + stateId: GlucJ + + - initialValue: GlucL0 + rightHandSide: (QfloGA * GlucH + QfloGJ * GlucJ - QfloGL * GlucL + GammaHGP + - GammaHGU) / VolGL + stateId: GlucL + + - initialValue: GlucK0 + rightHandSide: (GlucH - GlucK) * QfloGK / VolGK - GammaKGE / VolGK + stateId: GlucK + + - initialValue: GlucPV0 + rightHandSide: QfloGP / VolGPV * (GlucH - GlucPV) - VolPI / (TdifGP * VolGPV) + * (GlucPV - GlucPI) + stateId: GlucPV + + - initialValue: GlucPI0 + rightHandSide: (GlucPV - GlucPI) / TdifGP - GammaPGU / VolPI + stateId: GlucPI + + - initialValue: MIHGP0 + rightHandSide: (MIHGPinf - MIHGP) / tauInsu + stateId: MIHGP + + - initialValue: Func20 + rightHandSide: ((MC0HGP - 1.0) / 2.0 - Fun2) / tauCgon + stateId: Fun2 + + - initialValue: MIHGU0 + rightHandSide: (MIHGUinf - MIHGU) / tauInsu + stateId: MIHGU + + - initialValue: InsuB0 + rightHandSide: QfloIB / VolIB * (InsuH - InsuB) + stateId: InsuB + + - initialValue: InsuH0 + rightHandSide: (QfloIB * InsuB + QfloIL * InsuL + QfloIK * InsuK + QfloIP * + InsuPV - QfloIH * InsuH) / VolIH + stateId: InsuH + + - initialValue: InsuJ0 + rightHandSide: QfloIJ / VolIJ * (InsuH - InsuJ) + stateId: InsuJ + + - initialValue: InsuL0 + rightHandSide: (QfloIA * InsuH + QfloIJ * InsuJ - QfloIL * InsuL + GammaPIR + - GammaLIC) / VolIL + stateId: InsuL + + - initialValue: InsuK0 + rightHandSide: (QfloIK / VolIK) * (InsuH - InsuK) - GammaKIC / VolIK + stateId: InsuK + + - initialValue: InsuPV0 + rightHandSide: (QfloIP/VolIPV) * (InsuH - InsuPV) - VolPI / (VolIPV * TdifIP) + * (InsuPV - InsuPI) + stateId: InsuPV + + - initialValue: InsuPI0 + rightHandSide: (1 / TdifIP) * (InsuPV - InsuPI) - GammaPIC / VolPI + stateId: InsuPI + + - initialValue: Potn0 + rightHandSide: KappaPotnPtgt * (Ptgt - Potn) + stateId: Potn + + - initialValue: Pinh0 + rightHandSide: KappaPinhPrp * (Pprp - Pinh) + stateId: Pinh + + - initialValue: InitialRinsu0 + rightHandSide: KappaRinsu * (Rinsu0 - Rinsu) + KappaRinsuPotn * Potn - Secr + stateId: Rinsu + + - initialValue: Cgon0 + rightHandSide: (GammaPCR - GammaPCC) / VolC + stateId: Cgon + +parameters: + + - nominalValue: 0.59 + parameterId: QfloGB + + - nominalValue: 0.35 + parameterId: VolGBV + + - nominalValue: 0.45 + parameterId: VolBI + + - nominalValue: 2.1 + parameterId: TdifB + + - nominalValue: 5.07333 + parameterId: GlucH0 + + - nominalValue: 0.388889 + parameterId: GammaBGU + + - nominalValue: 1.26 + parameterId: QfloGL + + - nominalValue: 1.01 + parameterId: QfloGK + + - nominalValue: 1.51 + parameterId: QfloGP + + - nominalValue: 4.37 + parameterId: QfloGH + + - nominalValue: 0.0555556 + parameterId: GammaRBCU + + - nominalValue: 1.38 + parameterId: VolGH + + - nominalValue: 1.01 + parameterId: QfloGJ + + - nominalValue: 1.12 + parameterId: VolGJ + + - nominalValue: 0.111111 + parameterId: GammaJGU + + - nominalValue: 0.25 + parameterId: QfloGA + + - nominalValue: 2.51 + parameterId: VolGL + + - nominalValue: 0.66 + parameterId: VolGK + + - nominalValue: 1.04 + parameterId: VolGPV + + - nominalValue: 6.74 + parameterId: VolPI + + - nominalValue: 5.0 + parameterId: TdifGP + + - nominalValue: 0.194444 + parameterId: GammaBPGU + + - nominalValue: 7.03 + parameterId: beta0PGU + + - nominalValue: 6.52 + parameterId: beta1PGU + + - nominalValue: 0.338 + parameterId: beta2PGU + + - nominalValue: 5.82 + parameterId: beta3PGU + + - nominalValue: 2.7 + parameterId: beta0HGP + + - nominalValue: 0.388852 + parameterId: beta1HGP + + - nominalValue: 65.0 + parameterId: tauCgon + + - nominalValue: 1.21 + parameterId: beta2HGP + + - nominalValue: 1.14 + parameterId: beta3HGP + + - nominalValue: 1.66 + parameterId: beta4HGP + + - nominalValue: 0.887748 + parameterId: beta5HGP + + - nominalValue: 25.0 + parameterId: tauInsu + + - nominalValue: 1.42 + parameterId: beta6HGP + + - nominalValue: 1.41 + parameterId: beta7HGP + + - nominalValue: 0.62 + parameterId: beta8HGP + + - nominalValue: 0.504543 + parameterId: beta9HGP + + - nominalValue: 0.861111 + parameterId: GammaHGP0 + + - nominalValue: 2.0 + parameterId: beta0HGU + + - nominalValue: 0.549306 + parameterId: beta1HGU + + - nominalValue: 5.66 + parameterId: beta2HGU + + - nominalValue: 5.66 + parameterId: beta3HGU + + - nominalValue: 2.44 + parameterId: beta4HGU + + - nominalValue: 1.4783 + parameterId: beta5HGU + + - nominalValue: 0.111111 + parameterId: GammaHGU0 + + - nominalValue: 0.394444 + parameterId: beta0KGE + + - nominalValue: 0.394444 + parameterId: beta1KGE + + - nominalValue: 0.198 + parameterId: beta2KGE + + - nominalValue: 25.5556 + parameterId: beta3KGE + + - nominalValue: 1.834 + parameterId: beta4KGE + + - nominalValue: 0.0872 + parameterId: beta5KGE + + - nominalValue: 0.45 + parameterId: QfloIB + + - nominalValue: 0.26 + parameterId: VolIB + + - nominalValue: 0.99 + parameterId: VolIH + + - nominalValue: 0.9 + parameterId: QfloIL + + - nominalValue: 0.72 + parameterId: QfloIK + + - nominalValue: 1.05 + parameterId: QfloIP + + - nominalValue: 3.12 + parameterId: QfloIH + + - nominalValue: 0.94 + parameterId: VolIJ + + - nominalValue: 0.72 + parameterId: QfloIJ + + - nominalValue: 1.14 + parameterId: VolIL + + - nominalValue: 0.18 + parameterId: QfloIA + + - nominalValue: 0.4 + parameterId: FracLIC + + - nominalValue: 0.3 + parameterId: FracKIC + + - nominalValue: 0.51 + parameterId: VolIK + + - nominalValue: 0.74 + parameterId: VolIPV + + - nominalValue: 20.0 + parameterId: TdifIP + + - nominalValue: 0.15 + parameterId: FracPIC + + - nominalValue: 3.27 + parameterId: beta1PIR + + - nominalValue: 7.33333 + parameterId: beta2PIR + + - nominalValue: 2.879 + parameterId: beta3PIR + + - nominalValue: 3.02 + parameterId: beta4PIR + + - nominalValue: 1.11 + parameterId: beta5PIR + + - nominalValue: 0.00794 + parameterId: KappaRinsu + + - nominalValue: 44310.0 + parameterId: Rinsu0 + + - nominalValue: 4025.0 + parameterId: KappaRinsuPotn + + - nominalValue: 0.0482 + parameterId: KappaPotnPtgt + + - nominalValue: 0.931 + parameterId: KappaPinhPrp + + - nominalValue: 0.00747 + parameterId: EMME1 + + - nominalValue: 0.0958 + parameterId: EMME2 + + - nominalValue: 91.0 + parameterId: InsuPV0 + + - nominalValue: 11.48 + parameterId: Cgon0 + + - nominalValue: 0.91 + parameterId: GammaMCC + + - nominalValue: 11.31 + parameterId: VolC + + - nominalValue: 2.93 + parameterId: beta0PCR + + - nominalValue: 2.1 + parameterId: beta1PCR + + - nominalValue: 4.18 + parameterId: beta2PCR + + - nominalValue: 0.621325 + parameterId: beta3PCR + + - nominalValue: 1.31 + parameterId: beta4PCR + + - nominalValue: 0.61 + parameterId: beta5PCR + + - nominalValue: 1.06 + parameterId: beta6PCR + + - nominalValue: 0.471419 + parameterId: beta7PCR + + - nominalValue: 0.0 + parameterId: Func20 + +assignments: + + - assignmentId: InsuH0 + formula: InsuPV0/(1-FracPIC) + + - assignmentId: InsuK0 + formula: InsuH0*(1-FracKIC) + + - assignmentId: InsuB0 + formula: InsuH0 + + - assignmentId: InsuJ0 + formula: InsuH0 + + - assignmentId: InsuPI0 + formula: InsuPV0-((QfloIP*TdifIP/VolPI)*(InsuH0-InsuPV0)) + + - assignmentId: InsuL0 + formula: 1/QfloIL*(QfloIH*InsuH0-QfloIB*InsuB0-QfloIK*InsuK0-QfloIP*InsuPV0) + + - assignmentId: GammaBPIR + formula: QfloIL/(1-FracLIC)*InsuL0 - QfloIJ*InsuJ0-QfloIA*InsuH0 + + - assignmentId: GammaPIC0 + formula: InsuPI0/(((1-FracPIC)/FracPIC)*(1/QfloIP)-TdifIP/VolPI) + + - assignmentId: Pprp0 + formula: (GlucH0^beta1PIR) /((beta2PIR^beta1PIR)+beta3PIR*(GlucH0^beta4PIR)) + + - assignmentId: Ptgt0 + formula: (Pprp0^beta5PIR) + + - assignmentId: Pinh0 + formula: Pprp0 + + - assignmentId: Potn0 + formula: Ptgt0 + + - assignmentId: InitialRinsu0 + formula: ((KappaRinsu*Rinsu0)+ KappaRinsuPotn * Potn0)/(KappaRinsu+EMME1* Potn0) + + - assignmentId: Secr0 + formula: EMME1*Ptgt0*InitialRinsu0 + + - assignmentId: GlucPV0 + formula: GlucH0 - GammaBPGU/QfloGP + + - assignmentId: GlucK0 + formula: GlucH0 + + - assignmentId: GlucBV0 + formula: GlucH0 - GammaBGU/QfloGB + + - assignmentId: GlucJ0 + formula: GlucH0-GammaJGU/QfloGJ + + - assignmentId: GlucL0 + formula: (QfloGA*GlucH0+QfloGJ*GlucJ0+GammaHGP0-GammaHGU0)/QfloGL + + - assignmentId: GlucBI0 + formula: GlucBV0-(GammaBGU*TdifB)/VolBI + + - assignmentId: GlucPI0 + formula: GlucPV0-GammaBPGU*TdifGP/VolPI + + - assignmentId: MIPGU0 + formula: beta0PGU+beta1PGU*tanh(beta2PGU*(1-beta3PGU)) + + - assignmentId: MCHGP0 + formula: beta0HGP * tanh(beta1HGP * 1) - Func20 + + - assignmentId: MC0HGP0 + formula: beta0HGP * tanh(beta1HGP * 1) + + - assignmentId: MIHGP0 + formula: beta2HGP - beta3HGP * tanh(beta4HGP * (1-beta5HGP)) + + - assignmentId: MIHGPinf0 + formula: MIHGP0 + + - assignmentId: MGHGP0 + formula: beta6HGP-beta7HGP*tanh(beta8HGP*(1-beta9HGP)) + + - assignmentId: MIHGU0 + formula: beta0HGU * tanh(beta1HGU) + + - assignmentId: MIHGUinf0 + formula: MIHGU0 + + - assignmentId: MGHGU0 + formula: beta2HGU+beta3HGU*tanh(beta4HGU*(1-beta5HGU)) + + - assignmentId: GammaKGE0 + formula: piecewise(beta0KGE+beta1KGE*tanh(beta2KGE*(GlucK0-beta3KGE)), GlucK0 Pinh, + EMME1 * Ptgt * Rinsu) + + - assignmentId: SecrN + formula: Secr / Secr0 + + - assignmentId: Pprp + formula: (GlucH^beta1PIR)/((beta2PIR^beta1PIR)+beta3PIR*(GlucH^beta4PIR)) + + - assignmentId: Ptgt + formula: Pprp^beta5PIR + + - assignmentId: CgonN + formula: Cgon / Cgon0 + + - assignmentId: GammaPCC + formula: GammaMCC * Cgon + + - assignmentId: GammaPCR + formula: GammaBPCR * MGPCR * MIPCR + + - assignmentId: MGPCR + formula: beta0PCR - beta1PCR * tanh(beta2PCR * (GlucNH - beta3PCR)) + + - assignmentId: MIPCR + formula: beta4PCR - beta5PCR * tanh(beta6PCR * (InsuNH - beta7PCR)) diff --git a/doc/examples/Sorensen/utils.py b/doc/examples/Sorensen/utils.py new file mode 100644 index 0000000..c851be9 --- /dev/null +++ b/doc/examples/Sorensen/utils.py @@ -0,0 +1,100 @@ +import amici +import amici.plotting + +import importlib +import matplotlib.pyplot as plt +import numpy as np +import os +import sys + + +def simulate_and_plot_sorensen(sbml_dir: str): + """Plots and simulates the Sorensen model. + + This function is highly specific to the Sorensen model. + + Parameters: + ----------- + sbml_dir: + path to the SBML model. + """ + model = compile_model(sbml_dir) + + # Simulate. + model.setTimepoints(np.linspace(0, 300, 3001)) + solver = model.getSolver() + simulation = amici.runAmiciSimulation(model, solver) + + # Reproduce the figure. + fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5)) + + amici.plotting.plotObservableTrajectories( + simulation, + observable_indices=[0], + model=model, + ax=ax1, + ) + ax1.set_xlabel('Time [min]') + ax1.set_ylabel('Blood glucose concentration [mg/dL]') + ax1.set_title('Blood Glucose Concentration') + ax1.get_legend().remove() + + amici.plotting.plotObservableTrajectories( + simulation, + observable_indices=[1], + model=model, + ax=ax2, + ) + ax2.set_xlabel('Time [min]') + ax2.set_ylabel('Plasma insulin concentration [mU/L]') + ax2.set_title('Plasma Insulin Concentration') + ax2.get_legend().remove() + + plt.show() + + +def compile_model(sbml_dir: str): + """Compiles the Sorensen model and specifies the observables. + + Parameters: + ----------- + sbml_dir: + path to the SBML model. + """ + # observables are defined, as specifying a whole PEtab problem including + # observables would mean a huge overhead. + + # Define the variables that will be plotted to reproduce the figure from + # the original thesis. + observables = { + 'blood_glucose': { + 'name': 'Blood glucose concentration [mg/dL]', + # GlucPV: model state for glucose in the peripheral vascular space. + # Multiplied by 18.0156 to convert [mmol/L] to [mg/dL]. + # Multiplied by 0.84 as described in Panunzi et al. 2020 [2]. + 'formula': 'GlucPV * 18.0156 * 0.84', + }, + 'plasma_insulin': { + 'name': 'Plasma insulin concentration [mU/L]', + # InsuPV: model state for insulin in the peripheral vascular space. + # Multiplied by 0.144 to convert [pM] to [mU/L]. + 'formula': 'InsuPV * 0.144', + } + } + + # Compile the model with AMICI. + sbml_importer = amici.sbml_import.SbmlImporter(sbml_dir) + sbml_importer.sbml2amici( + model_name='Fig71_Sorensen1985', + output_dir=os.path.join('amici_models', sbml_dir), + observables=observables, + ) + + # Import the compiled model. + sys.path.insert( + 0, + os.path.abspath(os.path.join('amici_models', 'Fig71_Sorensen1985')) + ) + model_module = importlib.import_module('Fig71_Sorensen1985') + + return model_module.getModel() diff --git a/doc/examples/Sorensen/yaml2sbml_Sorensen.ipynb b/doc/examples/Sorensen/yaml2sbml_Sorensen.ipynb new file mode 100644 index 0000000..d9426a3 --- /dev/null +++ b/doc/examples/Sorensen/yaml2sbml_Sorensen.ipynb @@ -0,0 +1,197 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Glucose Insulin Metabolism Model by Sorensen (1985).\n", + "\n", + "In this notebook, the model of glucose and insulin metabolism, from the thesis of Thomas\n", + "J. Sorensen, 1985 [[1]](#References), is used as an example of how the `yaml2sbml.YamlModel` class can be used to easily extend a pre-existing `yaml2sbml`\n", + "model.\n", + "\n", + "Specifically, the model is edited to reproduce a figure from the original\n", + "thesis (Fig 71) [[1]](#References). The implementation that is loaded here is based on the implementation\n", + "provided by Panunzi et al., 2020 [[2]](#References). \n", + "\n", + "In the end, the model is exported to [SBML](http://sbml.org/) and simulated via [AMICI](https://github.com/AMICI-dev/AMICI)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Extend the current YAML model\n", + "\n", + "The Sorensen model has already been encoded in the YAML format. Within this notebook, the model is loaded, extended in order to model an intravenous infusion administration of glucose, using the `yaml2sbml` Model editor. The extended model is then plotted, to reproduce a figure from Sorensens PhD thesis.\n", + "\n", + "### Load the model" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import yaml2sbml\n", + "from yaml2sbml import YamlModel\n", + "\n", + "yaml_model = YamlModel.load_from_yaml('Sorensen1985.yaml')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Define the glucose infusion term\n", + "\n", + "An intravenous infusion administration of glucose is added to the model. After 17 minutes of initialization, glucose is infused at a rate of 64.81 mM/min for a three-minute period.\n", + "\n", + "In the YAML model, this is represented by a novel term `IVG` (\"intravenous glucose\"), that is set by a step function in an assignment rule." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# describe the intravenous glucose infusion\n", + "yaml_model.add_assignment('IVG', formula='piecewise(64.81, 17 <= t && t < 20, 0)')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Add the glucose infustion term to the model\n", + "\n", + "`IVG` is now added to the ODE of the glucose concentration of heart and lung space `GlucH`. Therefore, the current ODE is overwritten." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# Get the current ODE.\n", + "gluch_ode_old = yaml_model.get_ode_by_id('GlucH')\n", + "gluch_rhs_old = gluch_ode_old['rightHandSide']\n", + "\n", + "# Modify the ODE. `IVG` is divided by the volume term `VolGH` to\n", + "# model concentration.\n", + "gluch_rhs_new = gluch_rhs_old + ' + (IVG / VolGH)'\n", + "\n", + "# Set the new ODE.\n", + "yaml_model.add_ode(state_id='GlucH', \n", + " right_hand_side=gluch_rhs_new, \n", + " initial_value=gluch_ode_old['initialValue'], \n", + " overwrite=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Export to SBML" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# Write to SBML\n", + "yaml_model.write_to_sbml('Fig71_Sorensen1985.xml', overwrite=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Simulation in AMICI\n", + "\n", + "The model is now setup to reproduce to the figure. The utility function `simulate_and_plot_sorensen` simulates and plots the Sorensen model in [AMICI](https://github.com/AMICI-dev/AMICI)." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3sAAAFNCAYAAAC5cXZ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABpOklEQVR4nO3dd5hcZdnH8e+9szXZTd+E9ARSSAgQIHSkd2kqIlbkRbGDXbBiQbEgiqCCggJSRER6ERGCSDOBkAaBEAghpLdN3TJzv3+cM7uTZXZ3dndmz8zs73NxrplT5sx9ZsI8e5+nmbsjIiIiIiIixaUk6gBEREREREQk+5TsiYiIiIiIFCEleyIiIiIiIkVIyZ6IiIiIiEgRUrInIiIiIiJShJTsiYiIiIiIFCElexIZM/uzmf0oB+e9xMz+0s1zuJlNyFZMkntmtsDMjow6DhGR7jCzx83sE1HHESUzO9LM3kpZ1+97N5nZu8xsUdRxSM9Tsic5Y2ZvmNl2M9tiZhvM7H4zG50HcdWY2S/D+Laa2ZtmdoeZHRh1bF1hZgeY2QNmttHM1pvZc2Z2btRxtSX83I/t5jnecaPA3fdw98e7FZyISA9oVT6uCn/TqqOOKxNm9nEze7In37M7v+8WuMDM5odl/ltm9jcz2zPLYWZF60S3G+fZ6aa1u//H3Sd397xSeJTsSa6d6u7VwHBgFfCbKIMxswrg38CewClAP2AKcBtwUoShdYmZHUxwPTOBCcBg4DMU4LUkmVlp1DGIiPSAZPm4LzAD+HbE8RSrXwMXAhcAg4BJwF3AuyOMqVtUTkpnKNmTHuHuO4A7gKltHWNmnzSzxWHt1D1mNiJl3yFm9j8z2xQ+HpKyb7yZzTSzzWb2CDCknVA+CowCznD3+e4ed/et7n6Hu1/SRlw7NalpfVfTzPYws0fCuFeZ2TfD7RVm9iszeztcfhUmm5jZEDO7L6U27j9mVhLuG2FmfzezNWb2upld0M71/By4wd1/6u5rPTDb3c/K8HN1M/u0mb0axnK1mVmr174UfrYLzWzfjmIMm9HebmY3hq9bYGYzwn03AWOAe8M72l83s3FhHOeZ2ZsEySvhndeV4Xf+hJntEW4/H/gw8PXwHPeG25trDDv47I8M7+x+xcxWm9kKy+OaUBEpbu6+HHgQmNZ6n5ntZmb/NrN1ZrbWzG42swEp+79hZsvD39pFZnZMuP2S8Df0L+G+eWY2ycwuDn/3lpnZ8SnnOTflt36JmX0q0/jD396vmtnc8Pf6r2ZWGe5rr6zbqebJ2una0er3vc0yJs3rJgKfAz7o7v9293p33+buN7v7ZeEx/cNzrTGzpWb27ZQYP25mT5rZLyxoofS6mZ2Ucv5BZvansJzZYGZ3pew7xczmhNf+lJnt1dFnZmZ9Cf4tjAjLty1heXuJBS2Q/mJmdcDHLWjV83R4/hVmdpWZlYfnfyJ8qxfDc3zA3tk0dooFf99sDD/D01p9F1db0CJrs5k9a2a7ZfDPQfKQkj3pEWbWB/gA8Ewb+48GfgKcRVALuJSgtg0zGwTcD1xJUHP1S+B+MxscvvwWYDZBkvdD4Jx2QjkWeNjdt3bzkpJx1wD/Ah4CRhDUrj0a7v4WcBAwHdgbOICWO7dfAd4CaoFhwDcBDwuYe4EXgZHAMcAXzeyENO/dBziYIIluK742P9cUpwD7A3uFx50Qvvb9wCXAxwhqQE8D1mUY42nh+wwA7gGuAnD3jwJvEt7RdvefpbzmCIJa1uR5HgQmAkOB54Gbw3NcGz7/WXiOU9NcenufPcAuQP8w/vOAq81sYJrziIjklAXdG04GXki3m+A3fATB7+Nogt9lzGwy8Hlgf3evIfjtfCPltacCNwEDw3M/TPB330jgB8A1KceupqW1y7nAFRbe3MvQWcCJwHiCsuTj4fa0ZV0nztuWtGVMGscAb7n7c+2c6zcE5cGuBOXQxwg+g6QDgUUEf2P8DLjOrPmm6E1AH2APgrLqCgAz2we4HvgUwd8t1wD3WHjTMfSOzyz82+Qk4O2wfKt297fD408nKO8HEJSBceBLYVwHh9f6WQB3Pzx8zd7hOf6aesFmVkZQjv8zjPsLwM3hv6mks4HvE/z7WQxc2s5nKHlMyZ7k2l1mthHYBBxHUBOVzoeB6939eXevBy4GDjazcQRNLV5195vcvcndbwVeBk41szEEicp3wjt2TxD8gLVlCLAyuWJm08O7WnXWtY7LpwAr3f1yd9/h7pvd/dmUa/qBu6929zUEP5ofDfc1EiRfY929MWxL7+G11Lr7D9y9wd2XAH8g+NFtbSDB/8Mr2omvvc816TJ33+jubwKPESRIAJ8gSKj+F9YYLnb3pRnG+KS7P+DucYLCcO92P8XAJWEt63YAd78+/DzrCf642dvM+mdwnuR1t/XZQ/D5/yD87B8AtgDqyyAiPSlZPj5J0BT/x60PCH93HwnLtzUENzuPCHfHgQpgqpmVufsb7v5aysv/4+4Pu3sT8DeChOsyd28kSJTGWVhL6O73u/tr4W/9TIIk4F2duJYr3f1td19PUAZPD7e3VdZ1V6ZlzGDaKSPNLEZQdl0cljdvAJezc3mx1N3/EL7XDeH1DDOz4QSJ2afdfUN4fTPD15wPXOPuz4YtiG4A6gluQia19Zm15Wl3v8vdE+6+PWzF80z4d9EbBAnlER2cI+kgoJrg30ODu/8buA/4YMox/3D358J/PzdnEJ/kKSV7kmtnuPsAoJLgDuRMM9slzXEjCGqdAHD3LcA6gjuQO+0LLU3Zt6FVTV3rY1OtI/ihTr7PnDC+9xIUmp01GnitjX2t414aboMg6V0M/DNsMnNRuH0sQfONjcmF4E7osDTn3wAkUq+noxhafa5JK1OebyMoANq7tkxibH3OSuu4j8Gy5BMzi5nZZWb2Wthk5Y1wV3tNdFO199kDrAsLsNQYC2JwBBEpGme4+wB3H+vun03e6EplZsPM7DYLmmrWAX8h/B1098XAFwluhq0Oj0v9nVuV8nw7sDZMWJLrEP7umdlJZvaMBU0tNxLUNGb6ewttlyNtlXXdlWkZs1OZn8YQoIx3lhdpy0h33xY+rSYoI9e7+4Y05x0LfKVVOTmancuhtj6ztixLXbGgWe59FnR3qCO4WdCZMnKZuydStrV53RnGJ3lKyZ70iPDO1p0EdyIPS3PI2wQ/jgCE7dYHA8tb7wuNCfetAAaGx6fua8ujwPGtju/IVoJmGkmpyeoygqYf6bSOe0y4jfAO4lfcfVeC5ihftqCvxTLg9fAPgORS4+4ntz55WOg8Dbyvndjb+1w7sgxI10Y/4xjb0NZd3dTtHyJosnIsQfOaceF2S3NsOm1+9iIiBeTHBL93e7p7P+AjtPwO4u63uPthBL93Dvy0s28QNi38O/ALYFh4A/SB1PfpqnbKOggSiLbK1mx5FBhlbfTpA9YS1D62Li8yLSMHWUofylb7Lm1VTvbxoGVSRzIpIwF+R9DKaWL4b+ObZP6dvQ2MDrtlJGV63VJglOxJj7DA6QRND19Kc8itwLlhs8oKggLu2bBpwgPAJDP7kJmVmtkHCAZ6uS9sVjgL+L6ZlZvZYQT9FNpyI0GC+A8zmxbWIFUSjITWljnAe82sjwWdyc9L2XcfMNzMvmjBoCA11jKFw63At82s1syGAN8luCub7Lg9IWz3v4kgCU4AzwGbLeh0XxXGN83M9m8jtq8TdNT+WrIPo5ntbWbJfnntfa4d+SPwVTPbL/z+JpjZ2C7E2Noq2k6Qk2oImrysI/hjoHXzpo7O0eZnLyJSQGoImplvMrORwNeSO8xsspkdHf627yCorUukP027yglatqwBmiwYgOT49l+SmXbKOgjK1g+FZciJZN4EMWPu/irwW+BWCwYoKbdgIJSzzeyisKbzduDSsPweC3yZDMoLd19B0Lf8t2Y20MzKzCzZV+4PwKfN7MCw/OxrZu+2oJ9/R1YBg63jbgs1QB2wxcx2JxiJu/V52ionnyVItr8exn0kwd9Orfv0SxFQsie5dq+ZbSH4QboUOMfdF7Q+yN3/BXyH4O7iCoIapbPDfesI+sZ9heCP/68Dp7j72vDlHyLoQL0e+B5BQpeWB6OCHgUsJBj0pY6g4/X+BJ2l07kCaCD44byBcKCQ8HybCfoinkrQ5OHV8PwAPyJIROcC8wgGGUmONDaRYGCXLQS1c79198fCgucUgrbxrxPcdfwjQe1Wuut5Cjg6XJaY2XrgWoIEud3PtSPu/jeC7+wWYDPBUNWDOhtjGj8hSMQ2mtlX2zjmRoImJcsJvqvWA/tcR9BPZaOljH6Wor3PXkSkUHyfYGqGTQRl1p0p+yqAywh+g1cSDLRxcWffICzHLiBIejYQlKn3dCvqFmnLunDfhQRl50aCftZ3Zek9W7uAYACXq8P3eg14Dy39+79A0IJnCUH/yVsIBlfJxEcJagZfJhjk5osA7j4L+GT4vhsImrJ+PJMTuvvLBDcsl4Rl3Ig2Dv0qwXe1mSC5/Gur/ZcAN4Tn2OnvG3dvIPjsTyL49/Nb4GPhe0uRsez0kxUREREREZF8opo9ERERERGRIqRkT0REREREpAgp2RMRERERESlCSvZERERERESKkJI9ERERERGRIlQadQDdMWTIEB83blzUYYiISA+YPXv2WnevjTqOQqEyUkSkd2ivfCzoZG/cuHHMmjUr6jBERKQHmNnSqGMoJCojRUR6h/bKRzXjFBERERERKUJK9kRERERERIqQkj0REREREZEipGRPRERERESkCCnZExERERERKUJK9kRERERERIqQkj0REREREZEilLNkz8wqzew5M3vRzBaY2ffD7X82s9fNbE64TA+3m5ldaWaLzWyume2bq9hERERERESKXS4nVa8Hjnb3LWZWBjxpZg+G+77m7ne0Ov4kYGK4HAj8LnwUERERERGRTspZzZ4HtoSrZeHi7bzkdODG8HXPAAPMbHiu4usJz72+nu0N8ajDEBGRPGNm15vZajObn7JtkJk9Ymavho8Dw+1F2fJl5aYdPLtkXdRhiIgUtZz22TOzmJnNAVYDj7j7s+GuS8MC6wozqwi3jQSWpbz8rXBbQVq+cTtnXfM0F905N+pQREQk//wZOLHVtouAR919IvBouA47t3w5n6DlS8H77M2z+cC1z7BpW2PUoYiIFK2cJnvuHnf36cAo4AAzmwZcDOwO7A8MAr7RmXOa2flmNsvMZq1ZsybbIWfN5h1B4bXw7bqIIxERkXzj7k8A61ttPh24IXx+A3BGyvaiavkC8PybGwFYuELlpIhIrvTIaJzuvhF4DDjR3VeEBVY98CfggPCw5cDolJeNCre1Pte17j7D3WfU1tbmOPLuM4s6AhERKRDD3H1F+HwlMCx8XlQtX1pbVbcj6hBERIpWLkfjrDWzAeHzKuA44OXk3UgzM4K7lsn+CvcAHwv7JhwEbEop9ApOIhE8Gsr2RESkc9zdab+fe1qF0vollZI9EZHcyeVonMOBG8wsRpBU3u7u95nZv82sFjBgDvDp8PgHgJOBxcA24NwcxpZzDfEg21PNnoiIZGiVmQ139xXhjdHV4faMWr5A0PoFuBZgxowZnU4We0pjWEYCrKqrjzASEZHilrNkz93nAvuk2X50G8c78LlcxdPT6huDUThN2Z6IiGTmHuAc4LLw8e6U7Z83s9sIpiQq6JYvAJt3NDU/37itIcJIRESKWy5r9nq1+qawZi/iOEREJP+Y2a3AkcAQM3sL+B5Bkne7mZ0HLAXOCg8vqpYv0DKIGcCm7RqNU0QkV5Ts5Ugy2RMREWnN3T/Yxq5j0hxbVC1fYOeaPSV7IiK50yOjcfZGyf4IedthQkREJCL1TUFXhz7lMSV7IiI5pGQvR1I7n4uIiEiL+sagjBxaU6FkT0Qkh5Ts5UhTXHV6IiIi6SS7OgytqVSyJyKSQ0r2cqQpoZo9ERGRdJLNOGv7VVDflGBHOIK1iIhkl5K9HGkMa/aCfvUiIiKS1FKzVwFAnWr3RERyQslejjSpz56IiEhaLX32KgGNyCkikitK9nKkKaEaPRERkXSSzTiHVJcDUJcyFYOIiGSPkr0cUbInIiKSXrIZZ22yGecO1eyJiOSCkr0cUTNOERGR9JLJ3pDqINnbrJo9EZGcULKXIy0DtEQciIiISJ6pb0w249QALSIiuaRkL0eSUy9oCgYREZGd1ccTVJSW0K+qFFDNnohIrijZy5HkpOqNmlxdRERkJ/WNQbJXVRajtMTUZ09EJEeU7OVIMslT3z0REZGd1TclqCiLYWbUVJayWcmeiEhOKNnLkWTzzUaNyikiIrKT+qY4FaXBnyD9qsqo265mnCIiuaBkL0eSUy+oZk9ERGRn9U0JysNkTzV7IiK5o2QvR5JJXpP67ImIiOykKZ6grCSs2ass06TqIiI5omQvR5JJXoNq9kRERHbSFHdKYwYEyZ5q9kREckPJXo4k++o1qc+eiIjIThriCcpiLc041WdPRCQ3lOzlSLIZZzzhuGZWFxERadYUd8qSNXtVqtkTEckVJXs5kjq/nubaExERadGUSFBa0lKzt7UhrgHNRERyQMlejiSnXmj9XEREpLdrbNVnD2BLvZpyiohkm5K9HGlSzZ6IiEhaTYmd++wB6rcnIpIDSvZyJLU2r1FNU0RERJo1xZ3SkpY+ewB16rcnIpJ1SvZyJLVmT3PtiYiItGiIJygrbVWzp2RPRCTrlOzlSGMitRmnavZERESSmuJOWcnOffY2a2J1EZGsK21rh5l9OYPXb3X3a7IYT9FIHVVMc+2JiIi0aIonKA377CWTvbrtqtkTEcm29mr2vgZUAzXtLF/JdYCFaudmnKrZExERSWpMpM6zF9x3Vs2eiEj2tVmzB9zk7j9o78Vm1jfL8RSNxkSCspjRGHeNxikiUoTUAqbrmuIt8+xVV6jPnohIrrRZs+fuX29rn5m9r6NjerumuFNZFgPUZ09EpEipBUwXNca9eeqF0lgJfctjqtkTEcmB9mr22nMF8PdsBlJsmuIJ+oSFlyZVFxEpSmoB00WN8URzM06Amsoy9dkTEcmBro7GaR0f0rs1JlJr9tSMU0Sk2KgFTNc1JZzSlGSvX1WpavZERHKgq8mespcOxBNOVZjsaZ49EZFe54qoA8hX7k484c199iCs2VOfPRGRrGtv6oV5pE/qDBjW0YnNrBJ4AqgI3+cOd/+emY0HbgMGA7OBj7p7g5lVADcC+wHrgA+4+xudu5z80RhPtNTsqRmniEhvoxYwbUi2dkltxtmvspS1WxqiCklEpGi112fvlG6eux442t23mFkZ8KSZPQh8GbjC3W8zs98D5wG/Cx83uPsEMzsb+CnwgW7GEJmmuGr2RER6Mf3wtyHZjz05zx4ENXtL1m6NKiQRkaLVZrLn7ku7c2J3d2BLuFoWLg4cDXwo3H4DcAlBsnd6+BzgDuAqM7PwPAWnKZGgqlyjcYqIFKvutoDprRqbkjV7Lcme+uyJiORGe804N9POnUl379fRyc0sRtBUcwJwNfAasNHdk7/obwEjw+cjgWXhuZvMbBNBU8+1HV9GfnEP5tZTsiciUtS62wKmV0p2bdi5GWcwGqe7Y6YWsCIi2dJezV4NgJn9EFgB3ERwt/LDwPBMTu7ucWC6mQ0A/gHs3s14MbPzgfMBxowZ093T5UQ8EeTIasYpIlLUrgUeAh5095ejDqZQJMvE1gO0NCWcHY0trWJERKT7MhmN8zR3/627b3b3OndPNrnMmLtvBB4DDgYGmFkyyRwFLA+fLwdGA4T7+xMM1NL6XNe6+wx3n1FbW9uZMHpMU+tkTwO0iIgUo3OADcAlZva8mf3OzE7X3HrtS7Z2aT31AqAROUVEsiyTZG+rmX3YzGJmVmJmHwY67EVtZrVhjR5mVgUcB7xEkPSdGR52DnB3+PyecJ1w/78Ltb9esiCrLCsJ1wvyMkREpB3uvtLd/+zuZwMzaBlR+p9m9i8z0xx7aSRviLaeVB1gs5I9EZGsam80zqQPAb8OFwf+S8sAK+0ZDtwQ9tsrAW539/vMbCFwm5n9CHgBuC48/jrgJjNbDKwHzu7UleSRdzbjVM2eiEgxc/cE8HS4fNfMhgAnRBtVfkreEN1pgJbK4M+RTds1SIuISDa1N0DLB4F/hnPddarZJoC7zwX2SbN9CXBAmu07gPd39n3yUbImr7J5gBbV7ImIFBsz+w07D2TmBIOK/dvd/wvcHElgea65GWerPnugmj0RkWxrr2ZvDPC3cI68R4EHgecKtWllT0r20auuCD7e+qZ4lOGIiEhuzEqzbRDwCzP7q7v/qofjKQhNaSZV79/cZ081eyIi2dTeaJw/BX5qZjXAscD/Ab83s5cIRh972N1X9UyYhSVZkFWVxYiVGNsbleyJiBQbd78h3XYz+z3wFPCrHg2oQLQ1qTqoZk9EJNs67LPn7psJpk34B4CZTQVOIuiIrv4IaaT2R6gqi7GjUX32RER6C3ffrrni2tYy9cLO8+wB1KnPnohIVnWY7JnZvmk230UwYIukkRxprDRmVJaVqGZPRKSXCKcO+ijwVtSx5Kt42BukJCUhriwrobTEVLMnIpJlmYzG+VtgX2AuwaTq04AFQH8z+4y7/zOH8RWk1M7nlWUxdjQo2RMRKTZmtpmdB2gxYBswE/hUJEEVgOTUs7GUmj0zo19VmebZExHJskySvbeB89x9ATQ34/wB8HXgTkDJXiupnc+rymLs0AAtIiJFx91roo6hECVr9mKtZvqtqSxlswZoERHJqkySvUnJRA/A3Rea2e7uvkR9EtJL7XxeVR5ju2r2RESKmpntBYwjpVx19zu7eK4vAZ8gqDWcB5xLMHftbcBgYDbwUXdv6F7U0Ugk3tmME4J+e3XbVbMnIpJNmSR7C8zsdwSFDMAHgIVmVgHoVzmN1M7nlaUx9dkTESliZnY9sBdBF4fkiFxO0Pqls+caCVwATA0HerkdOBs4GbjC3W8LR/s8D/hdNuLvafFEsmZv52RPNXsiItmXSbL3ceCzwBfD9f8CXyVI9I7KSVQFrnmAlhKjsjymO5UiIsXtIHefmsXzlQJVZtYI9AFWAEcDHwr33wBcQqEme2kGaAEY0KeMRSs3RxGSiEjRymTqhe3A5eHS2pasR1QEmgdoiZVQVVbC6jrV7ImIFLGnzWyquy/s7oncfbmZ/QJ4E9hO0C9+NrDR3ZPVXm8BI7v7XlFJtFGzN6hvOeu3FmTLVBGRvFXS0QFmdoqZvWBm682szsw2m1ldTwRXqFoP0KJmnCIiRe1GgoRvkZnNNbN5Zja3Kycys4HA6cB4YATQFzixE68/38xmmdmsNWvWdCWEnGsZoKV1slfBxu2Nzc08RUSk+zJpxvkr4L3APHfXL3AGmgdoCade0AAtIiJF7TqCufXm0dJnr6uOBV539zUAZnYncCgwwMxKw9q9UcDydC9292uBawFmzJiRl2V2vI0BWgb3LccdNmxrYEh1RRShiYgUnQ5r9oBlwHwleplrTKnZ61tRytZ6dTgXESlia9z9Hnd/3d2XJpcunutN4CAz62PBkNfHAAuBx4Azw2POAe7uftjRSLRZs1cOoKacIiJZlEnN3teBB8xsJlCf3Ojuv8xZVAUudeqF/lVlbG2I0xhPUNZ6UiERESkGL5jZLcC97FxOdno0Tnd/1szuAJ4HmoAXCGrq7gduM7Mfhduuy0bgUYgnJ1VPU7MHsG5LAwzr6ahERIpTJsnepQQDsVQC5bkNpzg0pky90L+qDIC67Y0MVrMUEZFiVEWQ5B2fsq1LUy8AuPv3gO+12rwEOKBL0eWZ5nn2Wt3/HFStmj0RkWzLJNkb4e7Tch5JEUn2RyiNGf2qgo+4bkeTkj0RkSLk7udGHUMhaXuAlmSyV/+O14iISNdk0q7wATM7vuPDJKkp3jJAS7Jmb5Pm2hMRKSpmdn42jultmidVb9WMc2CfsBmnavZERLImk5q9zwBfNbN6gonUDXB375fTyApY6gAt/SqV7ImIFKmLzGxtO/sNuJBwdEwJJAdoKWlVs1cW9nNXM04RkezJZFL1mp4IpJi0HqAFlOyJiBShmcCpHRzzSE8EUkjaqtmDYJAW1eyJiGRPm8meme3i7ivbe3Emx/RGqQO0DGoeXUx9EEREion66nVN8zx7Je9M9gb1LWf9FiV7IiLZ0l6fvQcyeH0mx/Q6Tc3NOEsY1Lec8lgJKzftiDgqERGR6LU1zx6EyZ5q9kREsqa9Zpx7m1ldO/sNaG9/r9WUSGDWUpDt0r+SFUr2RERE2pxnD2BwdTnPv7mxZwMSESlibSZ77h7ryUCKSWPcKUuZQGiXfpWq2RMRESF1gJZ37htSXcH6rfXEE5625k9ERDonk9E4pZPiicROhdSogVX897X2BmwTEZFCZWYVwPuAcaSUq+7+g6hiymftDdAytF8lCQ/6uQ/tV9nToYmIFJ1M5tmTTmqMO6WxlkJs8i41rKqrZ4P6IYiIFKO7gdOBJmBryiJpNCd7aWruhtVUALCqToOaiYhkg2r2cqApkaAs1pJHTxkeTEn40so6DtltSFRhiYhIboxy9xOjDqJQJNwxA0tTszcsrM1bVbeDPenf06GJiBSdjGr2zCxmZiPMbExyyXVghax1X4OpI4Jkb86yjRFFJCIiOfSUme0ZdRCFIp7wtE04ISXZ26x+7iIi2dBhsmdmXwBWEUwMe3+43JfjuApaU9wpTUn2hlRXMHV4Px5/eU2EUYmISI4cBsw2s0VmNtfM5pnZ3KiDyldx97Rz7AEMqS6nxGCVBjUTEcmKTJpxXghMdvd1uQ6mWMT9naOIHTNlKL99/DXWbalncHVFRJGJiEgOnBR1AIUk0U7NXmmshCHVFeqzJyKSJZk041wGbMp1IMUknti5Zg/glL1GEE84dz6/PKKoREQkF9x9KTAAODVcBoTbJI14Iv3gLEnD+lWqGaeISJZkkuwtAR43s4vN7MvJJdeBFbKmNPMDTd6lhv3GDuTW597EwzmGRESk8JnZhcDNwNBw+UvYBULSSLjT3hR6w/qpZk9EJFsySfbeJOivVw7UpCzShnjcKU0zW+yHDhjDkrVbeWbJ+giiEhGRHDkPONDdv+vu3wUOAj4ZcUx5q6MJ04f2q2R1nWr2RESyocM+e+7+fQAzqw7Xt+Q6qEKXrmYP4OQ9h3PJvQu4fdYyDt5tcASRiYhIDhgQT1mPh9skjXT92lMNq6lk3dYG6pviVJTGejAyEZHik8lonNPM7AVgAbDAzGab2R65D61wxROJnSZVT6oqj3HG9JE8MG8Fm7Y1RhCZiIjkwJ+AZ83sEjO7BHgGuC7akPJXPO6UtDFAC8Au/YNBzFarKaeISLdl0ozzWuDL7j7W3ccCXwH+0NGLzGy0mT1mZgvNbEHYp4GwMFxuZnPC5eSU11xsZovD4atP6OpFRa0p0XZB9oH9R1PflODuFzVQi4hIMXD3XwLnAuvD5Vx3/1WkQeWxjmr2Rg7oA8BbG7b3VEgiIkUrk6kX+rr7Y8kVd3/czPpm8Lom4Cvu/ryZ1RDMQfRIuO8Kd/9F6sFmNhU4G9gDGAH8y8wmuXucApNuNM6kaSP7M21kP259bhkfO3hczwYmIiJZY2b93L3OzAYBb4RLct8gd1cH7TQS7dwQBRg1sAqAtzZsA9TlQUSkOzIajdPMvmNm48Ll2wQjdLbL3Ve4+/Ph883AS8DIdl5yOnCbu9e7++vAYuCADOLLOx11Pn//fqN5aUUdr6za3INRiYhIlt0SPs4GZqUsyXVJI+6etqtD0ogBVZipZk9EJBsySfb+D6gF7gyX2nBbxsxsHLAP8Gy46fNmNtfMrjezgeG2kQRz+iW9RfvJYd6KJ9ovyE7acxfM4P65K3owKhERySZ3PyV8HO/uu6Ys491916jjy1fxdiZVBygvLWGXfpVK9kREsqDDZM/dN7j7Be6+b7hc6O4bMn2DcBTPvwNfdPc64HfAbsB0YAVweWcCNrPzzWyWmc1as2ZNZ17aY4LRONv+aIfWVHLg+EHcP2+F5twTESlwZvZoJtskkHCnpL2J9giacgbNOEVEpDvazEjM7Ffh471mdk/rJZOTm1kZQaJ3s7vfCeDuq9w97u4JgoFekk01lwOjU14+Kty2E3e/1t1nuPuM2traTMLoce312Ut6914jWLx6C6+s0kwWIiKFyMwqw/56Q8xsoJkNCpdxFGjLlJ7QUc0ewKiBfVSzJyKSBe0N0HJT+PiLdo5pk5kZwdDTL4UjlSW3D3f3ZPvF9wDzw+f3ALeY2S8JBmiZCDzXlfeOWlvz7KU6cY9d+N7d87l/7ttM3mVyD0UmIiJZ9CngiwRl1mxa5tarA66KKKa8F0/QYc3e6IFV3D1nO43xBGWxTHqciIhIOm0me+4+O3w63d1/nbovnEZhZgfnPhT4KDDPzOaE274JfNDMpgNOMHLZp8L3W2BmtwMLCUby/FwhjsQJ4Tx7HRRktTUVzBg3iH8uXMWXj1eyJyJSaMKy8ddm9gV3/03U8RSKhDsd5W+jBvYh4bBy0w5GD+rTM4GJiBShTKZeOAf4dattH0+zbSfu/iQtdzlTPdDOay4FLs0gprzWlOi4PwLAsVOG8uMHXuatDdsYNVCFmYhIIXL335jZNGAqUJmy/cboospfmTXjDKZfeHP9NiV7IiLd0F6fvQ+a2b3A+Fb99R4jmDRW2pBJnz2AY6cMA+DfL6/OdUgiIpIjZvY94DfhchTwM+C0SIPKY5kM0DK+NpjOd8narT0RkohI0WqvZu8pgtEyh7DziJmbgbm5DKrQdTTPXtKutdXsOqQvjyxcpQnWRUQK15nA3sAL7n6umQ0D/hJxTHkrk5q9XfpV0qc8xpI1GsRMRKQ72uuztxRYChzcc+EUh0xr9gCOmTKUG55aypb6JqorMmlVKyIieWa7uyfMrMnM+gGr2Xl0aUkRz6Crg5mxa21fXlujmj0Rke7ocIgrMzvIzP5nZlvMrMHM4mZW1xPBFaqO5tlLdeyUYTTEE/znlfycM1BERDo0y8wGEEwnNBt4Hng60ojyWMI7rtkD2K22WjV7IiLdlElGchXwQeBVoAr4BHB1LoMqdJ2p2dtv7ED6V5Xxr5fUb09EpNCE0wz9xN03uvvvgeOAc9z93IhDy1sZd3UYUs3yjdvZ0ViQA3OLiOSFjKqf3H0xEAsnQ/8TcGJuwypsTfFERgUZQGmshCMn1zLzldUkEp7jyEREJJvc3UkZZdrd33B39WtvR9w7nmcPYLehfXGH1zVIi4hIl2WS7G0zs3Jgjpn9zMy+lOHreq3O1OwBHDm5lrVbGli4Qq1jRUQK0PNmtn/UQRSKRMKJZVBE7jqkGoDX1JRTRKTLMknaPhoe93lgK0Gn8/flMqhC15RwYpmUZKF3TawFYKb67YmIFKIDgafN7DUzm2tm88xMtXttyLQZ5/ghfTGD11arZk9EpKvaHf7RzGLAj939w8AO4Ps9ElWBy2RY6VRDqiuYNrIfMxet4XNHTchhZCIikgMnRB1AIUm4U5JBGVlVHmPsoD4sWqVWLyIiXdVuzZ67x4GxYTNOyVDcO9eME+CISbXMfnMDdTsacxSViIjkyI/cfWnqAvwo6qDyVaY1ewBTR/RjwdtK9kREuiqTZpxLgP+a2XfM7MvJJdeBFapEwnEn46kXko6YNJR4wnlq8docRSYiIjmyR+pK2Cpmv4hiyXtx73ievaSpw/uxdN02NutGqIhIl2SSkbwG3BceWxMu1bkMqpA1hSNqlnaizx7APmMGUFNRqn57IiIFwswuNrPNwF5mVhcumwkmVb874vDyVqITXR2mjugHwMsrN+cyJBGRotVun73QQnf/W+oGM3t/juIpePEw2cu0iUpSWayEQycMYeaiNbg71ok+fyIi0vPc/SfAT8zsJ+5+cdTxFIq4Z96Mc48R/QFY+HYd+48blMuwRESKUiY1e+kKMBVqbWhKJAA63WcP4IjJtby9aQeLV2uYaRGRQuHuF5vZSDM7xMwOTy5Rx5WvEgkyGqAFYGhNBYP7lrPg7U05jkpEpDi1WbNnZicBJwMjzezKlF39gKZcB1aoulqzB3D4pJYpGCYOq8lqXCIikhtmdhlwNrAQiIebHXgisqDyWDBAS2bHmhlTR/Rj/nIN0iIi0hXtNeN8G5gFnAbMTtm+GfhSLoMqZM199rqQ7I0cUMWEodXMfGUNn3jXrtkOTUREcuM9wGR3r486kELQmWacAPuMHsBVjy1ma30TfSsy6X0iIiJJbf5quvuLwItmdou7axisDCVr9jIdaay1IybVctMzS9neEKeqPJbN0EREJDeWAGWAkr0MJBKZzbOXtO/YgSQcXly2kUMmDMlhZCIixSeThhQHmNkjZvaKmS0xs9fNbEnOIytQ8W7U7EGQ7DU0JXjm9XXZDEtERHJnGzDHzK4xsyuTS9RB5atO1+yNGQjA7KUbchWSiEjRyqQ9xHUEzTZn09IXQdrQ0mevc/PsJR0wfhCVZSXMXLSGoyYPzWZoIiKSG/eEi2Qg3smavf5VZUwaVs3sN5XsiYh0VibJ3iZ3fzDnkRSJ7vTZA6gsi3HQroN5QvPtiYgUBHe/wcyqgDHuvqi75zOzAcAfgWkEA738H7AI+CswDngDOMvdCzL7SSQ6V7MHsN/Ygdw/d0XQBLSL5auISG+USfXTY2b2czM72Mz2TS45j6xAxcOpF7oyGmfSEZNqWbJ2K2+u25atsEREJEfM7FRgDvBQuD7dzLpT0/dr4CF33x3YG3gJuAh41N0nAo+G6wWps804AWaMHUTdjiYWrtConCIinZFJzd6B4eOMlG0OHJ39cApfd2v2IGUKhlfX8NHBY7MSl4iI5MwlwAHA4wDuPsfMujSkspn1Bw4HPh6eqwFoMLPTgSPDw24I3+sbXQ85Op2ZZy/pXRODgVmeXLyWaSP75yIsEZGi1GHNnrsflWZRoteGpnjX59lL2nVIX0YNrFJTThGRwtDo7q1n/U508VzjgTXAn8zsBTP7o5n1BYa5+4rwmJXAsC6eP3JBzV7nXjO0XyW771LDf15VuSgi0hkd/tya2TAzu87MHgzXp5rZebkPrTA1j8YZ63qyZ2YcMamWpxavpaGpq38viIhID1lgZh8CYmY20cx+AzzVxXOVAvsCv3P3fYCttGqy6e5O0MLmHczsfDObZWaz1qzJz8QonnBinazZg6B273+vb2B7g8aKExHJVCb31v4MPAyMCNdfAb6Yo3gKXrIZZ2ebqLR2xKRatjbENdS0iEj++wKwB8E8e7cAm+h6OfkW8Ja7Pxuu30GQ/K0ys+EA4ePqdC9292vdfYa7z6itre1iCLmT6MZctIdNrKUhnuBZTU0kIpKxTJK9Ie5+O2GTFHdvQlMwtCnhyT57XZt6IemQCUMoLTFmqimniEhec/dt7v4td98/XL7t7ju6eK6VwDIzmxxuOgZYSDC1wznhtnOAu7sdeATiYRnZlZq9A8OpiR59KW2eKyIiaWSSkWw1s8GETUbM7CCCu5aSRjb67AFUV5Sy39iBSvZERPKcmT0STpeQXB9oZg9345RfAG42s7nAdODHwGXAcWb2KnBsuF5w4t2o2assi3HU5KE8tGBlcw2hiIi0L5PROL9McEdxNzP7L1ALnJnTqApYNvrsJR0xuZafPbSI1XU7GNqvstvnExGRnBji7huTK+6+wcyGdvVk7j6HnUfATjqmq+fMF8nWL129IXrynsN5cP5KZr+5gf3HDcpmaCIiRSmT0TifB44ADgE+Bezh7nNzHVihasrCPHtJR4RTMDzx6tpun0tERHImYWZjkitmNpY2BlDp7ZI3RLvSjBPgqN2HUlFawgPzVnR8sIiIZDQa5+eAandf4O7zgWoz+2zuQytM8SzMs5c0dXg/amsq1JRTRCS/fQt40sxuMrO/AE8AF0ccU14K74d2qRknBF0cjphUy/1zV9AU12jVIiIdyaTP3idbN08BPpmziApccjTObNTsmRmHT6zlP6+uaU4iRUQkv7j7QwQjZv4VuA3Yz92702evaLUM0NL1c7x/xmhWb67n8UW6ESoi0pFMkr2YWUt7CzOLAeW5C6mwtdTsdW80zqTDJw1h47ZG5r61MSvnExGRnKgA1gN1wFQzOzziePJSPAs3RI+aXMvQmgpu+9+ybIUlIlK0Mhmg5SHgr2Z2Tbj+qXCbpNFSs5ed871rYi1m8MQra9lnzMDsnFRERLLGzH4KfABYQDhNEUGfvSciCypPJQdo6WozToDSWAln7jeKa55Ywqq6HQzTAGYiIm3KJCX5BvAY8JlweRT4ei6DKmTx5gFaspPtDepbzl6jBjDzFc0rJCKSp84AJrv7u9391HA5Leqg8lF3B2hJ+sD+o0m485dnlmYjLBGRopXJaJwJd/+du58ZLte4e4eTqpvZaDN7zMwWmtkCM7sw3D4onJPo1fBxYLjdzOxKM1tsZnPNbN/uX17PS/YXz8YALUlHTKplzrKNbNzWkLVziohI1iwByqIOohB0Z569VGMH9+WEqbtw49NL2VrflI3QRESKUiajcR4aJmWvmNkSM3vdzJZkcO4m4CvuPhU4CPicmU0FLgIedfeJBLWEF4XHnwRMDJfzgd914Xoil6zZ625BluqISbUkHJ5crCkYRETy0DZgjpldE960vNLMrow6qHzUPM9eN2v2AD51xK5s2t6ovnsiIu3IpM/edcCXgNlAhzV6Se6+AlgRPt9sZi8BI4HTgSPDw24AHidoKno6cKO7O/CMmQ0ws+HheQpGUxanXkjae1R/+lWWMnPRGk7Za0TWzisiIllxT7hIB7IxQEvSPmMGcuD4QfzxP0v48IFjqCyLdfucIiLFJpOOZZvc/UF3X+3u65JLZ97EzMYB+wDPAsNSEriVwLDw+Ugg9fbcW+G2gpKcIaEkC3ctk0pjJbxrYi1PvLoGd03BICKST9z9BuBWgpuis4Fbwm3SSjYGaEl14TETWbFpBzc9rb57IiLpZJLsPWZmPzezg81s3+SS6RuYWTXwd+CL7l6Xui+sxetU9mJm55vZLDObtWZN/s2xk0j2R8hergcETTlX1dWzaNXm7J5YRES6xcyOBF4FrgZ+C7yiqRfSS/Zrz0YzToBDJgzhiEm1XPXYYjZta8zKOUVEikkmyd6BwAzgx8Dl4fKLTE5uZmUEid7N7n5nuHmVmQ0P9w8HksNMLgdGp7x8VLhtJ+5+rbvPcPcZtbW1mYTRo7LZRCXV4ZOCa52pSWRFRPLN5cDx7n6Eux8OnABcEXFMeSme5emJAL5x4u7U7Wjk6scXZ++kIiJFIpPROI9Ksxzd0evCidivA15y91+m7LoHOCd8fg5wd8r2j4Wjch5E0Hy0oPrrQfabqCTt0r+S3XepYeYrSvZERPJMmbsvSq64+ytodM60msvILHZ1mDqiH+/fbxTXP/k6i1aq9YuISKpMRuPsb2a/TDadNLPLzax/Buc+FPgocLSZzQmXk4HLgOPM7FXg2HAd4AGC4asXA38APtuVC4paNkcaa+3wSbXMemODhpkWEckvs8zsj2Z2ZLj8AZgVdVD5KFetXy4+aQr9qsr45j/mNXenEBGRzJpxXg9sBs4KlzrgTx29yN2fdHdz973cfXq4PBAO8HKMu09092PdfX14vLv759x9N3ff090LsqBM9kfI5l3LpCMm1dIQT/DMkk6NjyMiIrn1GWAhcEG4LAy3SSvxHLV+Gdi3nG+ePIXZSzdwy3NvZvXcIiKFLJNkbzd3/567LwmX7wO75jqwQtXSjDP7554xbiBVZTE15RQRyS+lwK/d/b3u/l7gSkDzAKSRrHXLReuX9+07kkMnDObHD7zE0nVbs35+EZFClElKst3MDkuumNmhwPbchVTY4jksyCpKYxyy22AleyIi+eVRoCplvQr4V0Sx5LWmHDXjBDAzfn7m3pSWGF/66xyakk1tRER6sUySvc8AV5vZG2b2BnAV8OmcRlXAmvvs5aAgAzhici1L123j9bW6aykikicq3X1LciV83ifCePJWy/REuSkjRwyo4kfv2ZPn39zIbx9/LSfvISJSSDIZjXOOu+8N7AXs5e77uPuLuQ+tMCULMstRQXbkpKEA/Pvl1R0cKSIiPWRr6vyzZrYfagGTVrLPXmksN2UkwGl7j+D06SP49aOvMnvphpy9j4hIIchkNM4fm9kAd69z9zozG2hmP+qJ4ApR3D1ntXoAYwb3YeLQav798qqcvYeIiHTKF4G/mdl/zOxJ4K/A56MNKT/Fc1yzl/TDM6YxckAVn7/ledZvbcjpe4mI5LNMmnGe5O4bkyvuvgE4OWcRFbh4Ijf99VIdPWUozy5ZT92Oxpy+j4iIdMzd/wfsTtDt4dPAFHefHW1U+SnXXR2S+lWW8dsP78u6LQ18+fY5mo5BRHqtTJK9mJlVJFfMrAqoaOf4Xs3dczISZ6pjpwyjKeH855W1uX0jERHJiLs3uvv8cNGduDYkx0zJ9U1RgGkj+/OdU6fy+KI1/G6m+u+JSO+USVpyM/ComZ1nZucBjwA35DaswhVPeM6bp+wzegAD+pTx6EtqyikiIoWjuRlnjm+KJn3kwDGcuvcILv/nIs1RKyK9UiYDtPwU+BEwJVx+6O4/y3VghSrunvM7lqWxEo6aPJTHFq1uLjhFRETyXU8140wyM37y3j0ZN7gvF9z6Ams21/fI+4qI5IuM7q25+0Pu/tVweTjXQRWyRMIp6YFC7Ojdh7JhWyNzlmmkMRGRqJnZSDM7xMwOTy5Rx5SPcjkXbVuqK0q5+sP7sml7I1/86wu6SSoivUoPNaToPRLeM3csD59US2mJ8a+XNAWDiEiUzOynwH+BbwNfC5evRhpUnkrW7PXETdFUU4b34wen78F/F6/jd48v7tH3FhGJUmnUARSbuDs9UYb1rypj/3GD+PdLq/nGibvn/g1FRKQtZwCT3V1tBDsQRc1e0lkzRvPUa+u44l+vcuCug9l/3KAej0FEpKdlVLNnZlVmNjnXwRSDRA8M0JJ0zJShLFq1mWXrt/XI+4mISFpLgLKogygEzcleD9fsQdB/70dnTGPUwCouvPUFNm7T/HsiUvwymVT9VGAO8FC4Pt3M7slxXAUrnsjtpOqpjpkyDIB/v6ymnCIiEdoGzDGza8zsyuQSdVD5KKpmnEk1lWX85oP7sGZLPV+/Yy7u6r8nIsUtk5q9S4ADgI0A7j4HGJ+ziApcwumxmr3xQ/qy65C+/EtTMIiIROke4IfAU8DslEVa6cl59tqy16gBfOPE3fnnwlXc9MzSyOIQEekJmfTZa3T3TbbzD7NuhbUh0QOTqqc6ZspQbnhqKVvqm6iuUBdMEZGe5u6aezZDce/Zefbact5h43nqtXX86L6X2G/sQPYY0T/agEREciSTn9sFZvYhIGZmE83sNwR3LyWNeCL38+ylOmbKMBriCWYuWtNj7ykiImBmt4eP88xsbusl6vjyUSLCAVpSmRm/eP/eDOxbxoW3zWFHYzzSeEREciWTZO8LwB5APXArUAd8MYcxFbS498w8e0kzxg5kUN9yHl6wssfeU0REALgwfDwFODXNIq1EOUBLa4P6lvPzM/dm8eot/OLhRVGHIyKSEx22+3P3bcC3gG+ZWQzo6+47ch5ZgXLv2Zq90lgJx00ZxgPzVlDfFKeiNNZj7y0i0pu5+4rwUR2/MhT1AC2tHT6plo8eNJbr/vs6x0wZxsG7DY46JBGRrMpkNM5bzKyfmfUF5gELzexruQ+tMMV7cOqFpBOn7cLm+iaeem1dj76viEhvZmabzawuzbLZzOqiji8fRTnPXlsuPnl3xg7qw1f/9iKbdzRGHY6ISFZl0oxzqrvXEUwa+yDBSJwfzWVQhSye6Pk7lodMGEx1RSkPz1dTThGRnuLuNe7eL81S4+79oo4vHyUHaMmHZpxJfcpLufys6azYtJ0f3rcw6nBERLIqk2SvzMzKCJK9e9y9EY3G2aaEO7EeHmWsojTG0bsP5Z8LVzXfNRURkZ5hZmPSLVHHlY+SA7T0dAuYjuw3diCfOXI3bp/1Fo8s1HRGIlI8MklLrgHeAPoCT5jZWIJBWiSNhPd8M04ImnKu39rA/95Y3+PvLSLSy92fsjwKLCFoCSOtNM+zl0c1e0kXHjOJKcP7cfGd89i4rSHqcEREsqLDZM/dr3T3ke5+sgeWAkf1QGwFKYo+ewBHTKqlorSEh9SUU0SkR7n7ninLROAA4Omo48pHzfPs5V+uR3lpCZe/f282bmvgB2rOKSJFIpMBWvqb2S/NbFa4XE5QyydpBM04e74U61tRyuGTanl4wUrc1ZRTRCQq7v48cGDUceSjRMIpsWCeu3w0dUQ/PnPkbtz5/HIeX7Q66nBERLotk2ac1wObgbPCpQ74Uy6DKmQ9Pal6qhP32IUVm3Yw961Nkby/iEhvZGZfTlm+ama3Am9HHVc+ikd0Q7QzPn/0BCYMreZb/5jPlvqmqMMREemWTJK93dz9e+6+JFy+D+ya68AKVcIhqhuWx0wZSmmJ8aCacoqI9KSalKUCuA84PdKI8lQioq4OnVFRGuOn79uLtzdt56cPvhx1OCIi3ZJJsrfdzA5LrpjZocD23IVU2BKJ6O5aDuhTzqEThnDf3LfVlFNEpIe4+/eTC/AT4F533xF1XPkoHmEZ2Rn7jR3IuYeM56ZnlvLsEs1hKyKFK5Nk7zPA1Wb2hpm9AVwFfDqnURWwqJuonLr3CN7asJ0Xlm2MLAYRkd7EzG4xs35m1heYDyw0s69FHVc+int0XR0666snTGLMoD5cdOc8djTGow5HRKRLMhmNc4677w3sBezl7vu4+4u5D60wJTza+YNO2GMY5aUl3DNH3UVERHrIVHevI5iP9kFgPPDRSCPKU4mEU1IANXsQTLZ+2Xv35PW1W7niX69EHY6ISJdkMhrnj81sgLvXuXudmQ00sx/1RHCFKDnSWFRqKss4evJQ7p+3QhOsi4j0jDIzKyNI9u5x90ZAP8BpRN36pbMOmTCEs/cfzR+eWML85Rr8TEQKTybNOE9y943JFXffAJycs4gKXD70Rzh17xGs2VyvfgYiIj3jGuANgmmJnjCzsQQjV0sr8US0rV+64uKTpzC4uoKL75xHU3JWeBGRApFJshczs4rkiplVEYw2JmkkPPqRxo6ZMpS+5THueVFNOUVEcs3dr3T3ke5+sgeWAkd155xmFjOzF8zsvnB9vJk9a2aLzeyvZlaeleB7WDCIWdRRdE7/qjK+d+pU5i3fxI1PL406HBGRTsnkJ/dm4FEzO8/MzgMeAW7IbViFK6pJ1VNVlsU4fo9deHD+ShqadBdSRCSXzKzCzD5kZt80s++a2XeBb3bztBcCL6Ws/xS4wt0nABuA87p5/kgU0gAtqd6953COnFzL5f9cxNsbNSC5iBSOTAZo+SlwKTAlXH7o7j/LdWCFKp4ncwiduvdwNm1v5IlX1kQdiohIsbubYF69JmBrytIlZjYKeDfwx3DdgKOBO8JDbiDoH1hwCmmAllRmxg9Pn0bC4bt3L9D0RiJSMEozOcjdHyQYYSxjZnY9cAqw2t2nhdsuAT4JJDOQb7r7A+G+iwnuVMaBC9z94c68X75IOHlRkL1rYi2D+pbzjxeWc+zUYVGHIyJSzEa5+4lZPN+vgK8TTNIOMBjY6O5N4fpbwMh0LzSz84HzAcaMGZPFkLKj0AZoSTV6UB++dNxEfvzAyzy8YCUnThsedUgiIh3KZDTOzWZWFy47zCxuZpl0PP8zkK7wu8Ldp4dLMtGbCpwN7BG+5rdmFsv8MvJHPOHE8qAcK4uVcPr0ETyycBUbtzVEHY6ISDF7ysz2zMaJzCx5k3R2V17v7te6+wx3n1FbW5uNkLIqKCPzoJDsov87dDxThvfje/csoG5HY9ThiIh0KJNmnDXu3s/d+wFVwPuA32bwuieA9RnGcTpwm7vXu/vrwGLggAxfm1cSnj9NVN6/32ga4gnu1px7IiK5dBgw28wWmdlcM5tnZnO7eK5DgdPM7A3gNoLmm78GBphZsjXOKGB5d4OOQj6VkV1RGivhJ+/dk9Wb67n84UVRhyMi0qFOjYkVjjJ2F3BCN97z82FheL2ZDQy3jQSWpRzTbhMVM5tlZrPWrMm//miJPOmzBzB1RD+mDu/H32Yv6/hgERHpqpOAicDxwKkEXRhO7cqJ3P1idx/l7uMIWrz8290/DDwGnBkedg5BP8GCU+g1ewDTRw/gnIPHceMzS3nhzQ1RhyMi0q5MmnG+N2U508wuA3Z08f1+B+wGTAdWAJd39gR530Qlz0Yae/+MUcxfXsdLKzTlk4hINplZv/Dp5jaWbPoG8GUzW0zQh++6LJ+/R8QT+dGvvbu+cvwkhtVUcvGd82jU3Hsikscyqdk7NWU5gaAAO70rb+buq9w97u4J4A+0NNVcDoxOObRgm6jkW0F2+vSRlMWMv816K+pQRESKzS3h42xgVvg4O2W9W9z9cXc/JXy+xN0PcPcJ7v5+d6/v7vmjEExPFHUU3VdTWcb3T9+Dl1du5ronX486HBGRNnU4Gqe7n5utNzOz4e6+Ilx9DzA/fH4PcIuZ/RIYQdAc5rlsvW9P8jwryAb1LefYKcO4a85yLjppd8pL8yg4EZEClpKIjY86lkJRDM04k07YYxeOnzqMX/3rFd6953BGD+oTdUgiIu/QZrJnZr8B2pxIxt0vaO/EZnYrcCQwxMzeAr4HHGlm08PzvgF8KjzXAjO7HVhIME/R59w93pkLyRdxz58+e0lnzRjNg/NX8q+XVnHynhoqWkREolHoA7S0dslpe3DcL2fy7bvm8+dz98fyrPwXEWmvZq9bTVDc/YNpNrfZx8DdLyWYvL2g5cuk6qkOn1TLyAFV3Pj0G0r2REQkMsVUswcwYkAVXz1hMt+/dyH3zl3BaXuPiDokEZGdtJnsufsNPRlIsUgk8m/C2FiJ8eGDxvCzhxbx6qrNTBxW0/GLREREsiyeKK6aPYCPHTyOu15Yzg/uXcARE2vp36cs6pBERJplMhrnvWZ2T6vlJjO70MwqeyLIQpJw8i7ZA/jAjNGUx0q46ZmlUYciIlKUzGyomY1JLlHHk48SeTZidTbESowfv3dPNmxr5LKHXoo6HBGRnWQyWscSYAvB6Jl/AOoIRuScFK5Lirg7+ViODa6u4JS9hnPn88vZUt8UdTgiIkXDzE4zs1eB14GZBH3SH4w0qDwVz8PWL9mwx4j+nHfYeG59bhnPvb4+6nBERJplkuwd4u4fcvd7w+UjwP7u/jlg3xzHV3ASedwf4SMHj2VLfRP/eKEgZ7UQEclXPwQOAl4JR+Y8Bngm2pDyU9zza3qibPrisRMZOaCKb/5jHvVNBTnGnIgUoUySverU5ijh8+pwtSEnURWwuOfvXct9Rg9g2sh+3PjUG7i3OdCqiIh0TqO7rwNKzKzE3R8DZkQdVD4KbohGHUVu9Ckv5UdnTGPx6i1cM3NJ1OGIiACZJXtfAZ40s8fM7HHgP8BXzawvoEFcUrg77uTdaJxJZsb/HTqeV1dv4fFFa6IOR0SkWGw0s2rgCeBmM/s1sDXimPJSsTbjTDpq96G8e6/hXPXYYpas2RJ1OCIiHSd77v4AwSTnXwQuBCa7+/3uvtXdf5Xb8ApLIqwsy9dkD+DUvUcwvH8l1zzxWtShiIgUi9OB7cCXgIeA14BTI40oTyXycC7abPveqVOpKC3hW/+Yr1Y0IhK5TGr2cPd6d38xXHbkOqhCFQ+zvVhGn2o0ymIlnHfYeJ5Zsp4Xl22MOhwRkYIX3vyMA32Ae4G/APorP41ir9kDGFpTyUUn7c7TS9bx9+fVR15EopXHaUnhSYR38PK98/nZB4yhprKUa59QnwIRke4ys0+Z2UpgLjALmB0+SitxL7559tL54P5jmDF2IJfev5D1WzW8gYhER8leFiWTvXwdjTOpuqKUjxw0lgfnr2DpOnUrERHppq8C09x9nLvv6u7j3X3XqIPKR/k8YnU2lYRz722pb+JH9y+MOhwR6cXaTPbMbN/2lp4MslAkm3EWQn+Ecw8ZR2lJCdeodk9EpLteA7ZFHUQhyOcRq7Nt0rAaPnX4btz5/HL+u3ht1OGISC9V2s6+y8PHSoIhpF8EDNiLoHnKwbkNrfAkEsFjITRRGdqvkg/sP5rb/vcmnzliN0YP6hN1SCIihepi4CkzexaoT2509wuiCyk/xePFP0BLqs8fPYH75r7Nt/4xj4e+eDiVZbGoQxKRXqbNmj13P8rdjwJWAPu6+wx33w/YB1CP4zTizc04Iw4kQ589ajfMjKsfWxx1KCIihewa4N8EE6nPTlmklaBmL+ooek5lWYxL37Mnb6zbxlX/VlkrIj2vvZq9pMnuPi+54u7zzWxKDmMqWM199gqgZg9geP8qPnTAGG56ZimfOXI3xg7uG3VIIiKFqMzdvxx1EIUgnoBYSS/K9oBDJwzhvfuM5PczX+O06SOYNKwm6pBEpBfJ5Bd3rpn90cyODJc/EIw4Jq0kwj57VkBNVD575G6Ulhi/0R1HEZGuetDMzjez4WY2KLlEHVQ+SvSymr2kb717CjWVpVx857zmvxVERHpCJj+55wILCCZUvxBYGG6TVuIFVrMHQd+9jxw0ljuff4vX1myJOhwRkUL0QcJ+e7Q04dTUC2nEe8lonK0Nrq7gmydPYfbSDdz6vzejDkdEepEOk71wEvWrge8C3wGu0sTq6TVPql5gBdlnjtyNPuWl/OSBl6MORUSk4IRTLbReNPVCGolE75hnL50z9xvFwbsO5rIHX2Z1nf6MEpGe0WGyZ2ZHAq8CVwG/BV4xs8NzG1ZhCiv2Cq4gG1JdwWeO3I1/vbSKp19bF3U4IiIFxczeb2Y14fNvm9mdZrZP1HHlo7j3zpo9CLp4XPqeadQ3Jfj+fZp7T0R6RibNOC8Hjnf3I9z9cOAE4IrchlWYWubZiziQLjjvsPGMHFDFpQ8sVH8CEZHO+Y67bzazw4BjgeuA30ccU16KJ3rPPHvp7FpbzeePmsD9c1fw2Murow5HRHqBTJK9MndflFxx91eAstyFVLgKsc9eUmVZjK+dMJn5y+v4xwuaWUNEpBPi4eO7gWvd/X6gPMJ48lbCe28zzqRPH7EbE4ZW8+275rOtoSnqcESkyGWS7M1KMxqnOp6nkWiu2SvMguy0vUew16j+/Ozhl9m8ozHqcERECsVyM7sG+ADwgJlVkFn52uvEE05pL0/2yktL+Ml792T5xu384uFXog5HRIpcJoXRZwhG4LwgXBaG26SVZOvHQk32SkqM75+2B6s313PFI69GHY6ISKE4C3gYOMHdNwKDgK9FGlEecncSXrhlZDbtP24QHzloDH966nVmL90QdTgiUsQyGY2znmBwlu8RjMh5VbhNWmkejbOA7+fuM2YgHzxgDH9+6nXmL98UdTgiInnP3be5+53AJjMbQ9DVQcMbt5IsI3t7zV7SRSdNYUT/Kr5+x4vsaIx3/AIRkS7QaJxZlPDCbsaZ9I0Tdmdgn3K+fdd8DdYiItIBMzvNzF4FXgdmho8PRhtV/mlKdnVQsgdAdUUpP3nvnry2Ziu/flStaUQkNzQaZxa11OwVdkHWv08Z33r3FOYs28jNz2nyVxGRDvwQOAh4xd3HE4zI+Uy0IeWf5A1R1ey1OHxSLWfNGMW1Tyxh7lsbow5HRIqQRuPMouaavSIoyN6zz0gOnTCYyx54iWXrt0UdjohIPmt093VAiZmVuPtjwIyog8o3TUVyQzTbvvXuqQypLufrd8yloSkRdTgiUmQ0GmcWFUszTggmf73svXthZnztjhfVnFNEpG0bzawaeAK42cx+DWyNOKa8k1Cyl1b/qjJ+/J49eXnlZq5+bHHU4YhIkdFonFkUD2/IxYog2QMYPagP3zllCs8sWc+NT78RdTgiIvnqdGA78CXgIeA14NRII8pDqtlr2zFThnHG9BFc/dhiXlpRF3U4IlJEMhqN091/6e7vDZcrNBpnevHmzucRB5JFZ80YzVGTa7nsoZd5fa1uVIuItObuW9097u5N7n6Du18ZNuuUFKrZa9/3Tt2DAX3K+NodL9IYV3NOEcmONtMSM5tnZnPbWnoyyEKRbMZZLDV7EDbnfN9eVJTG+MKtz1PfpOGhRUQAzGyzmdWFj8nnyXVVz7TSXLNXRGVkNg3sW84PTp/G/OV1XDPztajDEZEiUdrOvlN6LIoi0ZzsFdldy2H9Kvn5mXtx/k2z+ckDL3PJaXtEHZKISOTcvSbqGApJsYxYnUsn7zmcd+85nF8/+ipH7T6UPUb0jzokESlwbdbsufvS1gtBh/M3w+fSSrIgsyK8a3n8Hrtw7qHj+PNTb/DQ/JVRhyMiEjkzqzSzL5rZVWZ2vpm1dwO111Oyl5kfnTGNgX3K+dJf52iydRHptvaacR5kZo+b2Z1mto+ZzQfmA6vM7MSeC7FwFGvNXtLFJ01hr1H9+fodL2o6BhERuIFgioV5wMkE89JKG+JFXkZmy8C+5fz0zL14ZdUWLv/noo5fICLSjvaGErkK+DFwK/Bv4BPuvgtwOPCTHoit4BTbaJytlZeWcNUH98WBT900m20NTVGHJCISpanu/hF3vwY4E3hX1AHlM9XsZe6oyUP58IFj+OOTr/PMEo31IyJd116yV+ru/3T3vwEr3f0ZAHd/OZMTm9n1ZrY6rBFMbhtkZo+Y2avh48Bwu5nZlWa2OBwAZt/uXFRUinE0ztbGDO7DlWfvw0sr6/ja3+birvn3RKTXakw+cXfd/epAXAO0dMq33j2FsYP68JXbX2TzjsaOXyAikkZ7aUnquL/bW+3L5C/8PwOtm3teBDzq7hOBR8N1gJOAieFyPvC7DM6fd4q9GWfSUbsP5aITd+f+eSs0AayI9GZ7p47ACeyl0Tjbppq9zulTXsovPzCdFZu28/17F0YdjogUqPaSvb3TFGDJ9T07OrG7PwGsb7X5dII+DoSPZ6Rsv9EDzwADzGx4Zy4kHySTvZJecNfy/MN35T37jOQX/3yFhxdowBYR6X3cPebu/cKlxt1LU573izq+fKNkr/P2HTOQzx45gTtmv6WyVkS6pL3ROGNpCrDkelkX32+Yu68In68EhoXPRwLLUo57K9z2DuGIZ7PMbNaaNWu6GEZuNDfj7AXJnpnxk/fuyd6jB3DBrS8we2nrvF5ERKSFBmjpmguOmci0kf246O9zWblpR9ThiEiBiax3mQedvTrd4cvdr3X3Ge4+o7a2NgeRdV1vacaZVFkW4/pzZjC8fyXn3TCLxau3RB2SiIjkKdXsdU15aQm/Pnsf6psSXHjbC82fo4hIJno62VuVbJ4ZPq4Oty8HRqccNyrcVlCKfTTOdAZXV3Dj/x1IaUkJ51z/HKvqdNdRRETeScle1+1WW80PTp/Gs6+vV195EemUnk727gHOCZ+fA9ydsv1j4aicBwGbUpp7FoxELxiNM50xg/vw53P3Z+O2Bj523XOs39oQdUgiIpJnNBpn97xv35G8Z5+R/Opfr/Dc6+o6ISKZyVlaYma3Ak8Dk83sLTM7D7gMOM7MXgWODdcBHgCWAIuBPwCfzVVcudTbmnGmmjayP3/42AzeWLeVj/zxWTZuU8InIiItksleaaz3lZHZYGb88IxpjBnUhwtve4ENurEqIhnIWbLn7h909+HuXubuo9z9Ondf5+7HuPtEdz/W3deHx7q7f87dd3P3Pd19Vq7iyqV4LxqNM51DJgzh2o/NYPHqLXz0uufYtF3zAomISKA3DWKWK9UVpVz1oX1Zu6Wer92huW5FpGO9rMFhbiVUkHHEpFqu+eh+vLyyjo9dr4RPREQCzTV7va2vQ5ZNG9mfi0+awr9eWsV1T74edTgikuf0i5tF6nweOGr3ofz2w/ux8O1NfPDaZ1i7pT7qkERECoKZjTazx8xsoZktMLMLw+2DzOwRM3s1fBwYdayd1dRL+7XnwrmHjuPEPXbhJw++zNOvrYs6HBHJY/rJzaJ42JpCnc/huKnD+MPHZrBk7RbO+v3TLN+4PeqQREQKQRPwFXefChwEfM7MpgIXAY+6+0Tg0XC9oCT7tatmr/vMjJ+/fy/GDe7DF259nhWbVMaKSHr6xc2i3joaZ1uOnDyUm847kDWb63n/757itTWah09EpD3uvsLdnw+fbwZeAkYCpwM3hIfdAJwRSYDd0NTc+iXiQIpETWUZ13x0Btsb4nz25uepb4pHHZKI5CH95GZRopcP0JLO/uMGcev5B1HflOD9v3+a2Us1XLSISCbMbBywD/AsMCxlSqKVwLA2XnO+mc0ys1lr1qzpmUAzlGhO9vSnR7ZMGFrNL96/Ny+8uZEf3rcw6nBEJA/pFzeL4r146oX2TBvZnzs+cwj9Kkv54B+e5Z4X3446JBGRvGZm1cDfgS+6e13qPg+GYEw7DKO7X+vuM9x9Rm1tbQ9EmrkmzbOXEyftOZxPHbErf3nmTW597s2owxGRPKNkL4s0Gmfbxg/py52fPZTpowZwwa0v8JtHX9WQ0SIiaZhZGUGid7O73xluXmVmw8P9w4HVUcXXVc01e5pnL+u+dvxkjphUy3fums9/F6+NOhwRySNK9rIonggeVbOX3qC+5dz0iQN4zz4jufyRV/jK7S+yo1F9DEREkszMgOuAl9z9lym77gHOCZ+fA9zd07F1l2r2cqc0VsJvPrQPu9b25dN/mc3i1ZujDklE8oSSvSxqmVQ94kDyWEVpjF+etTdfPm4Sd76wnPf97imWrd8WdVgiIvniUOCjwNFmNidcTgYuA44zs1eBY8P1gqKuDrnVr7KM687Zn4rSEv7vz7NYv7Uh6pBEJA8o2csid6fEgiGRpW1mxgXHTOS6c2awbP02TvnNkzy+qOBaJImIZJ27P+nu5u57ufv0cHnA3de5+zHuPtHdj3X3ghvtKh42f1GylzujB/Xh2o/NYGXdDj554yy2N6j1jEhvp2Qvi+IJV3+9TjhmyjDu/cJhDO9fybl//h+/+tcrzRPTi4hIcUk24yxTn72c2nfMQK44azrPv7mBz948m8ZkHxMR6ZWU7GVR3J0S3bHslLGD+/KPzx7KGdNH8qt/vcoH//CMJmAXESlCDWHSUaaJ9nLu3XsN50dnTOOxRWv4yu0v6kaqSC+mX9wsSiRcHc+7oKo86Md3+fv3ZsHyTZz0qye4b66mZxARKSaNTcmaPf3p0RM+fOBYvn7iZO558W2+d898jYAt0kvpFzeLGpoSVJTpI+0KM+N9+43igQvfxa611Xz+lhf46t9eZNP2xqhDExGRLGiMJ4iVmPrs9aDPHjmheQ6+H9y3UAmfSC+kzCSL6psSVJTqI+2OsYP78rdPH8wXjp7Anc+/xfFXzOSRhauiDktERLqpMZ6gXLV6Pe6iE3fn/w4dz5/++wbfvmt+83yHItI7lEYdQJSWrtvKbf9bttO2dDe9nLQb3+H5NzdQURrLUnS9V1mshK8cP5njpg7j63fM5ZM3zuKUvYZzyWl7MKS6IurwRESkCxriCQ3OEgEz4zunTKG8tITfz3yNhqYEl71vL9WwivQSvTrZW75xO9f95/V37kjz+9fWT2LrLnrHT92l23FJYK9RA7j3C4dxzczXuPLRxTy5eC1fOX4yHzpgjAopKQrxhNOUSNAU92BJJGhKeLDEw+fxlGPC7fGE05hw4okEjSmvjYevTSScuIePyW3uxBOEj8HS/NydeDzlNcljEztvaz73Tq+leVs83bnD5zd/4iBqa3SzpjdrjCcoV+uXSJgZ3zhxMhWlJfz60VfZsK2RKz84nT7lvfrPQJFeoVf/X37IbkN45dKTog5D2lEWK+HzR0/kxGm78O275vOdu+Zz8zNLueS0PTho18FRhyd5IpFwGuIJGuNB8tMYT9DQlGjZ1pS6P5Gy32lsStmW8tpgPXhty/7g/A1N8eb3aYynJGIpSVk8EewPHoPEbKfjEp62JUFPMoOYGSUlRmmJNT+PlRglZsRK2GlbLDymZX+4z2jeVllW0ryv+bWmfloSDNCiwVmiY2Z86bhJDK4u55J7FvCBa57huo/PYGhNZdShiUgO9epkTwrHhKE13PrJg3hw/kouvf8lzr72GU7ecxe+fNxkJgytjjq8Xs89SKZ2NCTY3hhne2OcHY1x6psS1CcfmxLUN8Wpb0x53pQI1+M77W+It97e9msbwoQq28yCmw0VsRLKSksoixllsRLKYyWUxUooKw3Wy0pKKC8toaokWI+VGGUxo7SkhNISozRmxEqC18dSjykJtpfGLDyu5fiykuCY0uR5Wh+Tsi/1nKWpSVlKQpaaqJWECVysxDCNHiw9qDGeULKXBz528DhGDqji87e8wBlX/ZerP7wv+4wZGHVYIpIjSvakYJgZJ+85nKN3H8o1M5dwzROv8dD8lZwxfSQXHDORcUP6Rh1i3muKJ9haH2dLQxNbdjSxpT5YtoaP2+qb2N4YJmwNTeFjgh2NcbYl1xsT7GiIs62xqXnf9sZ4txKu0hKjorSEirJY8FhaQkVpjIqy4HlVWYwBVWXhenBMeXhceWkJ5bEYZaXWkojFguSqvLSkJUELE7bUY8pLrfm1zdtiLYmZkiGR7NneGNcgZnnimCnD+NunD+ZTN83mrGue5psnT+Hjh4zTb55IEVKyJwWnsizGhcdO5CMHjeGaJ5Zw49NvcPeLb3Pa3iM477DxTBvZP+oQc6ahKcGm7Y1s2t7Apu2NbNwWLJu2N7JxeyObtgXbN+9oYnOYxCUTuS31TexoTGT8XlVlMarKY1SVxagsK6FPeSlVZTH6V5WxS78K+pSXUlkWC48recd6ZWmMymTylpKkVZTGmhO1irIguSrV3X6Rord5RxM1lfqzI19MG9mf+y84jK/+7UW+f+9C/rt4LZe+Z0+G9VOzTpFiol9dKViDqyv45slT+MRh4/n9zCX89X9v8o8XlnPA+EF8/JBxHDNlaN6PjppIOBu2NbB2SwNrt9SHS/h8c7C+bmsD67Y0sGFbA9sa4u2er19lKf37lFFTUUZ1ZSnD+lVSXVFK34pSqitiVFeU0bciRk1lsK1vRSk1zftL6VMeo095KRWlJZSoj5WIZNHm+ib6V5VFHYakGNCnnD98bAbX//cNfvbQyxz7y5l88+QpfGDG6JyUAdsb4izfuI3Vm4Oybs3meuq2N4atSOJsawi6ADQlErgHA58HfZtb+jjHUprHJ5uu7/xYEu5Ps725CXzQLD5WsnMT+WQz/FhJSdjUfucm86UpTfR3apbf/H4tze5VSyr5QsmeFLyh/Sr57qlTufDYifxt1jL+9N83+OzNz9OvspR37zWCM6aPYL+xA3u89sjdqdvexNubtvP2xu28vWkHb2/czoqN23l74w7e3rSdlZt20JSm+WNpiTG4upwh1RUMqa5gwtBqBvYpZ0BVGf37lNG/qowBfcqDx6oyBvQpo6ayTINgiEjWbGtoojEe/j55MA2Rt6zi3jIxkYf72Wl/8nnwujfXbeWQCUN68AokE2bGeYeN55jdh3LRnXO5+M553PT0Ui46aXcOn1TbpXO6O2+u38YLb25kzrKNLF69hSVrtvD2ph1pj09tPVJZVkJpSUnzaOdmhtEy+nlylN/kSMPJUY2b1+M7b49qWsESo1U/65T+2q2SxtJW/bxb+mO3JJHJvtfNA2uFA2CVpAyQFSw0d0NIfR4L10tSBuBKfU1JyjliJbS8JuW4WPh+qfsMI/yv+bsqsSBOg3BfsJ663cKdzdvDbclj23pekvK6dO9J8jwdxJJ8bbi5+XXh2dlpx077kuvW+pB3vL51vp+6nnpM8rPNFSV7UjT6V5XxiXftyscPGcd/X1vHXS8s564XlnPrc29SU1nKuyYO4bAJtUwfPYCJw6q7PVDAjsY4K8IE7u2N21uepyR1W1vVxJWWGLv0r2TEgCpmjB3I8AFVDKupYHCY1NXWBAle/6oy3RUUkUh95i/PM/OVNVk953A1Ecxb44b05dZPHsQ9L77Nzx9exMeuf44ZYwdy3mHjOW7qsHZvmG5raGLuW5t4/s0NPL90I3OWbWDtlgYA+pbHmDCshgN3HcyuQ/oyZnAfamsqqK2uoLamgn6VZTltSZJIkxQGIyS3mtYmfB5v9Tw5zU3LcS1T48RbPU+OvNxy/lbHpjlvumO3N8ZbnTcRTmMTJLvu4bQ44dQ3yalu3Am377wvqoRXMnPKXsO56kP75uz85lGP/d0NM2bM8FmzZkUdhuSxrfVNzHxlDTMXreHxV1azqq4egIrSEiYNq2H0oCpGDwwKnuqKUvpUlFIes2DuskQwBP/GbY2s39bAxq2NrNvawMq6oGZu/daGd7zfkOpyRgyoYkT/KoYPqGTkgCqG969ixIAgwRtSXaHaN5EuMrPZ7j4j6jgKRXfLyIfmr+StDdt2uoPdfIec8M55yjqtamCSd9CTry0pMY6bMoyBfcu7HJP0jPqmOLc++ybX/fd1lq3fztCaCo6ZMpSpI/ozqE85TYkEKzbtYOm6rby4bBOLVm1uHqRr1yF92WfMQPYdO4B9xwxk0rAalXsRc2+VCCbnXHXHU57vtK918hgmm633JWv4U5vbOkGimdyerPV3gte3NNENWwd4+u2+07aWVgKJMHdJ3Zbcnjyfp8aSch52ijd4j/AlzedMXU89hncc42m2pT8mXbqVPO+EodWcOG14O99gx9orH5XsSa/h7ry+divzlm9i3ltB4bR8w3be2ridhqb2By4pj5UwsG8ZA/uUs0v/Sob3r2LkgOBx+IBKRvSvYpf+lVSW5XcfQZFCpmSvc1RGSnfFE84jC1dx95zl/OfVtWypb9pp/8A+ZUwb2Z99Rg9g+pgB7DN6oJJ5kQi0Vz6qGaf0GmbGrrXV7FpbzenTRzZvTyR8p5ErG+MezJMWDsE/sE85fcpjalYpIiK9SqzEOHHaLpw4bRcSCWfV5h3UbW8iVmIM7Rc0wRSR/KZkT3q9khKjf1WZRokTERFpQ0mJBa1Zind2I5GipMmtREREREREipCSPRERERERkSKkZE9ERERERKQIKdkTEREREREpQkr2REREREREipCSPRERERERkSIUydQLZvYGsBmIA03uPsPMBgF/BcYBbwBnufuGKOITEREREREpdFHW7B3l7tNTZnu/CHjU3ScCj4brIiIiIiIi0gX51IzzdOCG8PkNwBnRhSIiIiIiIlLYokr2HPinmc02s/PDbcPcfUX4fCUwLJrQRERERERECl8kffaAw9x9uZkNBR4xs5dTd7q7m5mne2GYHCYTxC1mtqibsQwB1nbzHIVA11k8esM1gq6zmGTrGsdm4Ry9xuzZs9ea2dJunqY3/PuE3nGdveEaQddZTHrDNUJ2rrPN8tHc0+ZUPcbMLgG2AJ8EjnT3FWY2HHjc3Sf3wPvPSuk3WLR0ncWjN1wj6DqLSW+4xmLVW7673nCdveEaQddZTHrDNULur7PHm3GaWV8zq0k+B44H5gP3AOeEh50D3N3TsYmIiIiIiBSLKJpxDgP+YWbJ97/F3R8ys/8Bt5vZecBS4KwIYhMRERERESkKPZ7sufsSYO8029cBx/R0PMC1EbxnFHSdxaM3XCPoOotJb7jGYtVbvrvecJ294RpB11lMesM1Qo6vM/I+eyIiIiIiIpJ9+TTPnoiIiIiIiGRJr072zOxEM1tkZovN7KKo48kmM3vDzOaZ2RwzmxVuG2Rmj5jZq+HjwKjj7Awzu97MVpvZ/JRtaa/JAleG3+1cM9s3usg7p43rvMTMloff5xwzOzll38XhdS4ysxOiibpzzGy0mT1mZgvNbIGZXRhuL6rvs53rLLbvs9LMnjOzF8Pr/H64fbyZPRtez1/NrDzcXhGuLw73j4v0AiStYi0ji7F8BJWRRfabqjKySL7PvCgf3b1XLkAMeA3YFSgHXgSmRh1XFq/vDWBIq20/Ay4Kn18E/DTqODt5TYcD+wLzO7om4GTgQcCAg4Bno46/m9d5CfDVNMdODf/tVgDjw3/TsaivIYNrHA7sGz6vAV4Jr6Wovs92rrPYvk8DqsPnZcCz4fd0O3B2uP33wGfC558Ffh8+Pxv4a9TXoOUd32nRlpHFWD6GcauMfOexhfqbqjKySL7PfCgfe3PN3gHAYndf4u4NwG3A6RHHlGunAzeEz28AzogulM5z9yeA9a02t3VNpwM3euAZYIAF8zfmvTausy2nA7e5e727vw4sJvi3ndfcfYW7Px8+3wy8BIykyL7Pdq6zLYX6fbq7bwlXy8LFgaOBO8Ltrb/P5Pd8B3CMWTBEs+SN3lZGFnT5CCoj21Cov6kqI9MruO8zH8rH3pzsjQSWpay/Rfv/wAqNA/80s9lmdn64bZi7rwifrySYBqPQtXVNxfj9fj5snnF9ShOjgr/OsInCPgR3u4r2+2x1nVBk36eZxcxsDrAaeITgjutGd28KD0m9lubrDPdvAgb3aMDSkYL9t5iB3lI+QhH/pqZRVL+pSSojC//7jLp87M3JXrE7zN33BU4CPmdmh6fu9KB+uKiGYi3Ga0rxO2A3YDqwArg80miyxMyqgb8DX3T3utR9xfR9prnOovs+3T3u7tOBUQR3WnePNiKRNvW68hGK97pCRfebCiojKZLvM+rysTcne8uB0Snro8JtRcHdl4ePq4F/EPzjWpWs1g8fV0cXYda0dU1F9f26+6rwxyIB/IGWZgsFe51mVkbw436zu98Zbi667zPddRbj95nk7huBx4CDCZoSJedzTb2W5usM9/cH1vVspNKBgv+32JZeVD5CEf6mplOMv6kqI4vr+4ToysfenOz9D5gYjoZTTtAJ8p6IY8oKM+trZjXJ58DxwHyC6zsnPOwc4O5oIsyqtq7pHuBj4QhVBwGbUpo+FJxWbe/fQ/B9QnCdZ4ejN40HJgLP9XR8nRW2P78OeMndf5myq6i+z7auswi/z1ozGxA+rwKOI+h78RhwZnhY6+8z+T2fCfw7vEst+aMoy8heVj5Ckf2mtqUIf1NVRrYo6O8zL8rH1iO29KaFYPSiVwjazn4r6niyeF27EoxW9CKwIHltBG1+HwVeBf4FDIo61k5e160E1fmNBO2bz2vrmghGP7o6/G7nATOijr+b13lTeB1zwx+C4SnHfyu8zkXASVHHn+E1HkbQ/GQuMCdcTi6277Od6yy273Mv4IXweuYD3w2370pQEC8G/gZUhNsrw/XF4f5do74GLWm/16IrI4u1fAyvQWVk8fymqowsku8zH8pHC08sIiIiIiIiRaQ3N+MUEREREREpWkr2REREREREipCSPRERERERkSKkZE9ERERERKQIKdkTEREREREpQkr2REREREREipCSPZEsMLPBZjYnXFaa2fLw+RYz+20O3u/PZva6mX26k697IDm5ZzvH/Dy8hq92K0gRERFURopEqTTqAESKgbuvA6YDmNklwBZ3/0WO3/Zr7n5HZ17g7idncMzXzGxr18MSERFpoTJSJDqq2RPJITM70szuC59fYmY3mNl/zGypmb3XzH5mZvPM7CEzKwuP28/MZprZbDN72MyGZ/A+fzaz35nZM2a2JHzf683sJTP7c8pxb5jZEDMbF+77g5ktMLN/mllVzj4IERGRVlRGiuSekj2RnrUbcDRwGvAX4DF33xPYDrw7LMx+A5zp7vsB1wOXZnjugcDBwJeAe4ArgD2APc1seprjJwJXu/sewEbgfV28JhERkWxQGSmSZWrGKdKzHnT3RjObB8SAh8Lt84BxwGRgGvCImREesyLDc9/r7h6ee5W7zwMwswXhuee0Ov51d09umx0eIyIiEhWVkSJZpmRPpGfVA7h7wswa3d3D7QmC/x8NWODuB3f13OG56lO2J8/d1vEAcUBNVEREJEoqI0WyTM04RfLLIqDWzA4GMLMyM9sj4phERETygcpIkU5SsieSR9y9ATgT+KmZvUjQrOSQSIMSERHJAyojRTrPWmrIRaRQhKOH3dfZYaU7cf5L6JmhsUVERLJKZaRIC9XsiRSmTcAPOzthbCbM7OfARwDNIyQiIoVIZaRISDV7IiIiIiIiRUg1eyIiIiIiIkVIyZ6IiIiIiEgRUrInIiIiIiJShJTsiYiIiIiIFCEleyIiIiIiIkXo/wGKO2QXNpYQlgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "%matplotlib inline\n", + "from utils import simulate_and_plot_sorensen\n", + "\n", + "simulate_and_plot_sorensen('Fig71_Sorensen1985.xml')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Trajectories are plotted, which closely match\n", + "the original thesis figures.\n", + "\n", + "Note that this figure has a shift in the x-axis, as the simulation here\n", + "started at `time == 0`, but in the thesis started at approximately `time == -20`.\n", + "As there are no explicit time dependencies in the model, the figure is\n", + "otherwise similar." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# References\n", + "\n", + "[1] Sorensen, J. T. (1985).\n", + "\"A physiologic model of glucose metabolism in man and its use \n", + "to design and assess improved insulin therapies for diabetes.\"\n", + "https://dspace.mit.edu/handle/1721.1/15234\n", + "\n", + "[2] Panunzi, S., Pompa, M., Borri, A., Piemonte, V., & De Gaetano, A. (2020).\n", + "A revised Sorensen model: Simulating glycemic and insulinemic response to oral\n", + "and intra-venous glucose load. Plos one, 15(8), e0237215.\n", + "https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0237215" + ] + } + ], + "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.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/doc/release_notes.rst b/doc/release_notes.rst new file mode 100644 index 0000000..308e8ca --- /dev/null +++ b/doc/release_notes.rst @@ -0,0 +1,12 @@ +Release Notes +============= + + +0.2 series +.......... + + +0.2.0 (2021-02-03) +------------------ + +* Initial release on `pypi` diff --git a/setup.py b/setup.py index 970ad5d..9d58794 100644 --- a/setup.py +++ b/setup.py @@ -1,4 +1,5 @@ import setuptools +import os ENTRY_POINTS = { 'console_scripts': [ @@ -8,9 +9,13 @@ ] } +with open(os.path.join(os.path.dirname(__file__), + "yaml2sbml", "version.py")) as f: + version = f.read().split('\n')[0].split('=')[-1].strip(' ').strip('"') + setuptools.setup( name="yaml2sbml", - version="0.1.1", + version=version, author="Jakob Vanhoefer, Marta R. A. Matos", author_email="marta.ra.matos@gmail.com", description="A small package to convert ODEs specified in " @@ -26,10 +31,12 @@ "numpy>=1.19.4", "matplotlib>=3.1.0", "flake8>=3.7.2", - "nbmake>=0.1.0", ], + "nbmake>=0.1.0", + "scipy>=1.6.0"], extras_require={'examples': ["amici>=0.11.10", "numpy>=1.19.4", - "matplotlib>=3.1.0"]}, + "matplotlib>=3.1.0", + "scipy>=1.6.0"]}, python_requires='>=3.6', classifiers=[ "Programming Language :: Python :: 3.6+", diff --git a/tests/test_YamlModel.py b/tests/test_YamlModel.py index 43842b5..fb0f851 100644 --- a/tests/test_YamlModel.py +++ b/tests/test_YamlModel.py @@ -54,6 +54,9 @@ def test_time(self): self.assertTrue(model.is_set_time()) self.assertEqual(model.get_time(), time_var) + model.delete_time() + self.assertFalse(model.is_set_time()) + def test_parameter(self): """ Test all functionality regarding the 'parameters' keyword. @@ -68,13 +71,13 @@ def test_parameter(self): # test get_parameter_by_id self.assertIsInstance(model.get_parameter_by_id(parameter_id), dict) - # test over_write + # test overwrite with self.assertRaises(ValueError): model.add_parameter(parameter_id=parameter_id) # now over write parameter model.add_parameter(parameter_id=parameter_id, - over_write=True) + overwrite=True) # test get_parameter_ids self.assertListEqual(model.get_parameter_ids(), @@ -99,7 +102,7 @@ def test_ode(self): # test get_ode_by_id self.assertIsInstance(model.get_ode_by_id(state_id), dict) - # test over_write + # test overwrite with self.assertRaises(ValueError): model.add_ode(state_id, right_hand_side, initial_value) @@ -107,7 +110,7 @@ def test_ode(self): model.add_ode(state_id, right_hand_side, initial_value, - over_write=True) + overwrite=True) # test get_ode_ids self.assertListEqual(model.get_ode_ids(), @@ -131,12 +134,12 @@ def test_assignment(self): # test get_assignment_by_id self.assertIsInstance(model.get_assignment_by_id(assignment_id), dict) - # test over_write + # test overwrite with self.assertRaises(ValueError): model.add_assignment(assignment_id, formula) # now over write assignment - model.add_assignment(assignment_id, formula, over_write=True) + model.add_assignment(assignment_id, formula, overwrite=True) # test get_assignment_ids self.assertListEqual(model.get_assignment_ids(), @@ -161,12 +164,12 @@ def test_function(self): # test get_function_by_id self.assertIsInstance(model.get_function_by_id(function_id), dict) - # test over_write + # test overwrite with self.assertRaises(ValueError): model.add_function(function_id, arguments, formula) # now over write function - model.add_function(function_id, arguments, formula, over_write=True) + model.add_function(function_id, arguments, formula, overwrite=True) # test get_function_ids self.assertListEqual(model.get_function_ids(), @@ -195,7 +198,7 @@ def test_observable(self): # test get_observable_by_id self.assertIsInstance(model.get_observable_by_id(observable_id), dict) - # test over_write + # test overwrite with self.assertRaises(ValueError): model.add_observable(observable_id, observable_formula, @@ -206,7 +209,7 @@ def test_observable(self): model.add_observable(observable_id, observable_formula, noise_formula, - over_write=True, + overwrite=True, observable_name=observable_name) # test get_observable_ids @@ -234,7 +237,7 @@ def test_condition(self): # test get_condition_by_id self.assertIsInstance(model.get_condition_by_id(condition_id), dict) - # test over_write + # test overwrite with self.assertRaises(ValueError): model.add_condition(condition_id, condition_dict, @@ -244,7 +247,7 @@ def test_condition(self): model.add_condition(condition_id, condition_dict, condition_name=condition_name, - over_write=True) + overwrite=True) # test get_condition_ids self.assertListEqual(model.get_condition_ids(), diff --git a/yaml2sbml/YamlModel.py b/yaml2sbml/YamlModel.py index 319331c..f38b935 100644 --- a/yaml2sbml/YamlModel.py +++ b/yaml2sbml/YamlModel.py @@ -44,7 +44,7 @@ def load_from_yaml(yaml_file): # read in yaml_file with open(yaml_file, 'r') as f_in: yaml_contents = f_in.read() - new_model._yaml_model = yaml.full_load(yaml_contents) + new_model._yaml_model.update(yaml.full_load(yaml_contents)) # check, if the model is valid new_model.validate_model() @@ -53,14 +53,14 @@ def load_from_yaml(yaml_file): def write_to_yaml(self, yaml_dir: str, - over_write: bool = False): + overwrite: bool = False): """ Write model to yaml file given as directory yaml_dir. Arguments: yaml_dir: path/file, where the yaml should be written - over_write: + overwrite: Indicates, whether an existing yaml should be overwritten Returns: @@ -73,10 +73,10 @@ def write_to_yaml(self, raise ValueError('yaml_dir should contain path to the yaml ' 'and hence end with .yaml or .yml') - if (not over_write) and os.path.exists(yaml_dir): + if (not overwrite) and os.path.exists(yaml_dir): raise FileExistsError(f'Can not write yaml model. File {yaml_dir}' f' already exists. Consider to set ' - f'over_write=True.') + f'overwrite=True.') reduced_model_dict = self._get_reduced_model_dict() @@ -96,14 +96,14 @@ def write_to_yaml(self, def write_to_sbml(self, sbml_dir: str, - over_write: bool = False): + overwrite: bool = False): """ Writes the model as an SBML file to the directory given in sbml_dir. Arguments: sbml_dir: path/file, where the sbml should be written - over_write: + overwrite: Indicates, whether an existing yaml should be overwritten Returns: @@ -117,10 +117,10 @@ def write_to_sbml(self, raise ValueError('sbml_dir should contain path to the sbml ' 'and hence end with .xml or .sbml') - if (not over_write) and os.path.exists(sbml_dir): + if (not overwrite) and os.path.exists(sbml_dir): raise FileExistsError(f'Can not write SBML model. File {sbml_dir}' f' already exists. Consider to set ' - f'over_write=True.') + f'overwrite=True.') # generate SBML as string reduced_model_dict = self._get_reduced_model_dict() @@ -191,12 +191,15 @@ def _get_reduced_model_dict(self) -> dict: # functionalities regarding the time def is_set_time(self): - return bool(self._yaml_model['time']) + return 'variable' in self._yaml_model['time'].keys() def set_time(self, time_variable: str): self._yaml_model['time'] = {'variable': time_variable} + def delete_time(self): + self._yaml_model['time'] = {} + def get_time(self): if self.is_set_time(): return self._yaml_model['time']['variable'] @@ -206,7 +209,7 @@ def get_time(self): # functions adding a value def add_parameter(self, parameter_id: str, - over_write: bool = False, + overwrite: bool = False, nominal_value: float = None, parameter_name: str = None, parameter_scale: str = None, @@ -215,11 +218,11 @@ def add_parameter(self, estimate: float = None): """ Adds a parameter. Overwrites an existing parameter with the same id, - if over_write=True. + if overwrite=True. """ # if parameter exists: delete if overwrite if parameter_id in self.get_parameter_ids(): - if over_write: + if overwrite: self.delete_parameter(parameter_id) else: raise ValueError('Could not add parameter with id' @@ -240,14 +243,14 @@ def add_ode(self, state_id: str, right_hand_side: Union[float, str], initial_value: Union[float, str], - over_write: bool = False): + overwrite: bool = False): """ Adds a state/ODE. Overwrites an existing state/ODE with the same id, - if over_write=True. + if overwrite=True. """ - # if state exists: delete if over_write + # if state exists: delete if overwrite if state_id in self.get_ode_ids(): - if over_write: + if overwrite: self.delete_ode(state_id) else: raise ValueError(f'Could not add state/ODE with id {state_id}:' @@ -262,14 +265,14 @@ def add_ode(self, def add_assignment(self, assignment_id: str, formula: str, - over_write: bool = False): + overwrite: bool = False): """ Adds an assignment. Overwrites an existing assignment with the same id, - if over_write=True. + if overwrite=True. """ - # if assignment exists: delete if over_write + # if assignment exists: delete if overwrite if assignment_id in self.get_assignment_ids(): - if over_write: + if overwrite: self.delete_assignment(assignment_id) else: raise ValueError('Could not add assignment with id ' @@ -285,14 +288,14 @@ def add_function(self, function_id: str, arguments: str, formula: str, - over_write: bool = False): + overwrite: bool = False): """ Adds a function. Overwrites an existing function with the same id, - if over_write=True. + if overwrite=True. """ - # if function exists: delete if over_write + # if function exists: delete if overwrite if function_id in self.get_function_ids(): - if over_write: + if overwrite: self.delete_function(function_id) else: raise ValueError('Could not add function with id ' @@ -309,17 +312,17 @@ def add_observable(self, observable_id: str, observable_formula: str, noise_formula: str, - over_write: bool = False, + overwrite: bool = False, observable_name: str = None, observable_transformation: str = None, noise_distribution: str = None): """ Adds an observable. Overwrites an existing observable with the same id, - if over_write=True. + if overwrite=True. """ - # if observable exists: delete if over_write + # if observable exists: delete if overwrite if observable_id in self.get_observable_ids(): - if over_write: + if overwrite: self.delete_observable(observable_id) else: raise ValueError('Could not add observable with id ' @@ -338,12 +341,12 @@ def add_observable(self, def add_condition(self, condition_id: str, condition_dict: dict, - over_write: bool = False, + overwrite: bool = False, condition_name: str = None): """ Adds a condition. Overwrites an existing condition with the same id, - if over_write=True. + if overwrite=True. Arguments: condition_id: @@ -352,16 +355,16 @@ def add_condition(self, dict, of the form {: }. Corresponds to entries in the PEtab condition table. See details there. - over_write: + overwrite: bool, indicates if an existing condition should be overwritten condition_name: Condition name. Optional. Returns: """ - # if condition exists: delete if over_write + # if condition exists: delete if overwrite if condition_id in self.get_condition_ids(): - if over_write: + if overwrite: self.delete_condition(condition_id) else: raise ValueError('Could not add condition with id ' @@ -379,7 +382,7 @@ def _add_entry(self, block_key: str): """ Adds the entry in 'entry_dict' to the block indexed by 'block_key'. - If 'over_write=True', an existing value for that key is overwritten. + If 'overwrite=True', an existing value for that key is overwritten. Arguments: entry_dict: diff --git a/yaml2sbml/__init__.py b/yaml2sbml/__init__.py index d58a406..b11f69d 100644 --- a/yaml2sbml/__init__.py +++ b/yaml2sbml/__init__.py @@ -1,3 +1,9 @@ +"""Import of yaml2sbml's public API.""" + +# version +from .version import __version__ # noqa: F401 + +# API from yaml2sbml.yaml2sbml import yaml2sbml from yaml2sbml.yaml2PEtab import yaml2petab, validate_petab_tables from yaml2sbml.yaml_validation import validate_yaml diff --git a/yaml2sbml/version.py b/yaml2sbml/version.py new file mode 100644 index 0000000..7fd229a --- /dev/null +++ b/yaml2sbml/version.py @@ -0,0 +1 @@ +__version__ = '0.2.0'