{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "aa1351e1a3215a2",
   "metadata": {},
   "source": [
    "# Examples for symbolic calculations with CM\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9988a172",
   "metadata": {},
   "source": [
    "A `SymbolicCompartmentalModel` can operate in two modes: **numeric** (returning float values) and **symbolic** (returning closed-form analytical expressions via SymPy). In symbolic mode, every model property is expressed as a formula in terms of the named parameters, which is useful for gaining analytical insight, deriving limiting cases, and understanding model behaviour across the full parameter space without running numerical simulations."
   ]
  },
  {
   "metadata": {},
   "cell_type": "code",
   "outputs": [],
   "execution_count": null,
   "source": [
    "from symbolic_compartmental_model import SymbolicCompartmentalModel\n",
    "import sympy\n",
    "import matplotlib.pyplot as plt\n",
    "from pathlib import Path"
   ],
   "id": "9b1d81dae13a3110"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "outputs": [],
   "execution_count": null,
   "source": [
    "cm3 = SymbolicCompartmentalModel(n_states=2)\n",
    "k1 = cm3.add_parameter(symbol=\"k1\", lb=0.0, ub=10.0)\n",
    "k2 = cm3.add_parameter(symbol=\"k2\", lb=0.0, ub=10.0)\n",
    "cm3.contributed_turnovers = [[-k1, 0.0], [0.0, -k2]]\n",
    "cm3.observed_pool_weights = [0.5, 0.5]"
   ],
   "id": "64055b1e792b0d1a"
  },
  {
   "cell_type": "markdown",
   "id": "19ed3b90",
   "metadata": {},
   "source": [
    "We define a two-state CM where each state has an independent turnover rate (`k1` and `k2`), and both contribute equally to the observed signal. Because the off-diagonal entries are zero the two states do not exchange molecules — the labeling curve is a symmetric bi-exponential `0.5·exp(−k1·t) + 0.5·exp(−k2·t)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69f6a670bb04625f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-04T15:10:34.716178Z",
     "start_time": "2025-09-04T15:10:34.007317Z"
    }
   },
   "outputs": [],
   "source": [
    "print(\"\\n\\nmean age:\")\n",
    "display(sympy.simplify(cm3.as_symbolic(\"mean_age\")))\n",
    "print(\"\\n\\nmean residence time:\")\n",
    "display(sympy.simplify(cm3.as_symbolic(\"mean_residence_time\")))\n",
    "print(\"\\n\\nexpected decay rate:\")\n",
    "display(sympy.simplify(cm3.as_symbolic(\"expected_decay_rate\")))\n",
    "print(\"\\n\\nage CDF:\")\n",
    "display(sympy.simplify(cm3.as_symbolic(\"age_cdf\")))\n",
    "print(\"\\n\\nage PDF:\")\n",
    "display(sympy.simplify(cm3.as_symbolic(\"age_pdf\")))\n",
    "print(\"\\n\\nresidence time CDF:\")\n",
    "display(sympy.simplify(cm3.as_symbolic(\"residence_time_cdf\")))\n",
    "print(\"\\n\\nresidence time PDF:\")\n",
    "display(sympy.simplify(cm3.as_symbolic(\"residence_time_pdf\")))\n",
    "print(\"\\n\\ndecay rate:\")\n",
    "display(sympy.simplify(cm3.as_symbolic(\"decay_rate\")))"
   ]
  },
  {
   "metadata": {},
   "cell_type": "code",
   "outputs": [],
   "execution_count": null,
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(5, 3), dpi=150)\n",
    "cm3.plot(\"decay_rate\", x=[1.0, 3.0], ax=ax)\n",
    "fig.tight_layout()\n",
    "if Path(\"../results\").exists():\n",
    "    fig.savefig(\"../results/decay_rate.svg\")"
   ],
   "id": "abe382459157c13e"
  },
  {
   "cell_type": "markdown",
   "id": "a878dbad",
   "metadata": {},
   "source": [
    "`as_symbolic(key)` returns a SymPy expression for any of the following model properties:\n",
    "\n",
    "| Key | Meaning |\n",
    "|-----|---------|\n",
    "| `mean_age` | Expected age of a randomly chosen molecule (area under `f(t)`) |\n",
    "| `mean_residence_time` | Expected *remaining* lifetime from now (forward-looking) |\n",
    "| `expected_decay_rate` | Instantaneous turnover rate averaged over the age distribution |\n",
    "| `age_cdf` / `age_pdf` | Cumulative / probability density of molecular ages |\n",
    "| `residence_time_cdf` / `residence_time_pdf` | Distribution of remaining lifetimes |\n",
    "| `decay_rate` | Hazard function `h(t)` — instantaneous clearance rate given survival to age `t` |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1218db24606327a8",
   "metadata": {},
   "source": [
    "## Growing system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65c5f3136cdc4dd6",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-06T08:37:14.150611Z",
     "start_time": "2025-09-06T08:37:14.042320Z"
    }
   },
   "outputs": [],
   "source": [
    "cm3.growth_rate = 0.0\n",
    "print(\"\\n\\nmean age:\")\n",
    "display(sympy.simplify(cm3.as_symbolic(\"mean_age\")))\n",
    "print(\"\\n\\nmean residence time:\")\n",
    "display(sympy.simplify(cm3.as_symbolic(\"mean_residence_time\")))\n",
    "cm3.growth_rate = 1.5\n",
    "print(\"\\n\\nmean age:\")\n",
    "display(sympy.simplify(cm3.as_symbolic(\"mean_age\")))\n",
    "print(\"\\n\\nmean residence time:\")\n",
    "display(sympy.simplify(cm3.as_symbolic(\"mean_residence_time\")))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9e71bbb",
   "metadata": {},
   "source": [
    "The `decay_rate` (hazard function) `h(t)` describes how the probability of clearance changes as a molecule ages. For a single exponential `h(t)` is constant; for a multi-state model it varies with time, reflecting the shifting composition of the surviving population. At `t = 0` the fast-clearing state dominates, driving a high initial hazard; as fast-clearing molecules are depleted the hazard decreases toward the rate of the slower state."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b624b963",
   "metadata": {},
   "source": [
    "In a growing cell population, cell division continuously dilutes all molecules, adding an apparent clearance term `μ` to every state. Setting `growth_rate = μ` incorporates this into the model.\n",
    "\n",
    "This affects **`mean_residence_time`** (forward-looking — accounts for future dilution by growth) but **not `mean_age`** (retrospective — describes molecules that already exist). Setting `growth_rate = 0` recovers the non-growing case."
   ]
  }
 ],
 "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": "ipython3",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}