technical-notes/2-4-test-case-logistic-model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Technical note: Sensitivities of a logistic model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this note we look at a logistic model, defined with initial conditions that depend on model parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The model\n",
    "\n",
    "State:\n",
    "\\begin{align}\n",
    "\\frac{dy}{dt} = h y (1 - y)\n",
    "\\end{align}\n",
    "\n",
    "Initial condition:\n",
    "\\begin{align}\n",
    "y(t_0, K, h) = \\left[1 + \\exp(-h(t_0 + K \\log 10)) \\right]^{-1}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Case 1: Time starts at 0\n",
    "\n",
    "In the first case we use\n",
    "\\begin{align}\n",
    "K &= -4 \\\\\n",
    "h &= 1\n",
    "\\end{align}\n",
    "and set $t_0 = 0$ to get\n",
    "\\begin{align}\n",
    "y_0(K, h) = \\left[1 + \\exp(-h K \\log 10) \\right]^{-1}\n",
    "\\end{align}\n",
    "\n",
    "We run for 20 time units."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Case 2: A negative starting time\n",
    "\n",
    "The second case is based on a drug-modelling example, and uses a negative starting time!\n",
    "\n",
    "\\begin{align}\n",
    "K &= -4 \\\\\n",
    "h &= 1\n",
    "\\end{align}\n",
    "\n",
    "Now t_0 = -4.605170185988091 and we will run the simulation for $5.1 - t_0$ time units."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the syntax defined from Myokit 1.34.0 onwards, we can write:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "import myokit\n",
    "\n",
    "model = myokit.parse_model(\n",
    "'''\n",
    "[[model]]\n",
    "c.y = 1 / (1 + exp(-c.h * c.K * log(10)))\n",
    "\n",
    "[engine]\n",
    "time = 0 bind time\n",
    "\n",
    "[c]\n",
    "h = 1\n",
    "K = -4\n",
    "dot(y) = h * y * (1 - y)\n",
    "''')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAAsTAAALEwEAmpwYAAAg+UlEQVR4nO3deXhc9X3v8fd3ZrRY8iIv8i5bGGxjm92K2QJxWA0kccLNQlYgNJTb0JumS8LTJDRtnnubpO29bRMShxIH0qSQhSQ4xMYmG5QQwEu8yQvIeJFsSZY3WbKsZWa+948Zw0QeedXR0cx8Xs+jR3PO+Y30ec6M5qNz5sw55u6IiEjhioQdQEREwqUiEBEpcCoCEZECpyIQESlwKgIRkQIXCzvA6RozZoxXV1eHHUNEJKesXr16n7tXZluWc0VQXV3NqlWrwo4hIpJTzGxnX8u0a0hEpMCpCERECpyKQESkwKkIREQKnIpARKTABVYEZrbYzPaa2cY+lpuZ/buZ1ZnZejO7LKgsIiLStyC3CB4FFpxg+S3A9PTXvcA3A8wiIiJ9COxzBO7+vJlVn2DIQuC7njoP9ktmVmFmE9y9MahMIjJ4xRNJuhNJunqSdMWT9CSOffkbt+NJJ55wEkkn4U4imSSecJLuJJKQcMc9tTzpkExPH7uddHB3/I1l4KTmAenpzPmpbM6byzNl3u+NeX+0PHP+8af8P+7nnWRATfUorp2R9TNhZyXMD5RNAuozphvS844rAjO7l9RWA1OmTBmQcCJyatyd9q44e9u6aGnr4lBHD4eP9tCa8XW4s4cjXXE6uhMc6U5wtDvOka4EHd1xuuKpF/5EUtdGycbszdv3ve3cvCsCyzIv6zPB3R8GHgaoqanRs0VkACWTzt62LnbuP8LO/R3s2H+EnQc6aG7tfOPF/2hPIut9IwbDhxQxYkgR5cUxykuijBhSxMQRpQwpjlJWHKU0FqWkKEJJLEpJLEJJLEJxLEpR1CiORYhFIhRFjaJohFjUiEUiRCNGNGLE0t8jZsSiRsQgYqnpaCT1EnNsecTAzDBLvfhE3ridmnFsvpmlv6eWWa9XqmPTln4Je3M6c8ybU388//h1ZNlmDrAwi6ABqMqYngzsCSmLiABd8QRbm9pY19DK+vpDbNjdyo79R+jsSb4xJhYxqkaVMX54KZdUVTB2WAljh5dQOayEyqGlVJSlXvhHlBUxtDhGJBL+C52cWJhFsAS438yeAC4HWvX+gMjA6uxJ8Ptt+/nt1r38of4QWxrb6E6kXvRHlRdz4aQRXH3eGKrHlFM9uozq0eVMGFFKLKojz/NJYEVgZo8D84ExZtYA/B1QBODui4ClwK1AHdAB3B1UFhF50962Tn6zZS+/3LyXF17bx9GeBGXFUS6eXMHdb63m4skVXDhpBJNHDhkUuy0keEEeNfTBkyx34JNB/X4ReVN3PMnT6/fwvZd2smbXIQAmVQzhfTWTuX7WOK6YNoqSWDTckBKanDsNtYicuv3tXXz/5V3850s7aWnr4ryxQ/mrG2dww+xxnD9+mP7jF0BFIJKXtja1sfiF7fx07W6640nmz6zk41efwzXTx+jFX46jIhDJI4c7e/jqM1v4/su7KIlFeN/cydx9dTXnjR0WdjQZxFQEInnA3Vm2sYkvLqllX3sXd11Vzaeun05FWXHY0SQHqAhEctzuQ0d58Gcb+dWWvcyZOJxH7qzhoskVYceSHKIiEMlR7s53freDf16xFXf4/G2zuOuqah3jL6dNRSCSg3oSSR54cgNPrmng7TMr+YeFF1A1qizsWJKjVAQiOaa9K86ffX8Nz7/awl/eOIM/v+48HQkkZ0VFIJJD9rZ18vFHV7K5sY2v/o+LeP9bqk5+J5GTUBGI5IhtLe3cufgV9rd388idNbx95tiwI0meUBGI5IA1uw5yz6MriUaMH/zpFToqSPqVikBkkFu98wAffuRlxg8v5bGPz2Pq6PKwI0meURGIDGJ1e9u557FVTBgxhB/ddyVjhpaEHUnykA44Fhmkmg93cufiV4hFjMfunqcSkMCoCEQGoc6eBJ/47ioOdXTznbvmMWW0PiMgwdGuIZFB6O9/Xsv6hla+9dG5XDh5RNhxJM9pi0BkkPnByl08/ko9fzb/XG6eMz7sOFIAVAQig8jmxsN84alarpk+hr+6aWbYcaRAqAhEBomueIJP/2Atw0uL+NcPXEI0otNGyMDQewQig8T/e/Y1tjS1sfiuGkbrCCEZQNoiEBkEVu44wLee38YH51Vx3fnjwo4jBUZFIBKyrniCz/x4PZNHDuHzt80OO44UIO0aEgnZt557ne37jvDdj8+jvER/kjLwtEUgEqKd+4/w9d/UcdtFE7h2RmXYcaRAqQhEQuLu/N2SWoqjER58h3YJSXhUBCIheXZTM7/d2sKnb5zBuOGlYceRAqYiEAlBdzzJPy7bwnljh3LnlVPDjiMFTkUgEoL/enkn2/cd4XO3ziIW1Z+hhEvPQJEB1t4V52u/ruOqc0czf6beIJbwqQhEBti3/3s7+49085kF52Om00hI+FQEIgPoUEc3j/z369w8ZxyXVFWEHUcECLgIzGyBmW01szozeyDL8hFm9nMzW2dmtWZ2d5B5RML27Re209YV59M3zgg7isgbAisCM4sCDwG3ALOBD5pZ74OlPwlscveLgfnAv5hZcVCZRMLU2tHDd363g1svHM/544eHHUfkDUFuEcwD6tz9dXfvBp4AFvYa48AwS+0oHQocAOIBZhIJzWO/30F7V5w/v2562FFE/kiQRTAJqM+YbkjPy/R1YBawB9gAfMrdk71/kJnda2arzGxVS0tLUHlFAnOkK87i323nhlljmTVBWwMyuARZBNkOh/Be0zcDa4GJwCXA183suL8Sd3/Y3WvcvaayUofbSe55YmU9hzp6+J/zzws7ishxgiyCBqAqY3oyqf/8M90N/MRT6oDtwPkBZhIZcD2JJItf2M686lHMnToy7DgixwmyCFYC083snPQbwHcAS3qN2QVcD2Bm44CZwOsBZhIZcEs3NLL70FHuvXZa2FFEsgrs5OfuHjez+4HlQBRY7O61ZnZfevki4EvAo2a2gdSupM+6+76gMomEYfEL2zm3spzrzh8bdhSRrAK9Coa7LwWW9pq3KOP2HuCmIDOIhGl9wyHWNbTy9++aQ0QXo5dBSp8sFgnQ917ayZCiKO+5rPcBcyKDh4pAJCCtHT0sWbeHd186ieGlRWHHEemTikAkID9e00BnT5KPXDEl7CgiJ6QiEAmAu/P9l3Zy2ZQK5kwcEXYckRNSEYgE4MVt+3l93xE+coWuPiaDn4pAJAD/+fudjCwr4tYLJ4QdReSkVAQi/ayptZNnNzfz/poqSouiYccROSkVgUg/e/yVXSTd+dDlepNYcoOKQKQfxRNJnli5i2unVzJ1dHnYcUROiYpApB89/1oLzYe7tDUgOUVFINKPfvqHPYwsK+LtM3VeIckdKgKRftLW2cOK2ibecdFEimP605LcoWerSD95ZmMTXfEk775U5xWS3KIiEOknP/3DbqaOLuOyKRVhRxE5LSoCkX7Q2HqU37++n3dfMgkznW5acouKQKQfPLV2D+7wHu0WkhykIhDpBz/7w24unVJB9Rh9dkByj4pA5Cxt2nOYLU1t3K6tAclRKgKRs/SztbuJRYzbLpoYdhSRM6IiEDkLiaTz1NrdzJ85llHlxWHHETkjKgKRs/D7bftpPtylN4klp6kIRM7CU2t3M6wkxvWzdEoJyV0qApEz1JNIsmJTMzfOHqfrDkhOUxGInKEXt+2n9WgPt+gqZJLjVAQiZ2jZhkaGlsS4ZvqYsKOInBUVgcgZiCeSLK9t4vpZY7VbSHKeikDkDLy8/QAHO3q45QLtFpLcpyIQOQO/2NBIWXGU+TMrw44ictZUBCKnKZF0lm9s4rrztVtI8oOKQOQ0vbL9APuPdHOrjhaSPBFoEZjZAjPbamZ1ZvZAH2Pmm9laM6s1s+eCzCPSH5ZtbKS0KKLdQpI3YkH9YDOLAg8BNwINwEozW+LumzLGVADfABa4+y4z08czZVBLJJ1l6d1CZcWB/fmIDKggtwjmAXXu/rq7dwNPAAt7jfkQ8BN33wXg7nsDzCNy1lbvPEhLW5eOFpK8EmQRTALqM6Yb0vMyzQBGmtlvzWy1mX0s2w8ys3vNbJWZrWppaQkorsjJLd3QSEkswnXna+NV8keQRZDtwq3eazoGzAVuA24GvmBmM467k/vD7l7j7jWVldovK+FIJp1nNjYxf2Yl5SXaLST5I8giaACqMqYnA3uyjHnG3Y+4+z7geeDiADOJnLH1u1tpOtzJggvGhx1FpF8FWQQrgelmdo6ZFQN3AEt6jXkKuMbMYmZWBlwObA4wk8gZW1HbRDRiXDdzXNhRRPpVYNu37h43s/uB5UAUWOzutWZ2X3r5InffbGbPAOuBJPCIu28MKpPI2VixqZkrpo1iRFlR2FFE+lWgOzrdfSmwtNe8Rb2m/wn4pyBziJytbS3t1O1t56NXTA07iki/0yeLRU7Bs5uaAbhxtnYLSf5REYicghW1TVw4aQQTK4aEHUWk36kIRE5i7+FO/lB/iJu0NSB5SkUgchK/3LwXd7hpjg4blfykIhA5iRWbmpgyqowZ44aGHUUkECoCkRNo74rzYt1+bpo9DrNsH5YXyX0qApETeG5rC92JpHYLSV5TEYicwIpNTYwqL2bu1JFhRxEJjIpApA/d8SS/3rKXG2aNJRrRbiHJXyoCkT68vH0/bZ1xbpqt3UKS31QEIn1YUdvMkKIob50+JuwoIoFSEYhk4e78cnMz10wfQ2lRNOw4IoFSEYhksXH3YRpbO3VuISkIJy0CM7vfzHTIhBSUZzc1ETG4fpaKQPLfqWwRjAdWmtkPzWyB6VM1UgBWbGqmZuooRpUXhx1FJHAnLQJ3/zwwHfg2cBfwmpn9HzM7N+BsIqGoP9DBlqY27RaSgnFK7xG4uwNN6a84MBL4sZl9NcBsIqHQtQek0Jz0CmVm9r+AO4F9wCPA37h7j5lFgNeAzwQbUWRgPbupmeljh1I9pjzsKCID4lQuVTkGuN3dd2bOdPekmb0jmFgi4TjU0c0rOw7wp9dOCzuKyIA5aRG4+4MnWLa5f+OIhOs3W/eSSLp2C0lB0ecIRDI8u6mZscNKuHhyRdhRRAaMikAkrSue4LmtLVw/axwRnWROCoiKQCTtxW37OdKd0LWJpeCoCETSnt3UTFlxlCvPHR12FJEBpSIQAZJJ51ebm3nbjEqdZE4KjopABFjbcIjmw13crEtSSgFSEYgAyzc2EYsYbz9/bNhRRAacikAKnruzvLaJq84bw4ghRWHHERlwKgIpeFub29ixv4Ob5+hoISlMKgIpeMs3NmOmk8xJ4VIRSMFbXtvE3CkjGTusNOwoIqEItAjSF7LZamZ1ZvbACca9xcwSZvbeIPOI9FZ/oINNjYdZcIGOFpLCFVgRmFkUeAi4BZgNfNDMZvcx7ivA8qCyiPRleW0TgA4blYIW5BbBPKDO3V93927gCWBhlnF/DjwJ7A0wi0hWz2xsYtaE4VSNKgs7ikhogiyCSUB9xnRDet4bzGwS8B5g0Yl+kJnda2arzGxVS0tLvweVwrS3rZPVuw6yQFsDUuCCLIJsp2/0XtP/CnzW3RMn+kHu/rC717h7TWVlZX/lkwL37KZm3OHmC3S0kBS2U7lC2ZlqAKoypicDe3qNqQGeMDNIXQntVjOLu/vPAswlAsDy2maqR5cxc9ywsKOIhCrIIlgJTDezc4DdwB3AhzIHuPs5x26b2aPA0yoBGQitR3t4sW4f97z1HNL/iIgUrMCKwN3jZnY/qaOBosBid681s/vSy0/4voBIkH61uZl40rlJ7w+IBLpFgLsvBZb2mpe1ANz9riCziGR6en0jkyqGcGlVRdhRREKnTxZLwTnU0c3zr7Zw20UTdElKEVQEUoCW1zYRTzrvvGhi2FFEBgUVgRScn69rZOroMi6YNDzsKCKDgopACsq+9i5e3LaPd140UUcLiaSpCKSgLNvQSNLhnRdrt5DIMSoCKSg/X9/I9LFDmTleHyITOUZFIAWjqbWTlTsOaGtApBcVgRSMX2xoxB3ecdGEsKOIDCoqAikYP1+3hzkThzOtcmjYUUQGFRWBFIT6Ax2srT+k3UIiWagIpCA8vb4RgNsu1G4hkd5UBJL33J2n1u7mkqoKXYlMJAsVgeS9jbsPs6WpjffOnRx2FJFBSUUgee+Hq+opiUX0/oBIH1QEktc6exI8tXY3Cy4Yz4ghRWHHERmUVASS11ZsauZwZ5z311SdfLBIgVIRSF770ap6JlUM4cppo8OOIjJoqQgkbzUc7OCFun28d+5kXYBG5ARUBJK3nly9G3d0tJDISagIJC8lk86P19Rz9Xmj9dkBkZNQEUheemn7fuoPHOV9c/UmscjJqAgkL/1oVQPDSmMsuGB82FFEBj0VgeSdw509LN3QyLsunkhpUTTsOCKDnopA8s4PV9bTFU/ygbdot5DIqVARSF5JJJ1HX9xBzdSRXDS5Iuw4IjlBRSB5ZUVtEw0Hj3LPW88JO4pIzlARSF5Z/LvtTB45hJvm6E1ikVOlIpC8sb7hECt3HOSuq6qJ6pPEIqdMRSB549svbGdoSUxvEoucJhWB5IWm1k5+sb6R99dUMaxUp5sWOR2BFoGZLTCzrWZWZ2YPZFn+YTNbn/560cwuDjKP5K/Hfr+DpDt3X10ddhSRnBNYEZhZFHgIuAWYDXzQzGb3GrYdeJu7XwR8CXg4qDySvzq64/zXy7u4afZ4nVdI5AwEuUUwD6hz99fdvRt4AliYOcDdX3T3g+nJlwCdJlJO25NrdtN6tId7rtEhoyJnIsgimATUZ0w3pOf15R5gWbYFZnavma0ys1UtLS39GFFyXVc8waLfbuOSqgpqpo4MO45ITgqyCLIdv+dZB5q9nVQRfDbbcnd/2N1r3L2msrKyHyNKrnv85V3sPnSUv75pJmY6ZFTkTMQC/NkNQOZxfJOBPb0HmdlFwCPALe6+P8A8kmeOdMX5+m/quHLaaK4+T5eiFDlTQW4RrASmm9k5ZlYM3AEsyRxgZlOAnwAfdfdXA8wieejRF3ewr72bv75ZWwMiZyOwLQJ3j5vZ/cByIAosdvdaM7svvXwR8CAwGvhG+g857u41QWWS/NHa0cOi57Zxw6yxzNV7AyJnJchdQ7j7UmBpr3mLMm7/CfAnQWaQ/PSt57fR1hnnr26aGXYUkZynTxZLztnb1sl3freDd108kVkThocdRyTnqQgk5zz06zq6E0n+8sYZYUcRyQsqAskpdXvb+a9XdvH+miqqx5SHHUckL6gIJGckk87f/mQDZcUxPn3j9LDjiOQNFYHkjMdX7uKVHQf43G2zGDusNOw4InlDRSA5oam1ky8v3cJV547mfXN1SiqR/qQikJzw4FMb6U4k+cfbL9SHx0T6mYpABr1nNjayYlMzf3njDKaO1hvEIv1NRSCDWuvRHr7wVC1zJg7nnrfqNNMiQQj0k8UiZyOZdB54cj0HjnTznbveQiyq/1tEgqC/LBm0/vWXr7JsYxMPLDifCyaNCDuOSN5SEcig9NTa3fz7r+v4QE0Vf6Irj4kESkUgg84fdh3kb368nnnnjOJL775ARwmJBExFIIPKnkNH+cR3VzN+eCmLPjKX4pieoiJB05vFMmgc6Ypzz2Or6OpJ8PgnLmdUeXHYkUQKgopABoVk0vn0D9aytekwi+96C9PHDQs7kkjB0Ha3hM7d+fIzW1ixqZkvvGM282eODTuSSEHRFoGEqjue5G9/uoEfr27gI1dM4a6rqsOOJFJwVAQSmoNHuvnT763mle0H+IsbpvOp66frCCGREKgIJBTbWtr5+KMraWzt5N/uuISFl0wKO5JIwVIRyIB7sW4f931vNUXRCI9/4nLmTh0VdiSRgqYikAH1xCu7+PzPNjKtspxv3/kWqkaVhR1JpOCpCGRA7NrfwZd+sYlnNzXzthmVfO1DlzK8tCjsWCKCikACdqQrzjd+W8d//Pd2YhHjMwtmcu8103QmUZFBREUggXB3lqzbwz8u3ULT4U7ec+kkHrjlfMYN17WGRQYbFYH0K3dnza6DfHnZFlbuOMiFk0bw0Icv1RvCIoOYikD6RXc8ydINjXznd9tZ19DK6PJivnz7hbyvpopoRJ8NEBnMVARyxtydtfWH+Pm6Rpas28O+9i6mVZbzpYVzuP2yyZSX6Oklkgv0lyqnpSueYNWOgzz/agvLNjax60AHxdEI82dW8qHLp3Dt9Eoi2gIQySkqAjmhju446+pbWbPrIKt2HOCl1w9wtCdBUdS4Ytpo7r/uPG6eM54RQ3QoqEiuUhEIAPFEkj2HOtm2r52tTW1saTzMlqY2XtvbTiLpAEyrLOd9NZN524xKrpg2Wrt+RPJEoH/JZrYA+DcgCjzi7l/utdzSy28FOoC73H1NkJkKUXc8yYEj3exr76KlvYuWti4aD3XS2HqUPa2d1B/ooP5AB/H0Cz7AhBGlnD9+GDfMGsfcqSO5pKqCkbpQjEheCqwIzCwKPATcCDQAK81sibtvyhh2CzA9/XU58M3097yWTDoJdxJJJ550EgknnkwSTzo9iSTx9HRXPElPwumOJ+mOJ+mKJ+iKJ+nsSXC0J8HR7tTXke4ER7vjtHXFaeuM094Z53BnD4c6ejjU0c2R7kTWHGOGljCxopRZE4ZxywXjqR5dTvWYcmaOG8aIMu3qESkUQW4RzAPq3P11ADN7AlgIZBbBQuC77u7AS2ZWYWYT3L2xv8M892oLX3p6E6lfleK9bjgct9wdkul5xxa5+xvLHCfp6dvuJD01nXR/476J5JvzExn/dfeX4liEsuIoQ0tiDC2JMby0iHHDS5k5fhgVQ4qpKCtiVHkxY4aWUDmshMqhJYwbUUJJLNrvWUQk9wRZBJOA+ozpBo7/bz/bmEnAHxWBmd0L3AswZcqUMwoztCTGjHFDMQwyDmo5dvPYefANyDwlfsQsNcaOjTfM3hwXMUuPN6KR1PJo5NiY1LzUmNTtqBmRiBE1Ixo1iiIRIhGjKJq6X1EkQixqxKIRiqMRimNGcTRKcSxCSSxCSVGE0liUIcVRSouilBVHKdLpGkTkLARZBNmOIez97/CpjMHdHwYeBqipqTmjf6nnTh3J3Klzz+SuIiJ5Lch/JRuAqozpycCeMxgjIiIBCrIIVgLTzewcMysG7gCW9BqzBPiYpVwBtAbx/oCIiPQtsF1D7h43s/uB5aQOH13s7rVmdl96+SJgKalDR+tIHT56d1B5REQku0A/R+DuS0m92GfOW5Rx24FPBplBREROTIebiIgUOBWBiEiBUxGIiBQ4FYGISIGzzFMq5AIzawF2nuHdxwD7+jFOfxmsuWDwZlOu06Ncpycfc01198psC3KuCM6Gma1y95qwc/Q2WHPB4M2mXKdHuU5PoeXSriERkQKnIhARKXCFVgQPhx2gD4M1FwzebMp1epTr9BRUroJ6j0BERI5XaFsEIiLSi4pARKTA5WURmNkCM9tqZnVm9kCW5WZm/55evt7MLhuATFVm9hsz22xmtWb2qSxj5ptZq5mtTX89GHSu9O/dYWYb0r9zVZblYayvmRnrYa2ZHTazv+g1ZsDWl5ktNrO9ZrYxY94oM3vWzF5Lfx/Zx31P+HwMINc/mdmW9GP1UzOr6OO+J3zcA8j1RTPbnfF43drHfQd6ff0gI9MOM1vbx30DWV99vTYM6PPL3fPqi9Qpr7cB04BiYB0wu9eYW4FlpK6QdgXw8gDkmgBclr49DHg1S675wNMhrLMdwJgTLB/w9ZXlMW0i9YGYUNYXcC1wGbAxY95XgQfStx8AvnImz8cAct0ExNK3v5It16k87gHk+iLw16fwWA/o+uq1/F+ABwdyffX12jCQz6983CKYB9S5++vu3g08ASzsNWYh8F1PeQmoMLMJQYZy90Z3X5O+3QZsJnV95lww4Ourl+uBbe5+pp8oP2vu/jxwoNfshcBj6duPAe/OctdTeT72ay53X+Hu8fTkS6Su/Deg+lhfp2LA19cxZmbA+4HH++v3nWKmvl4bBuz5lY9FMAmoz5hu4PgX3FMZExgzqwYuBV7OsvhKM1tnZsvMbM4ARXJghZmtNrN7sywPdX2RurpdX3+cYayvY8Z5+op66e9js4wJe919nNTWXDYne9yDcH96l9XiPnZ1hLm+rgGa3f21PpYHvr56vTYM2PMrH4vAsszrfYzsqYwJhJkNBZ4E/sLdD/davIbU7o+Lga8BPxuITMDV7n4ZcAvwSTO7ttfyMNdXMfAu4EdZFoe1vk5HmOvuc0Ac+H4fQ072uPe3bwLnApcAjaR2w/QW2voCPsiJtwYCXV8neW3o825Z5p32+srHImgAqjKmJwN7zmBMvzOzIlIP9Pfd/Se9l7v7YXdvT99eChSZ2Zigc7n7nvT3vcBPSW1uZgplfaXdAqxx9+beC8JaXxmaj+0iS3/fm2VMWM+1O4F3AB/29M7k3k7hce9X7t7s7gl3TwL/0cfvC2t9xYDbgR/0NSbI9dXHa8OAPb/ysQhWAtPN7Jz0f5N3AEt6jVkCfCx9NMwVQOuxTbCgpPc/fhvY7O7/t48x49PjMLN5pB6f/QHnKjezYcduk3qjcWOvYQO+vjL0+V9aGOurlyXAnenbdwJPZRlzKs/HfmVmC4DPAu9y944+xpzK497fuTLfV3pPH79vwNdX2g3AFndvyLYwyPV1gteGgXt+9fc74IPhi9RRLq+Sejf9c+l59wH3pW8b8FB6+QagZgAyvZXUJtt6YG3669Zeue4Hakm98/8ScNUA5JqW/n3r0r97UKyv9O8tI/XCPiJjXijri1QZNQI9pP4LuwcYDfwKeC39fVR67ERg6YmejwHnqiO13/jY82xR71x9Pe4B5/rP9PNnPakXqwmDYX2l5z967HmVMXZA1tcJXhsG7PmlU0yIiBS4fNw1JCIip0FFICJS4FQEIiIFTkUgIlLgVAQiIgVORSDSBzMbnXFWyqaMM2e2m9k3ws4n0l90+KjIKTCzLwLt7v7PYWcR6W/aIhA5TZa6DsLT6dtfNLPHzGxF+nz1t5vZV9PnrX8mfeoAzGyumT2XPmHZ8gE+e6vICakIRM7eucBtpE7/+z3gN+5+IXAUuC1dBl8D3uvuc4HFwP8OK6xIb7GwA4jkgWXu3mNmG0hdKOSZ9PwNQDUwE7gAeDZ9aqQoqdMciAwKKgKRs9cF4O5JM+vxN994S5L6GzOg1t2vDCugyIlo15BI8LYClWZ2JaROORzCRXRE+qQiEAmYpy4h+F7gK2a2jtTZJa8KNZRIBh0+KiJS4LRFICJS4FQEIiIFTkUgIlLgVAQiIgVORSAiUuBUBCIiBU5FICJS4P4/LFDRCceuYR0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "s = myokit.Simulation(model)\n",
    "d = s.run(20)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot()\n",
    "ax.set_xlabel('Time')\n",
    "ax.set_ylabel('y')\n",
    "ax.plot(d.time(), d['c.y'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Proposed API\n",
    "\n",
    "At the moment, we can't do this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "ename": "NotImplementedError",
     "evalue": "Sensitivities with respect to parameters used in initial conditions is not implemented (1 / (1 + exp(-c.h * c.K * log(10))) depends on c.h).",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNotImplementedError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-38-4b1b7711a214>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0ms\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmyokit\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mSimulation\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msensitivities\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'c.y'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m'c.h'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'c.K'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/dev/myokit/myokit/_sim/cvodessim.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(***failed resolving arguments***)\u001b[0m\n\u001b[1;32m    176\u001b[0m                     \u001b[0;32mfor\u001b[0m \u001b[0me\u001b[0m \u001b[0;32min\u001b[0m \u001b[0minits\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    177\u001b[0m                         \u001b[0;32mif\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdepends_on\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdeep\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 178\u001b[0;31m                             raise NotImplementedError(\n\u001b[0m\u001b[1;32m    179\u001b[0m                                 \u001b[0;34m'Sensitivities with respect to parameters used'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    180\u001b[0m                                 \u001b[0;34m' in initial conditions is not implemented ('\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNotImplementedError\u001b[0m: Sensitivities with respect to parameters used in initial conditions is not implemented (1 / (1 + exp(-c.h * c.K * log(10))) depends on c.h)."
     ]
    }
   ],
   "source": [
    "s = myokit.Simulation(model, sensitivities=(['c.y'], ['c.h', 'c.K']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reason for this error is that the simulation stores its _state_ and _default state_ as lists of floating point numbers, so no room for fancy expressions!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Requirements\n",
    "\n",
    "- The new API should be backwards compatible wherever possible. Users who do not use the new functionality should not have to change their scripts.\n",
    "- The existing system with a `state`, `default_state` and `reset()`+`pre()` is already slightly confusing to new users. The new system should not be significantly more confusing.\n",
    "- A good rule of thumb might be that if users aren't thinking about initial conditions they won't notice any changes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decisions\n",
    "\n",
    "- The methods `sim.state()` and `sim.set_state()` work with floats\n",
    "- Using `pre()` updates the `default_state` and `default_s_state`, all connection with constants is lost"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Questions \n",
    "\n",
    "- Should links to expressions be maintained? These would all have to be stored python-only. From the point of view of the C code there is no such thing as a dependency of the sensitivity state on the constants - this is implemented purely by particular values in the state-sensitivity matrix\n",
    "- Should set_constant(p) where a model's initial state depends on p affect the default s_state ?\n",
    "- Should calling set_default_state with expressions be allowed. And if so, should this update the default_s_state? (because setting `x0 = p` implies `dx0/dp = 1`)\n",
    "  - **Should some link between the `default_state` and the `default_s_state` be maintained?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Proposal 1A: default_s_state as floats\n",
    "\n",
    "**The simulation turns the initial state into floats, but maintains model constants in symbolic form**\n",
    "\n",
    "- When a simulation is created, the `default_s_state` is set based on the model's initial state, i.e. mostly zeroes, except where dx0/dp is non-zero in which case we evaluate this expression.\n",
    "- Calling `set_state` sets the state (to floats), calling `set_default_state` sets the default state (to floats). Same for `s_state` and `default_s_state`.\n",
    "- Calling `set_constant` affects constants, but this is interpreted as \"during simulations\", so does not affect initial state sensitivities.\n",
    "- Calling `reset()` resets to the default state and `s_state`, which are floats and so not affected by the values of any constants.\n",
    "\n",
    "\n",
    "## Proposal 1B\n",
    "\n",
    "- **Option**: We could make `set_default_state` - when called with expressions - also re-evaluate `default_s_state` to spare the user the effort of evaluating dx0/dp? This would mean you could do `sim.set_constant(p, 3)`, `sim.set_default_state(model.initial_values())`, `sim.reset()` as quite a roundabout way to re-evaluate the initial sensitivities.\n",
    "\n",
    "## Proposal 2: default_s_state as expressions\n",
    "\n",
    "**The simulation maintains an initial_state in symbolic form, and calling reset() re-evaluates it**\n",
    "\n",
    "- When a simulation is created, the `default_s_state` is set based on the model's initial state. The `s_state` is set to an evaluation of that state`.\n",
    "- Calling `set_constant` will affect constants for the rest of the simulation, but will not affect the `s_state`.\n",
    "- Calling `reset` re-evaluates the `default_state` and also the `default_s_state` if this is a list of expressions.\n",
    "  - If the user has manually changed the `default_s_state` to remove links to constants, then the `default_s_state` remains unchanged, even though the `default_state` implies certain connections.\n",
    "- Calling `set_default_state` does not affect the `default_s_state`, because the user may have changed that manually and we don't want to overwrite it?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use-cases: Parameter estimation\n",
    "\n",
    "Below is a series of use cases, all based on parameter estimation, in which we would like to use sensitivities but also re-use the same simulation object with `set_constant` to change parameter values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use case 1: Partially known, partially inferred initial conditions\n",
    "\n",
    "We have a model where some initial conditions are known, and some are derived from parameters we wish to infer.\n",
    "\n",
    "We want to:\n",
    "\n",
    "1. Create a simulation.\n",
    "2. Repeatedly\n",
    "  - Set parameters, affecting the state and state sensitivities.\n",
    "  - Run\n",
    "  \n",
    "There are two possible scenarios: If the initial sensitivities should all be zero except through the initial values depending on parameters, then the model can calculate them.\n",
    "\n",
    "If the initial sensitivities should not be zero, then we need to either (1) do some pre-pacing (losing the link to parameters) or (2) have the user set the initial sensitivities (which could be based on something done with the model, I guess?)\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I think the second case is far more common, and so it's not worth it having the complex machinery needed for the simulation to be able to automatically evaluate a symbolic initial state?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Examples\n",
    "\n",
    "**These are all relatively simple - you can stop reading if you're thinking about set_constant**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. A simple run"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "sim = myokit.Simulation(\n",
    "    model,\n",
    "    protocol,\n",
    "    sensitivities=(['c.y'], ['c.h', 'c.K'])\n",
    ")\n",
    "\n",
    "sim.time() --> 0\n",
    "sim.state() --> A list with 1 float\n",
    "sim.s_state() --> A matrix filled with zeroes\n",
    "\n",
    "sim.initial_time() --> 0\n",
    "sim.initial_state() --> A list with 1 expression\n",
    "sim.initial_s_state() --> A matrix filled with zeroes\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "log, sens = sim.run(20)\n",
    "\n",
    "log --> A data log\n",
    "sens --> A list of matrices (or lists of lists)\n",
    "\n",
    "sim.time() --> 20\n",
    "sim.state() --> The new state\n",
    "sim.s_state() --> A matrix filled with numbers\n",
    "\n",
    "sim.initial_time() --> 0\n",
    "sim.initial_state() --> A list with 1 expression\n",
    "sim.initial_s_state() --> A matrix filled with zeroes\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Manually setting the initial_state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "sim = myokit.Simulation(\n",
    "    model,\n",
    "    protocol,\n",
    "    sensitivites=(['c.y'], ['c.h', 'c.K'])\n",
    ")\n",
    "\n",
    "sim.set_initial_state([c.h + 5*c.K])\n",
    "\n",
    "sim.time() --> 0\n",
    "sim.state() --> A list with 1 float (evaluation of initial state)\n",
    "sim.s_state() --> A matrix filled with zeroes\n",
    "\n",
    "sim.initial_time() --> 0\n",
    "sim.initial_state() --> [c.h + 5*c.K]\n",
    "sim.initial_s_state() --> A matrix filled with zeroes\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "log, sens = sim.run(20)\n",
    "\n",
    "log --> A data log\n",
    "sens --> A list of matrices (or lists of lists)\n",
    "\n",
    "sim.time() --> 20\n",
    "sim.state() --> The new state\n",
    "sim.s_state() --> A matrix filled with numbers\n",
    "\n",
    "sim.initial_time() --> 0\n",
    "sim.initial_state() --> [c.h + 5*c.K]\n",
    "sim.initial_s_state() --> A matrix filled with zeroes\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Manually setting an s_state and initial_s_state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "sim = myokit.Simulation(\n",
    "    model,\n",
    "    protocol,\n",
    "    sensitivites=(['c.y'], ['c.h', 'c.K'])\n",
    ")\n",
    "\n",
    "sim.set_s_state(A matrix filled with numbers)\n",
    "sim.set_initial_s_state(A matrix filled with numbers)\n",
    "\n",
    "sim.time() --> 0\n",
    "sim.state() --> A list with 1 float\n",
    "sim.s_state() --> A matrix filled with numbers\n",
    "\n",
    "sim.initial_time() --> 0\n",
    "sim.initial_state() --> A list with 1 expression\n",
    "sim.initial_s_state() --> A matrix filled with numbers\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "log, sens = sim.run(20)\n",
    "\n",
    "log --> A data log\n",
    "sens --> A list of matrices (or lists of lists)\n",
    "\n",
    "sim.time() --> 20\n",
    "sim.state() --> The new state\n",
    "sim.s_state() --> A matrix filled with numbers\n",
    "\n",
    "sim.initial_time() --> 0\n",
    "sim.initial_state() --> A list with 1 expression\n",
    "sim.initial_s_state() --> A matrix filled with numbers\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Pre-pacing\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "sim = myokit.Simulation(\n",
    "    model,\n",
    "    protocol,\n",
    "    sensitivities=(['c.y'], ['c.h', 'c.K'])\n",
    ")\n",
    "\n",
    "sim.time() --> 0\n",
    "sim.state() --> A list with 1 float\n",
    "sim.s_state() --> A matrix filled with zeroes\n",
    "\n",
    "sim.initial_time() --> 0\n",
    "sim.initial_state() --> A list with 1 expression\n",
    "sim.initial_s_state() --> A matrix filled with zeroes\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "log, sens = sim.pre(20)\n",
    "\n",
    "log --> A data log\n",
    "sens --> A list of matrices (or lists of lists)\n",
    "\n",
    "sim.time() --> 0\n",
    "sim.state() --> The new state\n",
    "sim.s_state() --> A matrix filled with numbers\n",
    "\n",
    "sim.initial_time() --> 0\n",
    "sim.initial_state() --> An updated list with 1 expression\n",
    "sim.initial_s_state() --> A matrix filled with zeroes\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Using set_constant"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maybe in the setting of a fit that uses sensitivities, and where we infer the parameter that the sensitivity depends on?\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}