{ "cells": [ { "metadata": {}, "cell_type": "markdown", "source": [ "# Treloar (1974): Michaelis-Menten Enzyme Kinetic Model\n", "\n", "Treloar (1974) considered the reaction rate in the presence of an enzyme treated and untreated with puromycin (a protein synthesis inhibitor). Here the concentration of the enzyme is in parts per million and the rate of the reaction is measured in counts per minute squared. The data used for this analysis can be obtained in R from `datasets::Puromycin`.\n", "\n", "Using this data, we will explore how puromycin modified the reaction by fitting a Michaelis-Menten Model, which is a model for enzyme-catalysed reactions for the transformation of a substrate into a single product. The equation is\n", "$$ V_i = \\frac{\\bar{V} \\times C_i}{K_m + C_i} $$\n", "where $V_i$ is the velocity of the reaction, $C_i$ is the substrate concentration, $\\bar{V}$ is the limiting rate (or maximum velocity of the reaction), and $K_m$ is the Michaelis constant which the concentration of substrate at which the reaction rate is half of the total velocity.\n", "\n", "## Setup" ], "id": "26dae5e696b59d5a" }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-13T18:36:36.453109900Z", "start_time": "2026-04-13T18:36:36.397100500Z" } }, "cell_type": "code", "source": [ "import numpy as np\n", "import scipy as sp\n", "import pandas as pd\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt\n", "import delicatessen as deli\n", "from delicatessen import MEstimator\n", "\n", "print(\"Versions\")\n", "print(\"NumPy: \", np.__version__)\n", "print(\"SciPy: \", sp.__version__)\n", "print(\"Pandas: \", pd.__version__)\n", "print(\"Matplotlib: \", mpl.__version__)\n", "print(\"Delicatessen:\", deli.__version__)" ], "id": "5798f7082e94e6f5", "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Versions\n", "NumPy: 2.3.5\n", "SciPy: 1.16.3\n", "Pandas: 2.3.3\n", "Matplotlib: 3.10.8\n", "Delicatessen: 4.2\n" ] } ], "execution_count": 10 }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-13T18:52:35.245798500Z", "start_time": "2026-04-13T18:52:35.228895Z" } }, "cell_type": "code", "source": "d = pd.read_csv(\"data/puromycin.csv\")", "id": "b5e3899f44e8d5c3", "outputs": [], "execution_count": 23 }, { "metadata": {}, "cell_type": "markdown", "source": "As with most data analyses should, we begin with a visualization of the corresponding data points. Here, we will plot the rate by the concentration stratified by puromycin treatment. The following code applies this process.", "id": "9f231409387d6ff2" }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-13T18:52:37.420287900Z", "start_time": "2026-04-13T18:52:37.404438800Z" } }, "cell_type": "code", "source": [ "conc = np.asarray(d['conc'])\n", "rate = np.asarray(d['rate'])\n", "puromycin = np.where(d['state'] == 'treated', 1, 0)" ], "id": "5d6d0818e6a49cbe", "outputs": [], "execution_count": 24 }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-13T18:55:22.646139400Z", "start_time": "2026-04-13T18:55:22.566086600Z" } }, "cell_type": "code", "source": [ "plt.plot(conc[puromycin == 1], rate[puromycin == 1], 'o', label='Treated')\n", "plt.plot(conc[puromycin == 0], rate[puromycin == 0], 'x', label='Untreated')\n", "plt.xlabel(\"Concentration (ppm)\")\n", "plt.ylabel(\"Velocity (counts/min/min)\")\n", "plt.legend()\n", "plt.tight_layout()" ], "id": "a3656910e70f4a3c", "outputs": [ { "data": { "text/plain": [ "
" ], "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVEJJREFUeJzt3QmcjfX+wPHvzGAszaDC2JJd2bIlCZWdqxHVJN203dt2iyhLRrIUbomi+t8Kl3TRgqgkXCpLyhKyVZaImbGNLYwx8/u/vr86586ZhTmcM2fOcz7v1+t5Pef5Pc95zjO/8zPz9VvDRMQIAAAAgl54oB8AAAAAvkFgBwAA4BAEdgAAAA5BYAcAAOAQBHYAAAAOQWAHAADgEAR2AAAADkFgBwAA4BAFAv0A+UW5cuXkxIkTgX4MAACALKKiomT//v1yIQR2fwZ1+/btu2BmAQAABEr58uUvGNwR2Im4a+o0w6i1AwAA+a22TiugchOjENhloBlGYAcAAIIVgycAAAAcgsAOAADAIQjsAAAAHII+dl4oWrSoXHnllRIWFua/bwR5whgjhw4dklOnTpHjAADHILDLBQ3kHnjgAbn55pv9/40gTy1btkymTJliAz0AAIIdgV0uaFDXqlUrmTVrlmzbtk3OnTvn/28GflWgQAGpVauW3HXXXfZ48uTJ5DgAIOgR2F1AsWLFbE2dBnWfffZZ3nwryBM7duyw+7i4OJk5cybNsgCAoMfgiQu44oor7F5r6uA8ru9V+04CABDsCOwuwDVQguZXZ3J9rwyIAQA4AYEdAACAQ9DHDvmSjlQtUaKE3H777YF+FAAAshUWHi5VGtaX6FJXyvGDh2Tnug1i0tMlkAjsHPjlX2jqjhdeeEGGDRvm888lGAMAhIq6rVtJ14FPS4mYMu60o4lJMnf0ONm05KuAPReBnQO//JiYGPdrHfE5fPhwqVmzpjvt5MmTHtdHRERIWlqaz58DAAAn0r/rvV4dpVUpHunFS5ey6VP7DgpYcEcfuzz68vXLzu7L1/O+lpSU5N6OHTtma/Bcxzp3mwZ2HTp0kDVr1khKSorcdNNNdvDAwIEDZefOnXbajx9++EG6d+/uvmd4eLi8++677vM6mvSpp55ynx86dKjcf//90rVrV/t5uuncf6pChQp2upjk5GQ5fPiwzJ07VypVquRx77Fjx9rzuhrEmDFjGMwAAMi3LXBdBz5tgzp9nfmcpscO6JPlXF4hsAvRL3/06NE2kLvmmmtk48aNMmjQILnvvvvk0Ucfldq1a8u4ceNk+vTp0rJlS3fw9dtvv8mdd94p1157ra0FfOmll+yxeuWVV2zwtmDBAltjqNvKlSvtRMALFy6UEydOSIsWLaR58+Y2sPziiy+kYMGC9r39+vWzQeGDDz5og8zLL7+cvnUAgHypSsP6tgUup7/dml6ybIy9LhBois2DLz8nGb/8HWvWS156/vnnZfHixfZ1oUKF5LnnnpM2bdrIt99+a9N27dplg6xHHnlEvv76azstiPbNc9m9e7c0a9bMrtzw4Ycfyu+//y6nT5+WyMhIWzPo0rNnTxsUPvzwwx4reRw9etRO/Lxo0SLp06ePjBo1SubMmWPPa3DZvn37PMwNAAByR/vK+/I6XyOwC9EvX5thXapVq2ZX2NAgKyMN+Nav/1/A+fjjj9tatauuukqKFCliz2uT7fnUr1/f3l9r7DIqXLiwVK1aVVavXi3lypWzexft76fPx9xyAID85vjBQz69ztcI7EL0y9caNpfLLrvM7jt37iz79u3zuE774LkGYWhzqzabrlq1ygZqzz77rDRt2vS8n6P3Xrt2ra25y+zgwYM++mkAAMgbOquFDoDUvvLZNcfqjBdHkw7Y6wKBwC6Ev3yXLVu2yJkzZ2xNnDa7Zkf7xmmfubfeesudpjVuGZ09e9aOsM1o3bp1Nig8cOBAllo7l/3799sA8ZtvvrHHeo9GjRrZ9wIAkJ+Y9HQ7q4UOgNTXGf++/zGNWZh8MmZ8wOazY/BEHnz5+iVn/oLzw5fvooMZtDZOB0zoAIoqVapIgwYN5B//+Ic9Vj///LM0btxY2rVrJ9WrV7eDJ5o0aeJxH+13V69ePalRo4ZdY1cHTrz//vt2pOsnn3xi++xdffXVdrTsa6+9JuXLl7fv09c6kCM2NtZOy/Lmm2/ayYkBAMiPNi35yk5pcuyAZ8uTVtYEcqoTRY1dHn35WeaxSzpgg7pAfvkZDRkyxDaN6uhYDex0cIPWmOnIV/Wvf/3LBns68lWnMpkxY4YNwDp27Oi+xzvvvGMHRGj/uKioKPv6q6++siNrdQqT2bNn23Rt7l2yZIkcP37cvk+nOilbtqxMnTpV0tPTZfLkyXYgRfHixQOWHwAAnI/+/f5x6Tf5buUJ+XN2vZDeoqKijNJ95nOVKlUy06ZNs/tL+Yyw8HBTtXED06BjW7vX40D/3Gy++37ZyAPKAGWAMkAZkADEKZk3auzyKnpOT8/zKU0AAEBooY8dAACAQ1BjBwA+pCPk8mOfGwChgcAOAHxE137OMlAqMcmOjs8vA6UAOBtNsQDgo6BO57XSeSsz0mNN1/MA4OjATucu++677+y0F7q+qE5xoXOgZaRrj06cONHOhaYT3H700UdSunRpj2sqVqwon376qV1NQe/zz3/+M8tEuQDgz+ZXranTQWmZJyP/49hI7IA+OS4aDiA4hYWHS9XGDaRBx7Z2nx/+jQe0KVYnqn3jjTfk+++/t5PZ6pxpX375pVx77bVy6tQpe41OmqtLXd15551y7NgxG+TpfGg62a3SBeY/++wzSUxMlBtvvNHOhzZt2jRJTU2VwYMHB/LHAxAitE9dxubXzPSXfcmyMfY6RscDzlA3n3a9CGhgl3FyW3X//ffbSXJ1OSldXio6Oloeeughueeee2Tp0qX2mgceeEC2bdtml6DSheN1JQQNBNu0aWOXrdqwYYOdbFcnxH3hhRdsgAcA/qQDJXx5HYDg6Hohduq4rF0vArn6RODrDDNwrTRw5MgRu9cAr1ChQrJ48WL3Ndu3b5dff/1VmjVrZo91v2nTJhvUuSxcuNDeq3bt2nn+MyD/0P8MaI0v4G86+tWX1wHIv8LyedeLfBPYhYWFyfjx42X58uWyefNmmxYTEyMpKSm2CTYj7Uen51zX6HHm865z2dFgUZe2yriFSlDTq1cvSU5OzvV9du3aJb1795a8QjCGYKRTmmgTTE7Tmmh6ckKivQ6AM7pehOUQuGXsehHSgZ32tatTp47cfffdfv8sXQ9VB2y4Nl27FBdP+zlqYA6EKg3ctF+NSFiW4O6P4zC7NjTz2QHBLzqfd73IF4HdhAkT5C9/+YvccsstHkGWDojQUbGZF4MvU6aMPee6Ro8zn3edy86oUaNs/z3XVr58efGXoUN7SHx8XLbnNF3PB8qUKVPsSOR+/frJ/v377chjHZyiA1lctWdXX321rUk1xtgtY61fly5dbO2q1qpeddVVtib05Zdflt9++01Onjwp3377rR0g43L55ZfLf/7zH3teRzBv3LjRI5DX57n55pulT58+7s+rVKmSPafN6p9//rkdGa3fqw6QueKKK9zvLVq0qEydOtWe15+lb9++eZiTwB8Lgmu/mmMHDnpkx9GkAwHtbwMgtLpehOeHoO7222+XW2+9VXbv3u1xbu3atXL27Flp3bq1O02nQ9E/9qtWrbLHuq9bt66UKvW/uaPatm1rm2+3bNmS7WfqPTUAyLj5S1paugwfcW+W4E6PNV3PB5IG01WrVrV7Ddh0AItuqlu3brJ37147GEWbtTM2bWsgNWDAAHn44Ydt0KV9HDUo1D6PGqzVq1dPPvzwQ/niiy+kWrVq9j2FCxe236mOctba2bffflvee+89adKkiT2vTb4rV6606a7P08/XwP6///2vrF+/Xho3biwdOnSwwfsHH3zgfh4NKDWIjI2NtQNqNEBs2LBhnucnQpsGbyPbd5M3H3hcpvd/3u5f7NCdoA5wkJ1B0PXCBGp74403THJysmnZsqUpU6aMeytcuLD7mjfffNPs3r3b3HzzzaZhw4ZmxYoVdnOdDw8PNxs3bjRffPGFqVevnmnXrp1JSkoyL774Yq6fIyoqyijdZz5XqVIlM23aNLu/2J8zPj7OpJv5dp/dsT+2pUuXmnHjxmVJ79Wrl81zfT1lyhSza9cum4eu87NmzTIzZsxwH+v53r17Z7mH0vx2pVWsWNGkpqaasmXLely7aNGi834X8+fPNy+//PJ5n3vw4MH2+82YVr58efsM1atXN8WKFTNnzpwxd9xxh/t8yZIlze+//55tHvj6+2UjDygDlAHKQGiVgbqtW5lXNqw0r2xYYcZuWuXe9FjT9bwvP+98cUrmLaDTnTz++ON2/9VXnk0UWmOkzWrq6aeflvT0dPn4449ts6yOeHW9T+k5bcZ96623bO2dNvHpe59//nnJL0aOnGX3WkM3OD5OIiMLyvNDprvTA0mbUjUPXRISEmwN6IVo86s2pbroe7QJ96effvK4Tr+zw4cPu/viPffcc3LXXXfZ5m9tutXzrjkLc1K/fn1bo5hdzarWNhYpUsTeR6e/cdGmYh1BDQCAv7peZJnHLumA7U8bsvPY5abDvQYQ//jHP+yWkz179tjmvfxMgzhXUJeSkur3oE4HhWTum6hKlCjhMco48zx/2q9NA7ALOX36tMfxZZddJufOnbNT1KSlpXmc0/526tlnn7XNrdqHTqeo0SBc++9pgHc+eu/58+fbpt/MNBB1NfUCAJBXNHj7cek3dvSrDpTQPnXa/BroQVIBDexCSXyGoE73euzP4E5rq7SvWWba7yxzrdr5aH/E3CzPpv3ftMZOl3vTKWuy07x5c/nkk0/k/fffdwf22mcyY1/I7D5v3bp10r17d9sHM3PQqHbs2GHfp5NWa588VwCr985cGwwAgK9oEJffVpMJ+OCJUOAaKKHNr0UKd7P77AZU+JI2TWtg89prr9lmUn2tzdo9evSQsWPH5vo+Gky1bNlSypUr5zEKNbOff/5Zpk+fbker6mAYHU2rgyJ0PeBOnTq5r9GBLTrAolatWvKvf/0ry4hm/TwN0HSAjH6eBn86FY6OqJ0xY4YdPFGlShUbtE6ePNnWLmrN36RJk+wACm2y1cEc//73vz2amAEACAUEdnkY1Llq6HTv7+BOJxbWgEwDKF25Q/ufad82XXNX+ynmlvZV1CBNa8V0OpTz0eXeNLDTwFFrDOfOnWuDO20qVyNHjrS1b/r5y5Yts9OW6DUZvfLKK7ZWTmvx9PN0GhVtbtXaPq3J07WEtRlXm3CPHj3qDt60mVeXodMmW/15tdZQR+ACABBqTKhv/hwVO3RojxxHv2q6ng/0zx/KG6NiA/8dsJEHlAHKAGVAnDEqNhQMGzYjx3P5YVQsAABwDppiAQAAHILADgAAwCEI7AAAAByCwA4AAMAhGDxxAboSg82oAmSVE7m+V9f3DFyqsPDwfDcTPYDQQbRyAa51TnU+OJ3LDc6i36u60Bx9QG7Ubd0q69qRiUkyd/S4gK4dCSB0ENhdgK5qoJPp6uS+atu2bXZNVAR/TZ0Gdfq96vd76tSpQD8SHBDU9Xp11J9TSf1P8dKlbLouGE5wB8DfCOxyYcqUKXYfF+e/JcAQGBrUub5f4FKaX7WmToM6fZ35nDbFxg7oYxcMp1kWgD8R2OWC9r/SdUlnzpwpV155pV2/FMH/nWrzKzV18AXtU5ex+TUzDe5Klo2x1+W3BcMBOAuBnRc0CHCtewoALjpQwpfXAcDFYroTALhEOvrVl9cBwMUisAOAS6RTmujo15z6z2l6ckKivQ4A/InADgAukQZuOqWJSFiW4O6P4zD5ZMx4Bk4A8DsCOwDwAZ3KRKc0OXbgoEf60aQDTHUCIM/o8M6Qn3I/KipKjh8/LtHR0XLixIm8y30AjsPKEwACGacwKhYAfEibXpnSBECg0BQLAADgEAR2AAAADkFgBwAA4BAEdgAAAA5BYAcAAOAQBHYAAAAOQWAHAADgEAR2AAAADkFgBwAA4BAEdgAAAA5BYAcAAOAQBHYAAAAOQWAHAADgEAR2AAAADhHQwK5FixYyb9482bdvnxhjJDY21uN8sWLFZMKECbJ37145deqUbN68WR555BGPayIjI2XixIly6NAhOXHihHz00UdSunTpPP5JAAAAQjyw08Btw4YN8sQTT2R7/tVXX5UOHTrIvffeK9dcc42MHz/eBnFdunRxXzNu3Dh7fOedd0qrVq2kXLlyMnv27Dz8KQAAAPIPkx82FRsb65G2adMmEx8f75G2Zs0aM2LECPs6OjrapKSkmO7du7vP16xZ096radOmuf7sqKgo+x7dBzof2MgDygBlgDJAGaAMUAbkIuOUfN3HbuXKlXLbbbfZWjh18803S40aNeTLL7+0x40aNZJChQrJ4sWL3e/Zvn27/Prrr9KsWbOAPTcAAEAgFJB87Mknn5S3337b9sFLTU2V9PR0+dvf/ibffPONPR8TEyMpKSly7Ngxj/clJSXZcznRYFD75rlERUX58acAAADIG+H5PbC74YYbbB86rZ3r16+fvPHGG9K6detLuu+gQYPk+PHj7k0DRwAAACfIl33sChcubPvPderUyeO6d955xyxYsMC+vuWWW+z7ihcv7nHN7t27TZ8+fXL8rEKFCtl2atdWrlw5+tjlgzLARh5QBigDlAHKAGVAnNnHrmDBgrbJVJtfM0pLS5Pw8D8ee+3atXL27FmPGjztg1epUiVZtWpVjvfW9+jUKBk3AACAYFcg0NOdVKtWzX1cuXJlqV+/vhw5csTOXbds2TJ5+eWX5fTp03ZAhE5nct9990nfvn3t9dqMOmnSJDstir5Hj3XeOx10sXr16gD+ZAAAAIERsGrPVq1amexMmTLFni9TpoyZPHmy+e2338ypU6fM1q1bzdNPP+1xj8jISDNx4kRz+PBhc/LkSfPxxx/b93nzHEx3QtU/Vf+UAcoAZYAyQBmQfJoH3sQpYX++CGk6KlZr+6Kjo2mWBQAAQRun5Ns+dgAAAPAOgR0AAIBDENgBAAA4BIEdAACAQxDYAQAAOASBHQAAgEMQ2AEAADgEgR0AAIBDENgBAAA4BIEdAACAQxDYAQAAOASBHQAAgEMQ2AEAADgEgR0AAIBDENgBAAA4BIEdAACAQxDYAQAAOASBHQAAgEMQ2AEAADgEgR0AAIBDENgBAAA4BIEdAACAQxDYAQAAOASBHQAAgEMQ2AEAADgEgR0AAIBDENgBAAA4BIEdAACAQxDYAQAAOESBi3lTxYoVpVKlSlK0aFE5ePCgbN68Wc6ePev7pwMAAIDvAzsN5B577DG5++67pUKFChIWFuY+p0HdN998I2+//bZ8/PHHYozJ/RMAAAAg75piX3vtNdmwYYNUrlxZ4uPj5dprr5XixYtLoUKFJCYmRjp16iTLly+X4cOHy8aNG6Vx48a+eToAAAD4tsbu999/lypVqsiRI0eynNOm2KVLl9pNA7v27dvbpto1a9bk/ikAAABwybQ9NeTbTaOiouT48eMSHR0tJ06coFgBAICgjFMYFQsAAOAQXgd2pUuXlmnTpsm+ffskNTVVzp0757EBAAAgSAK7f//739KwYUMZMWKE3HHHHdKtWzePzRstWrSQefPm2SBRR9LGxsZmuaZWrVryySefyNGjR+XkyZPy3Xff2T58LpGRkTJx4kQ5dOiQrZ786KOPbPAJAADgD0OH9pD4+Lhsz2m6ng+aeexuuukmG5DpKNlLVaxYMXufyZMny5w5c7Kc1wEbOtp20qRJMnToUNu+XLt2bTlz5oz7mnHjxknnzp3lzjvvlGPHjtkgb/bs2fY5EbrCwsOlSsP6El3qSjl+8JDsXLdBTHp6oB8LAOAAaWnpMnzEvfb1yJGzPII6TX9+yPQAPt0fgydyvW3evNlcd911Xr0nN5uKjY31SJsxY4aZNm1aju+Jjo42KSkppnv37u60mjVr2ns1bdo0158dFRVl36N7X/9cbHmfB3VbtzJDFs01Yzetcm96rOl8H5RJygBlgDJAGRAf5EF8fJxJN/PtPrtjX27exCleN8X26dNHRo8ebScs9iedAFlr4n766Sf54osvJCkpSb799luP5tpGjRrZufQWL17sTtu+fbv8+uuv0qxZsxzvre/RESYZNzhD3datpNero6R46VIe6Xqs6XoeAIBLpTV1WjOnNXSnz8x219RlrMELBK8Du1mzZsnNN98sO3bssE2jhw8f9th8RfvJacA1cOBAG9i1a9fONtdqM2vLli3tNTo5ckpKim2CzUiDQD2Xk0GDBtlnd23axw/OaH7tOvBp+58WfZ35nKbHDuiT5RwAABdDg7iUlFSJjCxo94EO6i6qj53W2OWF8D//+OrAifHjx9vX2h/vxhtvlEcffVS+/vrri773qFGj5NVXX3UfawBJcBf8tE9diZgyOZ7XgK5k2Rh73Y416/P02eB82lla+91k94td+91ERITLsGEzAvJsAPxD/227gjrd63GggzuvAzud6iQv6ChXnU5ly5YtHulbt251D4xITEy0o2J1ebOMtXZlypSx53Kia9vqBmfRgRK+vA5wUmdqAL6V8d+2/pt3HatABne5Cuy0Rss10/GF+qP5auUGDeq+//57qVmzpkd6jRo1bB86tXbtWhugtW7d2jbRus5r/79Vq1b55DkQPHT0qy+vA7zh+kWe8Rd75l/8AJwhPpt/29n9Dsi3gV1ycrKULVvWrgur88npnHPZDXbQ9AIFCng13Um1atXcx5UrV5b69evbNWn37t0rL7/8su3Tp82uuhZthw4dpEuXLraPn9L+cToVijar6nv0eMKECbJy5UpZvXp1rp8DzqBTmhxNTLIDJbLrR6fTnRxNOmCvA/wh4y/2wX820RDUAc4TERGe7b9t17Gez9drxepghRUrVkhaWpp74EJOvOn71qpVK1m2bFm2kyA/8MAD9rXudbBDhQoV7IhXnc9OJzV20abYsWPHSo8ePezrhQsXyuOPP24HUOQWa8U6b1Rs5gEUf8xhFyZT+w6STUu+Cugzwvl0hJyr302Rwt5N3A4AlxKn5CqwczoCO+cFdzo6NuNAiuSERPlkzHiCOuRZE42rMzU1dgDyMk7xevCE0pqxevXq2SlJXKNXXebPn38xtwR8Rmvkflz6DStPIM/l187UAEKLV7Mft2/f3iQlJZm0tLQs27lz54JyRm9Wngj8d8BGHgR7Gchp1nl/zkbPRh5QBkKjDER5sfKE1zV2Ojjhww8/lOHDh8uBAwf8E2oCQJDJz52pAYQOr/vY6XxxDRo0kJ07d4pT0McOAAA4IU7x+r+QH330kXu6EQAAAARxjV2RIkVsU6zOabdp0yY7kXDmptpgQ40dAAAIyVGxOl9cu3bt5MyZM7bmLuNkxfo6GAM7AACAkKyxS0hIkNdff11Gjx6d7QoUwYgaO+fRyYmrNKxv14XVJcR0tYk/JikGACC4+LXGrlChQnaZL6cEdQiNCYp1qbG5o8cxQTEAwNG8HjwxdepUiYuL88/TAD5aUkzXi81IjzVdzwMA4FRe19hFRERI//79pX379rJx48Ysgyf69evny+cDvGp+1Zq6zOvEus5pU2zsgD52VQqaZQEATuR1YFe3bl1Zv369fV2nTh2PczTPIpC0T13G5tfMNLgrWTbGXrdjzR9lGACAkA7sbr31Vv88CXCJdKCEL68DAMCxfex+/fVXO5VJmzZtbHMskN/o6FdfXgcAgGMDu7/+9a+SkpIib775phw6dEhmzpwp99xzjxQvXty/Twjkkk5poqNfc+o/p+nJCYn2OgAAQjqw+/rrr+WZZ56RGjVqSPPmzeWHH36QJ598UhITE2XJkiXSu3dvqVy5sn+fFjgPDdx0ShOdnjFzcPfHcZh8MmY8AycAAI7l9XQnasuWLXaC4mbNmtlgbsaMGdK6dWv58ccf7TJjnTp18v2TArmwaclXMrXvIDl24KBH+tGkAzZdzwMA4FRerzxxoXVkdRoUnRVZa/GCBStPOA8rTwAAnMKvK0+osLAwqVatmpQuXVrCM8wXptOdzJ0792JuCfiUNr0ypQkAINR4Hdg1bdpU/vOf/0ilSpVsgJeRBnYFClxUrAgAAIBL5HUU9n//93+yZs0a6dy5syQkJDApMQAAQLAGdtWrV5c77rhDduzY4Z8nAgAAQN6Mil29erXtXwcAAIAgr7HT1SfGjh0rMTExdmqT1NRUj/OaBgAAgCCY7iQtLS1Lmg6a0IEUwTp4gulOAABASE53wuoSAAAA+ZPXgd2ePXv88yQAEMSGDu0haWnpMnLkrCzn4uPjJCIiXIYNmxGQZwMQOnIV2HXp0kUWLFgg586ds6/PZ/78+b56NgAIGhrUDR9xr32dMbjToE7Tnx8yPYBPByCUmAttaWlpplSpUu7XOW3nzp274L3y4xYVFWWU7gP9LGzkAWUgeMtAfHycSTfz7T67YzbygDJAGRA/xyk+XSs2WDF4AoCvuGroUlJSJTKyoK2py655FgD8EacQ2BHYAfCx02dm26BOg7sihbuRvwDy76hY1bhxY7nlllukdOnSEh7uOcdxv379LuaWAOCYGjtXUKd7PabGDkBe8TqwGzRokIwcOVK2b98uSUlJHmvFZnwNAKEm40AJDeZcx4rgDkBe8aoDX2JiounVq5ejOoAyeCLw3wEbeRDsZSCngRIMoAj8d8NGHkgIxSle19ilp6fLihUr/BNiAkCQ0nnqshso4TrW8wCQF7yKGp999lkzbtw4n0SgLVq0MPPmzTP79u2zkWhsbGyO17711lv2mt69e3uklyxZ0kyfPt0cO3bMJCcnm3fffdcUK1bMb5EwG3lAGaAMUAYoA5QByoA4pcbulVdekc8++0x++eUX2bJli6Smpnqc7969e67vVaxYMdmwYYNMnjxZ5syZk+N1Xbt2lRtuuEH27duX5dz7778vZcuWlbZt20rBggVlypQp8vbbb0vPnj29/MkAAACCm9eB3euvv25HxC5dulQOHz58SQMmvvjiC7udT7ly5WTChAnSvn17G1BmVKtWLenYsaMdpbt27Vqb9uSTT8rnn38uzzzzjCQkJFz0s4WisPBwqdKwvkSXulKOHzwkO9dtEJOeHujHAgAA/grsevXqZWvlNHjyt7CwMHnvvffk5ZdftrWDmTVr1kySk5PdQZ1avHix7QfYtGlTmTt3brb3LVSokERGRnrMDxPq6rZuJV0HPi0lYsq4044mJsnc0eNk05KvAvpsAAAgd7zuzXvkyBHZsWOH5IUBAwbY9Wm1ljA7MTExcuDAAY+0tLQ0+4x67nxTtuhEf64tuybeUAvqer06SoqXLuWRrsearucBAIADA7sXXnhBhg0bJkWKFBF/atiwofTu3Vvuv/9+n9971KhRdvZm11a+fHkJ5eZXranTPpf6OvM5TY8d0CfLOQAA4ICm2KeeekqqVq1qJyfevXt3lsETjRo18smDtWjRwq5ssWfPnv89bIECMnbsWOnTp49UrlxZEhMT7TUZRUREyOWXX27P5eTs2bN2g9g+dRmbXzPTgK5k2Rh73Y4168kyAACcFNjl1G/N17RvnfaXy2jhwoU2XUe+qlWrVknJkiVt7d66dets2q233mqXOVu9enWePGew04ESvrwOAAAEQWCnNWS7du2S4cOH++zDdbqTatWqeXxG/fr1bR+5vXv32n1GWjuoNXE//fSTPd62bZssWLBA3nnnHXn00UftdCcTJ06UmTNnMiI2l3T0qy+vAwAAgZPrjlMbN26UTZs2yYsvvihNmjTxyYfrNCU//PCD3dS4cePsa2+CR52vTgO8JUuW2JG6y5cvl7///e8+eb5QoFOa6OjXnKY10fTkhER7HQAAyN/C/pyp+IJ0ehCdBDg2Nlb+8pe/2PnrPv30U5k3b54sWrRIUlJSJFjpdCc6OlYHUpw4cUJCdVRs5gEUfwR7YTK17yCmPAEAIAjilFwHdtnNIXfbbbfZ7aqrrrL94TTImz9/vhw6FFzNdqEe2OU0j53W1H0yZjxBHQAATg/sMtJ+chrgaW2eTgzct29fefPNNyVYENj9gZUnAADIf/I8sMtIpxrRTdeSDRYEdgAAwAlxitezzt53333SqVMn9/GYMWPssl4rVqywTbI6kjWYgjoAAACn8Dqwe+655+T06dP29Q033CBPPPGE9O/f3/ar01GtAAAACJIJiitWrOiukevatat8/PHHdh45rbFbtmyZP54RAAAA/qixO3nypFxxxRX2dbt27exUJ+rMmTN+Xz8WAAAAPqyx00Du3XfflfXr10uNGjXspMCqdu3adu1YAAAABEmNnfap0zVaS5UqJd27d3cv+9WoUSOZMWOGP54RAAAAueD1dCfax+63336zK09kd07XeA02THcCAABCcrqTXbt2yZVXXpklXeeu03NAIA0d2kPi4+OyPafpeh4AAKfyOrALC9NKvqwuu+wyO4ACCKS0tHQZPuLeLMGdHmu6ngcAQEJ98MTYsWPtXptghw8fLqdOnXKfi4iIsEuJ/fDDD/55SiCXRo6cZfcaxLmOXUHd80Omu88DABDSgV2DBg3cNXZ169aVs2fPus/p6w0bNsgrr7zin6cELjK4GxwfJ5GRBQnqAAAhwevBE5MnT5bevXtfsPNeMGHwhDOdPjPbBnUpKalSpHC3QD8OAAD5b/DEgw8+6KigDs6kza+uoE73OQ2oAAAgpCcoLlq0qAwcOFBat24tpUuXlvBwz9iwatWqvnw+wGuZ+9S5jhV97AAATuZ1YKerTrRq1Uree+89SUhIyHY+OyBQshsokd2ACgAAnMjrwK5jx47SuXNnWblypX+eCLgEERHh2Q6UcB3reQAAnMrrwRM7d+6UTp06ybZt28QpGDwBAABCcvDEkCFD7Dx2RYoUuZRnBAAAQKCbYvv162cHSCQlJcnu3bslNTXV43yjRo18+XwAAADwV2A3d+5cb98CAACA/NjHzonoYwcAAEKyjx0AAAAc0hSblpZ23rnrChTw+pYAAADwAa+jsNtvv93juGDBgtKgQQPp1auXDB061BfPBFy0oUN7SFpaeraTEOvkxTqP3bBhM8hhAIBjGV9sPXr0MHPnzvXJvfJ6i4qKMkr3gX4WtkvLg/j4OJNu5tt9btLZyAPKAGWAMkAZEGfFKb750MqVK5sTJ06EQoY5dgsLDzdVGzcwDTq2tXs9DvQzXcyWOYgjqAv8d8JGHlAGKAOUAcmTOMUnHeIKFy4sTz31lOzbt88Xt0MA1G3dSroOfFpKxJRxpx1NTJK5o8fJpiVfBdV3knFt2MHxcRIZWTDbZcYAAJBQn+7kyJEjHoMnwsLC7DDcU6dOyb333ivz58+XYBPq051oUNfr1VG2KISF/2+gtElPt0Vkat9BQRfcqdNnZtugLiUlVYoU7hboxwEAwO9xitc1dn369PE4Tk9Pl4MHD8rq1avl6NGj3j8tAkoDOa2pyxzUuc5pcBc7oI/8uPSbPwO94KADJVxBne71mBo7AIDTeR3YTZs2zT9PgoCo0rC+R/NrZhrclSwbY6/bsWa9BAMN4rQZ1tX86jpWBHcAACe7qD52xYsXl4ceekiuueYae7x582aZPHmyrSZEcIkudaVPr8tvQV3mPncZjwEAkFAP7Bo1aiQLFy6U06dPy3fffWfT+vbtK4MHD5Z27drJ+vXBUauDPxw/eMin1wWazlOX3UAJ17GeBwDAqbz+Kzdu3DiZN2+eXH311dK9e3e7Va5cWT799FMZP368V/dq0aKFvZeOptUBGbGxsR4rWIwePVo2btwoJ0+etNdMnTpVypYt63GPkiVLyvTp0+XYsWOSnJws7777rhQrVszbHytk7Vy3wY5+zan/nKYnJyTa64KBTj6cU42cpjM5MQDAybwO7Bo3bixjxoyxS4u56Ot//vOf9pw3NADbsGGDPPHEE1nOFS1aVBo2bCgjRoyw+27duknNmjVtIJjR+++/L7Vr15a2bdvKX/7yF2nZsqW8/fbb3v5YIUsDN53SREe/Zg7uXKNiPxkzPqgGTgAAEMq8miQvMTHRtG3bNkt6u3bt7LmLnXxPxcbGnveaxo0b2+sqVqxoj2vVqmWPGzVq5L6mffv2Ji0tzZQtW9YvE/85davbupUZsmiuGbtplXuL/3KOTQ/0s7GRB5QBygBlgDIQymUgyp8TFM+aNUsmTZokzzzzjKxcudKmNW/eXF5++WWZMcO/a3DqoA2dXsU1rUqzZs1s8+vatWvd1yxevNhe07RpU5k7d65fn8dJdJ46ndJER7/qQAntU6fNr9TUAQAQPLwO7DSg0/5wOu2J9oNTqamp8tZbb8nAgQPFXyIjI20TsAaPrsn5YmJi5MCBAx7XabOwTqKs53JSqFAhe7+ME//hj6bXYJnSBAAA+KCPnQZxOkmxDlq47rrr7Hb55ZfbkbFnz54Vf9AA8oMPPrCrXDz22GOXfL9BgwbZqVlcG0uhAQCAkAzsdDkLDep0upMff/zRbvpa0/xR8+UK6ipVqmQHSGRcSiMxMVFKly7tcX1ERIQNNPVcTkaNGmV/DtdWvnx5nz83AABAvg/sZs6cKXfffXeW9Lvuusue80dQV716dWnTpo1tYs1o1apVNqDUUbMut956q4SHh9slznKiNYsaIGbcAAAAnMCrkRmHDx+2o1Ezp9esWdMcOnTIq3sVK1bM1K9f326qT58+9rWOei1QoICZO3eu2bNnj6lXr54pU6aMeytYsKD7Hp9//rlZu3atadKkibnxxhvN9u3bzfvvv++30SZs5AFlgDJAGaAMUAYoA/l1VKzXgd3JkydNnTp1sqRr2u+//+7VvVq1amWyM2XKFFOpUiWTE32f6x4lS5a0gdzx48fN0aNHzaRJk2zASGDHPzp+8VIGKAOUAcoAZUBCLLALc0V3ufXf//7X9qt76qmnPNInTpwo9erVsxMEBxvtG6iDKLS/Hc2yAAAgWOMUr6c7iY+Pt3PF1a9fX5YsWWLTWrduLU2aNLFrxQIAACBIBk/opMQ6MfDevXvtgIkuXbrIL7/8Ymvrli9f7p+nhF8NHdpD4uPjsj2n6XoeAADkf17X2Cld3/Xee+/1/dMgINLS0mX4iD++z5EjZ3kEdZr+/JDpfDMAAASJC3bEK1q0qFed/Ly9PtAbo2LFxMfHmXQz3+6zO2YjDygDlAHKAGWAMiDBEKdc+Ib79+83AwYMMDExMee9rk2bNnb6kYEDBwbVl09gJx7B3Okzswnq8kG5ZCMPKAOUAcoAZUD8MSq2Ro0a8tJLL0nnzp1tM+yaNWtk//79cubMGTtB8LXXXmv73Z07d86u6vCvf/1L0tPTJVgwKvZ/Tp+ZLZGRBSUlJVWKFO4WwG8FAAD4ZVTsTz/9JHfccYdUrFhR7rzzTmnRooXceOONUqRIETl06JCsX79e/va3v8mCBQuCKqCDJ+1T5wrqdK/HGfvcAQCA/C/kq3ppiqWPHf8O+D1AGaAMUAYoAxKKK084cQv1wC6ngRIMoAj8d8NGHlAGKAOUAcpAlBdxykVNdwJniYgIt1OaZG52dR3reQAAkP95vaSYEzF4AgAAOCFOoSoGAADAIQjsAAAAQjWw27VrlwwZMsROfQIAAIAgDuzGjx8v3bp1k507d8qXX34pcXFxUqhQIf88HQAAAPwX2L322mvSoEEDuf7662Xr1q0yYcIESUhIsHtNBwAAQOBc0hw5BQoUME899ZQ5ffq0OXfunFm/fr154IEHgmrenVCfx46NPKAMUAYoA5QByoCE9jx2BQoUkNtvv10eeOABadu2rXz77bcyadIkqVChgl1Xtk2bNtKzZ0/fhqAAAADIkdeBnTa3ajDXo0cPuy7stGnT5Omnn5bt27e7r5kzZ458//333t4aAAAAeRnYacC2aNEieeyxx2Tu3Lly7ty5bEfOzpw581KeCzkICw+XKg3rS3SpK+X4wUOyc90GMenp5BcAAPB+5YmrrrpK9uzZ46isC5aVJ+q2biVdBz4tJWLKuNOOJibJ3NHjZNOSrwL6bAAAIAhXnli6dKlcfvnlWdKLFy8uO3bs8PZ28CKo6/XqKCleupRnvpcuZdP1PAAACG1eB3ZXX321REREZEmPjIyU8uXL++q5kKn5VWvqtHJVX2c+p+mxA/pkOQcAAEJLrvvYdenSxf26ffv2cuzYMfexBnqtW7eW3bt3+/4JYfvUZWx+zUwDupJlY+x1O9asJ8cAAAhRuQ7sdKCEMsbI1KlTPc6lpqbaoK5fv36+f0LYgRK+vA4AAIR4YOdqftWlxJo0aSKHDx/253MhAx396svrAACAM3k93UmVKlX88yTIkU5poqNfdaBEdv3odLqTo0kH7HUAACB05Sqwe/LJJ+Xtt9+WlJQU+/p8dM1Y+JYGbjqliY5+1dcZg7s/5rALk0/GjGc+OwAAQlyu5rHT5tfGjRvLkSNH7OucaP+7qlWrSrAJ5nnskhMSbVDHPHYAADiTN3GK1xMUO1GwBHaKlScAAAgtUV7EKV73sUNgadMrU5oAAIDseD2j7UcffST9+/fPkv7ss8/KBx984O3tAAAAEKjArmXLlvL5559nSV+wYIE9BwAAgCAJ7C677DI5e/ZslnSdpFjbfgEAABAkgd2mTZskLi4uS/rdd98tW7Zs8dVzAQAA4CIYb7a//OUv5uzZs+bf//63ue++++w2depUmxYbG+vVvVq0aGHmzZtn9u3bZ1R27x82bJjZv3+/OXXqlFm0aJGpVq2ax/mSJUua6dOnm2PHjpnk5GTz7rvvmmLFinn1HFFRUfbzde9tfrCRB5QBygBlgDJAGaAMiB/zwMs4xfsP6NSpk1m+fLk5efKkOXjwoFmyZIlp2bKl1/fp0KGDGTFihOnatWu2gV3//v1tsHbbbbeZunXrmrlz55odO3aYyMhI9zWff/65Wb9+vbn++utN8+bNzU8//WTef/99Ajv+kfGLljJAGaAMUAYoA44oA34P7PyxZRfYaU1dv3793MfR0dHm9OnTJi4uzh7XqlXLvq9Ro0bua9q3b2/S0tJM2bJl/ZVhbOQBZYAyQBmgDFAGKAMmPwZ2Xvexc2nYsKH07NnTbtddd534WuXKlaVs2bKyePFid5pOzrd69Wpp1qyZPdZ9cnKyrF271n2NXp+eni5Nmzb1+TMBAADkZ15PUFyqVCmZOXOm3HzzzXL06FGbVqJECVm6dKkdQHHo0CGfPFhMTIzdJyUleaTrseuc7g8cOOBxPi0tzS595romO4UKFZLIyEiPGZ0BAACCndc1dhMmTLCBUO3ateWKK66wW506dexUJ6+//roEg0GDBtnaP9e2b98+CWVDh/aQ+PisI52Vput5AADgwMCuQ4cO8vjjj8u2bdvcaVu3bpUnnnhCOnbs6LMHS0xMtPsyZf634L3r2HVO96VLl/Y4HxERIZdffrn7muyMGjXKBqKurXz58hLK0tLSZfiIe7MEd3qs6XoeAAAEB6868B0/ftzUr18/S/p1111npxzx9eCJvn37enQezG7wRMOGDd3XtG3b1tGDJ8LCw03Vxg1Mg45t7V6PfXHf+Pg4k27m2312x2zkAWWAMkAZoAxQBsR5o2J1ypFly5Z5BE7lypUzS5cuNbNnz/bqXjrfnAaJuqk+ffrY1xUrVnRPd3LkyBHTpUsXU6dOHTNnzpxspztZu3atadKkibnxxhvN9u3bHTvdSd3WrcyQRXPN2E2r3Jsea7ov7u8K5k6fmU1Qlw++bzbygDJAGaAMUAbE34FdhQoVzLp160xKSor55Zdf7KavNbgqX768V/dq1aqVyc6UKVM8JihOSEiwNXU6QXH16tWzTFCsgZzWJB49etRMmjTJkRMUa/D2yoaV5pUNKzwCOz3WdF8Fd66gTveB/pnZyAPKAGWAMkAZoAyIV3FK2J8vvNamTRupVauWu4/dkiVLJFjpYBAdRKH97U6cOCH5TVh4uMQvnC3FS5eyrzMz6elyNOmAvNihu319sVx96lJSUiUysqA8P2S6jBw56xKfHgAA5FWc4vV0Jxnni8s4xxz8p0rD+lIixnMQSUYa7JUsG2Ov27Fm/SUFda5gznWsCO4AAAgOuQrsnnzySa+mQ4FvRZe60qfXXSioU649wR0AAA4L7J5++ulc3cwYQ2DnB8cPHvLpdZlFRIRn2+zqOtbzAAAg/7voPnZOQh87AADghDjloqtiChYsKDVq1LATAsO/dEDE3NHjbByeeXDEH8dh8smY8Zc0cAIAAAQ/rwO7IkWKyLvvviunTp2SzZs3y1VXXWXTdTmxAQMG+OMZISKblnwlU/sOkmMHDnrkh46G1XQ9DwAA4NUcMePHjzfff/+9ad68uTlx4oSpXLmyTb/tttvs/HbBON9MMMxj5++VJ9jIA8oAZYAyQBmgDEjQxyleT3fStWtXiYuLk9WrV9vBEi5ae1e1alXCZD/T5taLndIEAAA4m9dNsaVKlZIDBw5kSS9WrJhHoAcAAIB8HtitWbNGOnfu7D52BXMPP/ywrFq1yrdPBwAAgFzzuin2ueeekwULFsi1114rBQoUkN69e9vXN954o7Rq1crb2wEAACCva+xq165t9ytWrJDrrrvOBnWbNm2Sdu3a2abZZs2aybp163z1XAAAAPDXBMVpaWny/fff26lOZs6cKSdPnhSnyO8TFAMAgNAV5Y8JirWZVUe+jh07VhISEmTKlCly0003+eJ5AQAA4AO5DuyWL18uDz30kJQtW1aefPJJqVy5snz11Veyfft26d+/v5QpU8YXzwMAAIBArBWr89Y98MAD8te//lViYmLkiy++kNjYWAk2NMUCAAAnxCmXFNipokWLSs+ePWXUqFFSokQJO6gi2BDYAQAAJ8QpFx2FtWjRQh588EHp3r27pKenywcffCCTJk262NsBAADgEnkV2Gn/uvvvv99u1apVk5UrV8pTTz1lg7pTp05d6rMAAAAgLwK7zz//XNq0aSOHDh2SadOmyeTJk+Wnn366lM8GAABAIAK71NRUueOOO+TTTz+1Ta8AAAAI0sAuGEe7AgAAhJJcz2OHwBo6tIfEx8dle07T9TwAAAhtBHZBIi0tXYaPuDdLcKfHmq7nAQAATKhvUVFRRuk+0M9yvi0+Ps6km/l2n90xG3lAGaAMUAYoA5QBCek4JfhmEw5hI0fOsnutoRscHyeRkQXl+SHT3ekAACC0XfLKE04QbCtPnD4z2wZ1KSmpUqRwt0A/DgAACPaVJ5C3dHCEqx+dK6jTvavPXUREuAwbNoOvBQCAEEZgF2SDJ5Sr+dU1cMKVBgAAQhujYgEAAByCGrsgoU2trlq5zIMnXOcBAEBoY/AEgycAAIBDBk9QzRNktF9ddoMnAAAACOyCiGuwhDa/6jQnus9uNQoAABC6TKhvwbDyRE6rTLD6ROC/GzbygDJAGaAMUAYkn8Qp+brGLjw8XIYPHy47d+6UU6dOyS+//CLx8fFZrhs2bJjs37/fXrNo0SKpVq2aOHXwROZVJvRY0xk8AQAAJD//T2PQoEHm4MGDplOnTqZSpUqme/fu5vjx4+bJJ590X9O/f3+TnJxsbrvtNlO3bl0zd+5cs2PHDhMZGemXSJiNPKAMUAYoA5QBygBlQPIwD7yMU/LvlzN//nzz7rvveqR99NFH5r333nMf79+/3/Tr1899HB0dbU6fPm3i4jybLAnsAv99spEHlAHKAGWAMkAZkNBtil25cqW0bt1aqlevbo/r1asnN910kyxYsMAeV65cWcqWLSuLFy92v0eHA69evVqaNWsWsOcGAAAIhHw9QfHo0aPtnC3btm2TtLQ0iYiIkMGDB8t//vMfez4mJsbuk5KSPN6nx65z2SlUqJBERkZ6zA8DAAAQ7PJ1jd1dd90lPXv2lHvuuUcaNmwovXr1kmeeeUbuu+++S7rvoEGDbM2ea9u3b5/PnhkAACCQ8m179549e8zjjz/ukTZ48GCzdetW+7py5cq2zbl+/foe1yxbtsyMHz8+x/sWKlTItlO7tnLlyjF4Ih9832zkAWWAMkAZoAxQBsS5feyKFi0q6enpHmnaJKvToKhdu3ZJQkKC7YeXsVm1adOmsmrVqhzve/bsWbskR8YNAAAg2OXrPnbz58+3fer27NkjmzdvlgYNGkjfvn1l8uTJ7mvGjx9v57b7+eefbaA3YsQIO6fd3LlzA/rsAAAAgZBvqz0vu+wyM27cOLN7925z6tQp88svv5gRI0aYggULelw3bNgwk5CQYKc5WbRokalevbpXn+PveezCwsNN1cYNTIOObe1ejwOdt2zkAWWAMkAZoAxQBiQo8sCbOCXszxchTZtvdRCFjsD1dbNs3datpOvAp6VETBl32tHEJJk7epxsWvKVTz8LAACEdpySr/vYBTsN6nq9OkqKly7lka7Hmq7nAQAAfIXAzk/CwsNtTZ1WiOrrzOc0PXZAnyznAAAALhZRhZ9UaVjfNr/mFLhpesmyMfY6AAAAXyCw85PoUlf69DoAAIALIbDzk+MHD/n0OgAAgAshsPOTnes22NGvJtMEyy6anpyQaK8DAADwBQI7P9HATac0EQnLEtz9cRwmn4wZn2PgBwAA4C0COz/Seeqm9h0kxw4c9Eg/mnTApjOPHQAA8CUmKPbzBMU2k8PD7ehXHSihfeq0+ZWaOgAA4Os4JV+vFesUGsTtWLM+0I8BAAAcjqZYAAAAhyCwAwAAcAgCOwAAAIcgsAMAAHAIAjsAAACHILADAABwCAI7AAAAhyCwAwAAcAgCOwAAAIcgsAMAAHAIAjs/Gjq0h8THx2V7TtP1PAAAgK8Q2PlRWlq6DB9xb5bgTo81Xc8DAAD4kgn1LSoqyijd+/re8fFxJt3Mt/vsjtnIA8oAZYAyQBmgDFAGxEdxSgGfhojIYuTIWXavNXSD4+MkMrKgPD9kujsdAADAV8L+jPBCWlRUlBw/flyio6PlxIkTfvmM02dm26AuJSVVihTu5pfPAAAAoR2n0McuD2ifOldQp/ucBlQAAABcCgI7P3MNlNDmV62p0312AyoAAAB8IeQ7bfpr8EROAyUYQEGZ498dZYAyQBmgDFAGxA+DPBk84UcREeHZDpRwHet5AAAAX2HwRB4NngAAALgYDJ4AAAAIQbQFAgAAOASBHQAAgEMQ2AEAADgEgR0AAIBDENgBAAA4RL4P7MqVKyfvvfeeHDp0SE6dOiUbN26URo0aeVwzbNgw2b9/vz2/aNEiqVatWsCeFwAAIFDydWBXokQJWbFihaSmpkrHjh3l2muvlX79+klycrL7mv79+8tTTz0ljz76qDRt2lR+//13WbhwoURGRgb02QEAAAIh3y5rMmrUKPP111+f95r9+/ebfv36uY+jo6PN6dOnTVxcnF+W6mAjDygDlAHKAGWAMkAZkDzMA2/ilHxdY3fbbbfJmjVr5IMPPpCkpCRZt26dPPzww+7zlStXlrJly8rixYvdabqCxOrVq6VZs2Y53rdQoUJ2FueMGwAAQLDL14FdlSpV5LHHHpOff/5Z2rdvL2+99Za8/vrrct9999nzMTExdq9BX0Z67DqXnUGDBtkA0LXt27fPzz8JAABAiAd24eHhtpZu8ODB8sMPP8g777xjN+1PdylGjRpl14V1beXLl/fZMwMAAARKvg7sEhISZMuWLR5pW7dulauuusq+TkxMtPsyZcp4XKPHrnPZOXv2rJw4ccJjAwAACHb5OrDTEbE1a9b0SKtRo4b8+uuv9vWuXbts8Ne6dWv3ee0vp6NjV61alefPCwAAEGj5dnRP48aNzdmzZ82gQYNM1apVTY8ePczJkyfNPffc476mf//+5siRI6ZLly6mTp06Zs6cOWbHjh0mMjKSUbH54DtkIw8oA5QBygBlgDIgeTYqNl8Hdrp17tzZbNy40U5hsmXLFvPwww9nuWbYsGEmISHBXrNo0SJTvXp1f2YYG3lAGaAMUAYoA5QByoDJqzzwJk4Jc0V3oUybb3V0rA6koL8dAAAI1jglX/exAwAAQO4R2AEAADgEgR0AAIBDENgBAAA4BIEdAACAQxDYAQAAOASBHQAAgEMQ2AEAADgEgR0AAIBDENgBAAA4BIGdHw0d2kPi4+OyPafpeh4AAMBXCOz8KC0tXYaPuDdLcKfHmq7nAQAAfMmE+hYVFWWU7n197/j4OJNu5tt9dsds5AFlgDJAGaAMUAYoA+KjOKWAT0NEZDFy5Cy71xq6wfFxEhlZUJ4fMt2dDgAA4Cthf0Z4IS0qKkqOHz8u0dHRcuLECb98xukzs21Ql5KSKkUKd/PLZwAAgNCOU+hjlwe0T50rqNN9TgMqAAAALgWBnZ+5Bkpo86vW1Ok+uwEVAAAAvhDynTb9NXgip4ESDKCgzPHvjjJAGaAMUAYoA+KHQZ4MnvCjiIhwW0P34ksfStXGDSS61JVy/OAhe+w6DwAA4CsEdn40bNgMqdu6lcQvnC0lYsq4048mJsnc0eNk05Kv/PnxAAAgxFBl5Eca1PV6dZQUL13KI12PNV3PAwAA+AqBnZ+EhYdL14FP2yZvfZ35nKbHDuiT5RwAAMDFIqrwkyoN69vm15wCN00vWTbGXgcAAOALBHZ+ogMlfHkdAADAhRDY+YmOfvXldQAAABdCYOcnO9dtsKNfTXp6tuc1PTkh0V4HAADgCwR2fqKBm05posvxZg7u/jgOk0/GjM8x8AMAAPAWgZ0f6Tx1U/sOkmMHDnqkH006YNOZxw4AAPhS2J9LUIS0qKgoOX78uERHR8uJEyd8fn8dAaujX10rT2jzKzV1AADA13EKK0/kAQ3idqxZnxcfBQAAQhhNsQAAAA5BYAcAAOAQBHYAAAAOQWAHAADgEAR2AAAADhFUgd2AAQPEGCPjxunEv3+IjIyUiRMnyqFDh+wQ4I8++khKly4d0OcEAAAIhKAJ7Bo3biyPPPKIbNjguQSXBnldunSRO++8U1q1aiXlypWT2bNnB+w5AQAAAsnk961YsWJm+/btpnXr1mbp0qVm3LhxNj06OtqkpKSY7t27u6+tWbOmUU2bNs31/aOioux7dB/on5WNPKAMUAYoA5QBygBlQC4yTgmKGrs33nhDPvvsM1myZIlHeqNGjaRQoUKyePFid9r27dvl119/lWbNmuV4P32PzuKccQMAAAh2+X7libi4OGnYsKE0adIky7mYmBhJSUmRY8eOeaQnJSXZczkZNGiQvPDCC355XgAAgEDJ1zV2FSpUkNdee0169uxpAzhfGTVqlF1vzbWVL1/eZ/cGAAAIlHxdY6dNrWXKlJF169a50woUKCAtW7aUf/zjH9K+fXs7KrZ48eIetXb6nsTExBzve/bsWbtlRpMsAADIb7yJT/J1YKd96urUqeORNmXKFNm2bZuMGTNG9u7dawO01q1bu0fC1qhRQypVqiSrVq3yOsP27dvn458AAADANzRe0andgjawO3nypGzevNkj7ffff5fDhw+70ydNmiSvvvqqHDlyRI4fPy4TJkyQlStXyurVq3P9Ofv377fNsRfKrPNltAaFl3IPkOf5HeWcfA8VlHXyPL+WS41XLiRfB3a58fTTT0t6erp8/PHHtll24cKF8vjjj3t9n9xk1oVoUEdgl7fI87xHngcG+U6ehwLKec5yG18EXWB3yy23eBzroArtb6cbAABAKMvXo2IBAACQewR2PqC1hjovni+nZAF5nt9Qzsn3UEFZJ8+DWdifS1AAAAAgyFFjBwAA4BAEdgAAAA5BYAcAAOAQBHa5pHPj7dq1S06fPi3ffvutNGnS5LzX33HHHbJ161Z7/caNG6Vjx46++L5Cijd5/vDDD8vXX39tJ6rWbdGiRRf8jnBpeZ5RXFycGGNkzpw5ZKuf81yXUJw4caKde/PMmTOyfft2fr/kQb737t3brnp06tQp2bNnj50YX+dORe60aNFC5s2bZyfz198VsbGxF3xPq1atZO3atbac//zzz9KrVy+yO5d08ATbefLgrrvuMmfOnDH333+/ueaaa8y//vUvc+TIEVOqVKlsr2/WrJlJTU01zzzzjKlVq5YZPny4SUlJMbVr1yaf/ZTn06dPN4899pipX7++qVmzppk8ebJJTk425cqVI8/9lOeurVKlSmbv3r3mq6++MnPmzCG//ZjnBQsWNN9995359NNPzY033mjzvmXLlqZevXrkux/zvUePHub06dN2r3netm1bs2/fPjN27FjyPZd53qFDBzNixAjTtWtXo2JjY897/dVXX21OnjxpXnnlFft39IknnrB/V9u1a0eeywXzm6DuQnnw7bffmgkTJriPw8LCzG+//WYGDBiQ7fUzZ8408+fP90hbtWqVeeuttyiQfsrzzFt4eLg5duyY+etf/0qe+zHPNZ+XL19uHnzwQTNlyhQCOz/n+SOPPGJ++eUXU6BAAcp1Hua7Xrt48WKPNA04vvnmG76Hi8j/3AR2o0ePNps2bfJImzFjhlmwYAF5Lhf4+5fbar1QVbBgQWnUqJEsXrzYnabVyHrcrFmzbN+j6RmvV7rUWU7X49LzPLOiRYva+2izLPyX588//7wcOHBAJk+eTDbnQZ7fdtttsmrVKnnjjTckMTFRNm3aJIMGDZLwcH6V+zPfdf1xfY+rubZy5crSqVMn+fzzz734xuEN/o5evKBbUiyvXXnllVKgQAFJSkrySNfjWrVqZfuemJiYbK/XdPgnzzMbM2aM7YOUOcCG7/K8efPm8tBDD8l1111HtuZRnlepUkVuvfVWef/9921gUa1aNXnzzTdtsDJ8+HC+Bz/l+4wZM+z7li9fLmFhYTa/33rrLRk1ahR57ic5/R3VPqaFCxe2/e6QPf6bB8cZMGCA3H333XL77bezGoifXHbZZfLee+/J3/72Nzl8+LC/PgaZaM2c1pD+/e9/l3Xr1skHH3wgL774ojz66KPklR9pJ/7nnnvODrho2LCh/d3SuXNniY+PJ9+R71BjdwGHDh2Sc+fOSZkyZTzS9VibQrKj6d5cj0vPc5d+/frJwIEDpU2bNraZCv7J86pVq9rmqPnz57vTXM2BqampUrNmTdm5cyfZ78M8VwkJCTZ/09PT3Wk6+r5s2bK2FknPwbdlXY0YMcL+R2bSpEn2+Mcff5RixYrJ22+/bQNrbcqFb+X0d/TYsWPU1l0ANXYXoL8odbh169at3WlaFa/H2tclO5qe8XrVtm3bHK/Hpee5evbZZ2XIkCHSoUMH+374L8912oc6derYZljXplMZLF261L7eu3cv2e/jPFcrVqywza96nUuNGjVstwOCOv+UdVef3YzBtEpLS3O/F77H39FLwwiTXAyN16Hu9913nx12/X//9392aHzp0qXt+alTp5qXXnrJY7qTs2fPmr59+9qpN4YOHcp0J37O8/79+9vpC7p162bKlCnj3ooVK0b59lOeZ94YFev/PK9QoYId7f3666+b6tWrm06dOpnExETz3HPPUc79mO/6O1zzPS4uzk7D0aZNG/Pzzz/bGRD4G5q7PNffxTodlW6qT58+9nXFihXtec1vzffM052MGTPG/h3V6ayY7kRyW94I7HKTBzqHzu7du23woEPlr7/+eve5pUuX2j9qGa+/4447zLZt2+z1OmS7Y8eO/ALwY57v2rXLZEd/IVPG/ZPnBHaByfMbbrjBTp+kgYlOfTJo0CA77Qzl3H/5HhERYZ5//nkbzJ06dcr8+uuvZuLEiaZ48eLkey7zu1WrVtn+jnbls+413zO/Z926dfY70rLeq1cv8lsunNdhrugOAAAAwY0+dgAAAA5BYAcAAOAQBHYAAAAOQWAHAADgEAR2AAAADkFgBwAA4BAEdgAAAA5BYAcAAOAQBHYAkEeGDh0q69evz5PPuvXWW2XLli0SHh74X/PXXHONXT9Y11wF4H8s0UEeUAYoAxdVBnQ9Xl23dMeOHXbZnz179ph58+aZW2+9NejKlC5L17t3b5/dT8XGxmZZL/Pyyy/Pk59nzZo15p577gl4vrq2Dz/80MTHxwf8OdjIA3F4HgT+v3IAglKlSpVk7dq1tmbo2Weflbp160qHDh1k6dKl8sYbb4gTae1XWJiuxHhxfv/9dzly5Ij4W/PmzaVq1ary8ccfS34xZcoUeeyxxyQiIiLQjwI4XsCjSzbygDIQfGXgs88+M3v37jVFixbNci7j4ugVK1Y0c+fONSdOnDDHjh0zs2bNMqVLl3afHzp0qFm/fr259957ba3Z0aNHzYwZM8xll13mviYsLMw8++yzdhF2rRnURdife+459/kKFSrY+yYnJ5vDhw/bz6tUqZL7vC4wPmfOHNOvXz+zf/9+c+jQIbuIe4ECBex5XXw8M03XRcf1nl26dDGbN282qamp9r6NGzc2X375pTl48KB93mXLlpkGDRq4P09/joz0OOPPmvHnGjJkiM1H/bn0XPv27d3n9bPU7bffbv773/+a33//3fzwww/mhhtuOO93M2HCBPPBBx94pLk+++9//7utWdV7aZ5FR0dnySdd8P7AgQP2+3rrrbdMwYIF3ddoXmkt7bhx48yRI0dMYmKiefjhh205mDx5sjl+/Lj9njp06ODx+XqP06dPB2VtLht5IMGVBwF/ADbygDIQZGWgZMmSJi0tzQwcOPC812ngsm7dOvP111+bhg0bmuuvv958//33NjjIGHBoMPDRRx+Z2rVrm5tuuskGXyNHjnRfM3r0aBuw3XfffaZKlSqmefPm5qGHHrLnNDjToOvdd981derUMbVq1TLTp083W7dudQckGrBoAPbmm2+amjVrms6dO5uTJ0/agMT182iwo02F2rysmyuwS0lJMcuXLzfNmjUzNWrUMEWKFDG33HKL6dmzp72Xft4777xjEhIS3MHolVdeaQMyfb/eS4+zC+z69OljnysuLs7eW39O/bxq1ap5BHZbtmwxnTp1MtWrV7cBmwaKEREROea7Bn/9+/f3SNPP1uB68eLFpn79+qZFixbmp59+snnlukbzSb8LDayvvfZa+5lJSUke34V+dxrwDR482D6n7jXg1UBf81PT3njjDRv0al5lfIZVq1bZ5wh0+WUjD8TZeRDwB2AjDygDQVYGmjRpYgOOrl27nve6Nm3a2D/6WqPmSrvmmmvse7XWS4/1D70GWRlr6MaMGWODAH2t6VrT4wrkMm8aYGkQlzFNAzqtkWrbtq07YNFgKDw83H2N1lZpAHO+PnYamKl69epdMIDVYEcDxvP1scsc2P32229m0KBBHtesXr3a1iZmDOwefPDBLPmnQWVOz6O1jFoDmvmz9bsoV66cO01rB8+dO+cOZDWftDYzY0D2yCOP2GBPf0ZXYKeBuuu85qkGjFOnTnWn6f1U06ZNPZ7h448/trV6gS6/bOSBODgP6GMHwGu57WfmGg3522+/udO2bt0qycnJ9pzL7t275eTJk+7jhIQEKV26tPsehQsXliVLlmT7GfXr15dq1arJiRMn3Jv2Y9P3aD8zl82bN0t6enq2n3E+KSkpsnHjRo80fd/bb78tP/30kxw9elSOHz8ul112mVx11VWSW1FRUVK+fHlZsWKFR7oeZ8wblfHz9bldz5CTIkWKyJkzZ7Kk79mzR/bv3+8+XrVqle3zVrNmTXfahg0b5PTp0x7X6LNWrFgx2+fRPD18+LBs2rTJnZaUlJTtM+p9GRkL+FcBP98fgAP9/PPP9g96rVq1fHK/1NRUj2Pt4uaapiNjkJEdDah0EEfPnj2znDt48GCuPuN8svv8qVOnyhVXXCG9e/eWX3/91QZ/GgAVKlRI/CHjs//R/e+PgRw5OXTokJQsWdIvz5L5eVzPlDktu2e8/PLLZceOHX57LgDMYwfgImiN28KFC+WJJ57ItgamePHi7to5rempUKGC+5zWRmnQoXOs5TaIPHXqlLRu3Trb8+vWrZPq1avLgQMHbNCQcdOatNw6e/Zsrkds6qjT119/XRYsWGB/Dg3sSpUq5dX9tGZx37599l6Z753bvMmJzpV37bXXZknXGsWyZcu6j2+44QZJS0uT7du3e9SAam1nxmv0WbXm9VLVqVMnz+bxA0IVTbEALooGdRq4fPfdd9KtWzfbHKo1eE8++aStvVKLFy+2TXTvv/++NGjQQJo0aSLTpk2TZcuW2Vq23NCgacyYMfLPf/5T/vrXv0qVKlWkadOm8uCDD9rzem+tofrkk0/kpptukquvvlpatWolr732mm3qzC1tDm7ZsqWUK1fO1sZdKNjUZ9Gf9/rrr7fPoMFn5vtpMFqmTBkpUaJEtvd5+eWXZcCAAXLXXXdJjRo1ZNSoUXLdddfZZ78UGnRrXmSmzbNa21ivXj17XoPTDz74wN10qrTWcdKkSTYA79ixowwbNkwmTpzorim8lOlx9PvQMgHAfwjsAFyUXbt2ScOGDe28dWPHjpUff/xRFi1aZIMZna/MJTY21tbwff311/aP+s6dOyUuLs6rzxoxYoT9jOHDh9tawFmzZrn7b2lTqQZk2n9s9uzZ9rwGJlrr5E2N3fPPP2+DQq3p00DxfB566CFb66i1he+9954NkLTGMKN+/fpJ27ZtbU1XTrVU+r5XX33V/mwaAOs8gLfddpv88ssvcik00Kxdu7YNFjPS+2oeff755/Lll1/avnKPP/64xzXal1EDV/2+NJ/nzZsnL7zwglyqHj162M/U7wmA/2gP6Ev7bxgAIN/RGs7o6Gh59NFH3cuZde3a1dacnm8SYa1dvP322336LAULFrTB4j333CMrV6706b0BeKLGDgAc6MUXX7QDOy5lpQxf0b59L730EkEdkAcYFQsADnTs2DHbZy8/cA1mAeB/NMUCAAA4BE2xAAAADkFgBwAA4BAEdgAAAA5BYAcAAOAQBHYAAAAOQWAHAADgEAR2AAAADkFgBwAA4BAEdgAAAOIM/w+yxJFlsuwxKQAAAABJRU5ErkJggg==" }, "metadata": {}, "output_type": "display_data", "jetTransient": { "display_id": null } } ], "execution_count": 27 }, { "metadata": {}, "cell_type": "markdown", "source": [ "We can see some visual differences, but we can do better and accompany this plot with some numeric estimates. We will do that by using the Michaelis-Menten model described previously.\n", "\n", "## Overall Michaelis-Menten Model\n", "\n", "To begin, we will fit an overall model to the data points. The following is the estimating functions for the Michaelis-Menten model, which is based on the score function of the corresponding model\n", "$$ \\psi(O_i; \\theta) =\n", "\\begin{bmatrix}\n", " (V_i - \\hat{V}_i) \\times \\left( \\frac{C_i}{K_m + C_i} \\right) \\\\\n", " (V_i - \\hat{V}_i) \\times \\left( \\frac{-\\bar{V_m} C_i}{(K_m + C_i)^2} \\right) \\\\\n", "\\end{bmatrix} $$\n", "where $\\hat{V_i}$ is the predicted value of the velocity from the model, the first function is for the max velocity ($\\bar{V}$) and the second function is for $K_m$. The following is code to implement this model.\n", "\n", "Note that the Michaelis-Menten model is not a very stable non-linear model. To help the root-finding process, we can select good starting values. For $\\bar{V}$, we use the maximum rate observed in the data. For $K_m$, we can look at the previous plot. We see that halfway to the maximum velocity (i.e., 100) occurs around 0.1 for the concentration. Thus we use this as the starting value.\n", "\n", "Finally, note that the data consists of only 23 observations. This is relatively few, thus the empirical sandwich variance estimator may not be a reliable estimator of the variance due to so few observations. To address this issue, we apply a finite-sample correction." ], "id": "e9d4845ca73481a7" }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-13T19:01:48.915710300Z", "start_time": "2026-04-13T19:01:48.903368Z" } }, "cell_type": "code", "source": [ "def psi_mmekm(theta):\n", " # Extracting parameters\n", " v_max, k_m = theta[0], theta[1]\n", "\n", " # Predicted velocity and residual calculations\n", " v_hat = v_max * conc / (k_m + conc)\n", " residual = rate - v_hat\n", "\n", " # Estimating functions for the parameters\n", " ef_vmax = (conc / (k_m + conc)) * residual\n", " ef_km = (-v_max * conc / (k_m + conc)**2) * residual\n", "\n", " # Returning stacked estimating functions\n", " return np.vstack([ef_vmax, ef_km])" ], "id": "596364d91f32ca21", "outputs": [], "execution_count": 28 }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-13T19:06:56.830725800Z", "start_time": "2026-04-13T19:06:56.822637Z" } }, "cell_type": "code", "source": [ "estr = MEstimator(psi_mmekm, init=[np.max(rate), 0.1], finite_correction='hc1')\n", "estr.estimate()" ], "id": "54d0947b8af2772a", "outputs": [], "execution_count": 29 }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-13T19:07:00.845011500Z", "start_time": "2026-04-13T19:07:00.811105900Z" } }, "cell_type": "code", "source": "estr.print_results()", "id": "981e16754cd0772b", "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==============================================================\n", " Estimation Method: M-estimator\n", "--------------------------------------------------------------\n", "No. Observations: 23 | No. Parameters: 2\n", "Solving algorithm: lm | Max Iterations: 5000\n", "Solving tolerance: 1e-09 | Allow P-Inverse: 1\n", "Derivative Method: approx | Deriv Approx: 1e-09\n", "Small N Correction: hc1 | Distribution: t-stat\n", "==============================================================\n", " Theta StdErr Z-score LCL UCL P-value S-value \n", "--------------------------------------------------------------\n", " 190.81 11.42 16.71 167.07 214.55 0.00 42.79 \n", " 0.06 0.01 5.03 0.04 0.09 0.00 14.13 \n", "==============================================================\n" ] } ], "execution_count": 30 }, { "metadata": {}, "cell_type": "markdown", "source": [ "The previous analysis ignores the puromycin treatment, which we saw in the plot looks like it has an impact of the kinetics of the interaction. Therefore, we repeat the analysis but fit the models stratified by puromycin treatment\n", "\n", "## Stratified Models\n", "\n", "With estimating equations, we can easily estimate stratified models by simply limiting the contributions to the estimating functions to those within the strata of interest. We use the previous code and this method to estimate the stratified models" ], "id": "bab5ae6503aade6b" }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-13T19:09:36.582750600Z", "start_time": "2026-04-13T19:09:36.572597200Z" } }, "cell_type": "code", "source": [ "def psi_strata(theta):\n", " # Subset model parameters\n", " theta_0 = theta[:2]\n", " theta_1 = theta[2:]\n", "\n", " # Model among those treated with puromycin\n", " ef_puro = psi_mmekm(theta_0) * (puromycin == 1)\n", "\n", " # Model among those untreated with puromycin\n", " ef_nopur = psi_mmekm(theta_1) * (puromycin == 0)\n", "\n", " # Stacking estimating functions\n", " return np.vstack([ef_puro, ef_nopur])" ], "id": "9737b24ff3012fd5", "outputs": [], "execution_count": 31 }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-13T19:09:57.106519Z", "start_time": "2026-04-13T19:09:57.094107200Z" } }, "cell_type": "code", "source": [ "estr = MEstimator(psi_strata, init=[np.max(rate), 0.1, np.max(rate), 0.1],\n", " finite_correction='hc1')\n", "estr.estimate()" ], "id": "a7fc95ae8c1aa5fe", "outputs": [], "execution_count": 33 }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-13T19:10:03.333271800Z", "start_time": "2026-04-13T19:10:03.303198900Z" } }, "cell_type": "code", "source": "estr.print_results()", "id": "7d71961b44d30f07", "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==============================================================\n", " Estimation Method: M-estimator\n", "--------------------------------------------------------------\n", "No. Observations: 23 | No. Parameters: 4\n", "Solving algorithm: lm | Max Iterations: 5000\n", "Solving tolerance: 1e-09 | Allow P-Inverse: 1\n", "Derivative Method: approx | Deriv Approx: 1e-09\n", "Small N Correction: hc1 | Distribution: t-stat\n", "==============================================================\n", " Theta StdErr Z-score LCL UCL P-value S-value \n", "--------------------------------------------------------------\n", " 212.68 5.78 36.80 200.59 224.78 0.00 61.12 \n", " 0.06 0.01 6.80 0.04 0.08 0.00 19.14 \n", " 160.28 6.66 24.07 146.34 174.22 0.00 49.73 \n", " 0.05 0.01 4.77 0.03 0.07 0.00 12.89 \n", "==============================================================\n" ] } ], "execution_count": 34 }, { "metadata": {}, "cell_type": "markdown", "source": "These results suggest that there is a difference between the maximum velocity. However, there is random variability that is not directly incorporated. To formalize this comparison, we can stack in an additional estimating function for the difference between the maximum velocity between the strata. As the empirical sandwich variance estimator automates the delta method, we can easily obtain statistical inference for the differences between the maximum velocity. The following code applies this process.", "id": "e6708d24fd1c68ce" }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-13T19:14:41.659744900Z", "start_time": "2026-04-13T19:14:41.643442900Z" } }, "cell_type": "code", "source": [ "def psi_test_diff(theta):\n", " # Subset model parameters\n", " model_params = theta[:4]\n", " diff = theta[4]\n", "\n", " # Stratified models\n", " ef_models = psi_strata(theta=model_params)\n", "\n", " # Estimating function for the difference\n", " v_max_puromycin = model_params[0]\n", " v_max_nopurom = model_params[2]\n", " ef_diff = (v_max_puromycin - v_max_nopurom) - diff * np.ones(rate.shape[0])\n", "\n", " # Stacking estimating functions\n", " return np.vstack([ef_models, ef_diff])" ], "id": "1181033f541ca684", "outputs": [], "execution_count": 35 }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-13T19:14:59.478759600Z", "start_time": "2026-04-13T19:14:59.468139500Z" } }, "cell_type": "code", "source": [ "estr = MEstimator(psi_test_diff, init=[np.max(rate), 0.1, np.max(rate), 0.1, 0.],\n", " finite_correction='hc1')\n", "estr.estimate()" ], "id": "8fd8aef7f2d287c5", "outputs": [], "execution_count": 37 }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-13T19:15:14.612521Z", "start_time": "2026-04-13T19:15:14.576839700Z" } }, "cell_type": "code", "source": "estr.print_results(subset=[-1])", "id": "e46d84557f7dad3c", "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==============================================================\n", " Estimation Method: M-estimator\n", "--------------------------------------------------------------\n", "No. Observations: 23 | No. Parameters: 5\n", "Solving algorithm: lm | Max Iterations: 5000\n", "Solving tolerance: 1e-09 | Allow P-Inverse: 1\n", "Derivative Method: approx | Deriv Approx: 1e-09\n", "Small N Correction: hc1 | Distribution: t-stat\n", "==============================================================\n", " Theta StdErr Z-score LCL UCL P-value S-value \n", "--------------------------------------------------------------\n", " 52.40 9.06 5.78 33.37 71.44 0.00 15.80 \n", "==============================================================\n" ] } ], "execution_count": 39 }, { "metadata": {}, "cell_type": "markdown", "source": [ "Here, we see the estimated difference between the maximum velocities. As evident by the 95% confidence intervals (33.4, 71.4), this difference is quite substantial. Thus, we may reasonably conclude that puromycin does indeed modify the reaction rate.\n", "\n", "The advantage of estimating equations (and `delicatessen`) in this context is that we are able to easily estimate stratified models, compare different parameters between the stratified models, and obtain valid statistical inference for the parameter comparisons. This is not the case when fitting stratified models generally.\n", "\n", "## References\n", "\n", "Michaelis L, & Menten ML. (1913). Die kinetik der invertinwirkung. *Biochem Z*, 49(333-369), 352.\n", "\n", "Srinivasan B. (2022). A guide to the Michaelis–Menten equation: steady state and beyond. *The FEBS journal*, 289(20), 6086-6098.\n", "\n", "Treloar MA. (1974). *Effects of Puromycin on Galactosyltansferase of Golgi Membranes*. Doctoral dissertation, University of Toronto" ], "id": "23464f48c462e2af" } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.6" } }, "nbformat": 4, "nbformat_minor": 5 }