{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f6aa3188",
   "metadata": {},
   "source": [
    "# Calcium imaging — estimating physiological parameters\n",
    "\n",
    "**Companion to slides 11–12 · Neurofunctional Techniques**\n",
    "\n",
    "In this virtual experiment, we load a cell with different amounts of calcium indicator (different **positive** indicator buffering capacities, $\\kappa_B$), measure the calcium transient following **one action potential**, and use straight-line fits to infer physiological parameters.\n",
    "\n",
    "The model assumes a single well-mixed compartment, rapid linear buffering, and linear calcium extrusion. Simulated measurements include adjustable noise (3% by default). Noise is fixed while sliders move, so differences between conditions reflect parameter changes rather than new random samples.\n",
    "\n",
    "### Slide 11 — Decay time\n",
    "\n",
    "$$\\tau=\\frac{1+\\kappa_S+\\kappa_B}{\\gamma}$$\n",
    "\n",
    "Plotting $\\tau$ against $\\kappa_B$ gives a line with **slope** $1/\\gamma$ and **x-intercept** $-(1+\\kappa_S)$. At $\\kappa_B=0$, the **extrapolated y-intercept** is the unperturbed decay time $\\tau_0$. No experimental measurement is taken at $\\kappa_B=0$.\n",
    "\n",
    "### Slide 12 — Inverse amplitude\n",
    "\n",
    "$$A=\\frac{\\Delta[Ca^{2+}]_T}{1+\\kappa_S+\\kappa_B}\\qquad\\Longrightarrow\\qquad\\frac{1}{A}=\\frac{1+\\kappa_S+\\kappa_B}{\\Delta[Ca^{2+}]_T}$$\n",
    "\n",
    "Plotting $1/A$ against $\\kappa_B$ gives a line with **slope** $1/\\Delta[Ca^{2+}]_T$ and the **same x-intercept** $-(1+\\kappa_S)$. Its **extrapolated** y-intercept gives $1/A_0$.\n",
    "\n",
    "**Try this:** Change one parameter at a time. Predict what will happen to each graph *before* moving the slider. Use **Show fitted values** to check your predictions.\n",
    "\n",
    "**Comparison:** Click **Save current simulation** to retain up to two reference conditions. The dashed curves stay fixed while you move the sliders. The x-axis extends into negative κB values only to show the **extrapolated** x-intercept; negative buffering is not physically accessible. **There is no measured point at $\\kappa_B=0$**: this is the indicator-free limit inferred by extrapolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "815e4d9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from io import BytesIO\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import display, clear_output\n",
    "\n",
    "# All concentrations are in µM; all times are in seconds internally.\n",
    "ks = widgets.IntSlider(value=50, min=0, max=200, step=5, description='κS:', continuous_update=False, style={'description_width':'75px'}, layout=widgets.Layout(width='360px'))\n",
    "gamma = widgets.IntSlider(value=700, min=100, max=1500, step=50, description='γ (s⁻¹):', continuous_update=False, style={'description_width':'75px'}, layout=widgets.Layout(width='360px'))\n",
    "influx = widgets.FloatSlider(value=20, min=1, max=50, step=1, description='Influx (µM):', continuous_update=False, style={'description_width':'90px'}, layout=widgets.Layout(width='360px'))\n",
    "kmin = widgets.IntSlider(value=10, min=1, max=40, step=1, description='Min κB:', continuous_update=False, style={'description_width':'75px'}, layout=widgets.Layout(width='360px'))\n",
    "kmax = widgets.IntSlider(value=100, min=20, max=200, step=10, description='Max κB:', continuous_update=False, style={'description_width':'75px'}, layout=widgets.Layout(width='360px'))\n",
    "npoints = widgets.IntSlider(value=6, min=3, max=12, step=1, description='Data points:', continuous_update=False, style={'description_width':'90px'}, layout=widgets.Layout(width='360px'))\n",
    "noise = widgets.FloatSlider(value=3.0, min=0, max=10, step=0.5, description='Noise (%):', continuous_update=False, style={'description_width':'90px'}, layout=widgets.Layout(width='360px'))\n",
    "show_fit = widgets.Checkbox(value=False, description='Show fitted values and intercepts', indent=False)\n",
    "reset = widgets.Button(description='Reset defaults', icon='refresh')\n",
    "save = widgets.Button(description='Save current simulation', icon='plus', button_style='info')\n",
    "clear_refs = widgets.Button(description='Clear saved simulations', icon='trash')\n",
    "refs = []  # Up to two saved parameter sets\n",
    "status = widgets.HTML(value='<small>No saved simulations (maximum two).</small>')\n",
    "out = widgets.Output()\n",
    "\n",
    "# Render equations as actual images using Matplotlib mathtext.\n",
    "# This is reliable in Google Colab: HTML widgets do not always render LaTeX.\n",
    "def equation_image(formula, width=5.7):\n",
    "    fig = plt.figure(figsize=(width, 0.80), dpi=150)\n",
    "    fig.patch.set_facecolor('white')\n",
    "    fig.text(0.02, 0.5, formula, ha='left', va='center', fontsize=19)\n",
    "    buffer = BytesIO()\n",
    "    fig.savefig(buffer, format='png', dpi=150, bbox_inches='tight', pad_inches=0.12)\n",
    "    plt.close(fig)\n",
    "    return widgets.Image(value=buffer.getvalue(), format='png',\n",
    "                         layout=widgets.Layout(max_width='100%', height='auto'))\n",
    "\n",
    "info = widgets.VBox([\n",
    "    widgets.HTML('<b>Slide 11 · Decay time</b>'),\n",
    "    equation_image(r'$\\tau = \\frac{1}{\\gamma}\\,\\kappa_B + \\frac{1+\\kappa_S}{\\gamma}$'),\n",
    "    widgets.HTML('<small>y = τ &nbsp;&nbsp; x = κ<sub>B</sub></small>'),\n",
    "    widgets.HTML('<b>Slide 12 · Inverse amplitude</b>'),\n",
    "    equation_image(r'$\\frac{1}{A} = \\frac{1}{\\Delta [Ca^{2+}]_T}\\,\\kappa_B + \\frac{1+\\kappa_S}{\\Delta [Ca^{2+}]_T}$'),\n",
    "    widgets.HTML('<small>y = 1/A &nbsp;&nbsp; x = κ<sub>B</sub></small>'),\n",
    "    widgets.HTML('<small>Experimental points: κ<sub>B</sub> &gt; 0. '\n",
    "                 'The values at κ<sub>B</sub> = 0 and the negative x-intercept '\n",
    "                 'are obtained by extrapolation.</small>')\n",
    "], layout=widgets.Layout(width='510px', border='1px solid #d0d7de',\n",
    "                         padding='14px', gap='9px'))\n",
    "# Fixed random deviates ensure that varying model parameters doesn't regenerate noise.\n",
    "noise_rng = np.random.default_rng(20261030)\n",
    "noise_tau = noise_rng.normal(size=12)\n",
    "noise_amp = noise_rng.normal(size=12)\n",
    "\n",
    "def measurements(k_s, g, dca, k_lo, k_hi, count, noise_pct, deviations_tau, deviations_amp):\n",
    "    kb = np.linspace(k_lo, k_hi, count)\n",
    "    tau_ideal = 1000 * (1 + k_s + kb) / g\n",
    "    inv_ideal = (1 + k_s + kb) / dca\n",
    "    # Independent, multiplicative measurement errors; keep positive measurements.\n",
    "    tau_observed = tau_ideal * np.maximum(0.05, 1 + noise_pct / 100 * deviations_tau[:count])\n",
    "    inv_observed = inv_ideal * np.maximum(0.05, 1 + noise_pct / 100 * deviations_amp[:count])\n",
    "    return kb, tau_observed, inv_observed\n",
    "\n",
    "def draw(*_):\n",
    "    k_s, g, dca = ks.value, gamma.value, influx.value\n",
    "    kb, tau_ms, inv_amp = measurements(k_s, g, dca, kmin.value, kmax.value,\n",
    "                                         npoints.value, noise.value, noise_tau, noise_amp)\n",
    "    fit_tau = np.polyfit(kb, tau_ms, 1)\n",
    "    fit_inv = np.polyfit(kb, inv_amp, 1)\n",
    "    xint_tau = -fit_tau[1] / fit_tau[0]\n",
    "    xint_inv = -fit_inv[1] / fit_inv[0]\n",
    "\n",
    "    # Include the negative x-intercepts of both current and saved conditions.\n",
    "    all_ks = [k_s] + [r['ks'] for r in refs]\n",
    "    left_edge = -max(1.18 * (1 + max(all_ks)), 0.28 * kmax.value)\n",
    "    right_edge = kmax.value * 1.07\n",
    "    xx = np.linspace(left_edge, right_edge, 400)\n",
    "    fig, axes = plt.subplots(1, 2, figsize=(13.2, 5.0), constrained_layout=True)\n",
    "    colors = ['#9b59b6', '#d97706']\n",
    "\n",
    "    for i, r in enumerate(refs):\n",
    "        c = colors[i]\n",
    "        label = f\"Saved {i+1}: κS={r['ks']}, γ={r['gamma']}, influx={r['influx']} µM\"\n",
    "        # Saved curves extend to their x-intercepts, independently of the current sliders.\n",
    "        old_kb_fit, old_tau_fit, old_inv_fit = measurements(r['ks'], r['gamma'], r['influx'], r['kmin'], r['kmax'], r['npoints'], r['noise'], r['noise_tau'], r['noise_amp'])\n",
    "        y1 = np.polyval(np.polyfit(old_kb_fit, old_tau_fit, 1), xx)\n",
    "        y2 = np.polyval(np.polyfit(old_kb_fit, old_inv_fit, 1), xx)\n",
    "        axes[0].plot(xx, y1, '--', color=c, lw=1.8, alpha=.85, label=label)\n",
    "        axes[1].plot(xx, y2, '--', color=c, lw=1.8, alpha=.85, label=label)\n",
    "        old_kb, old_tau, old_inv = measurements(\n",
    "            r['ks'], r['gamma'], r['influx'], r['kmin'], r['kmax'], r['npoints'],\n",
    "            r['noise'], r['noise_tau'], r['noise_amp'])\n",
    "        axes[0].scatter(old_kb, old_tau, color=c, s=25, alpha=.8)\n",
    "        axes[1].scatter(old_kb, old_inv, color=c, s=25, alpha=.8)\n",
    "\n",
    "    axes[0].scatter(kb, tau_ms, color='#087f8c', s=50, zorder=4, label='Current measurements')\n",
    "    axes[1].scatter(kb, inv_amp, color='#087f8c', s=50, zorder=4, label='Current measurements')\n",
    "    if show_fit.value:\n",
    "        axes[0].plot(xx, np.polyval(fit_tau, xx), '-', color='#087f8c', lw=2, label='Current linear fit')\n",
    "        axes[1].plot(xx, np.polyval(fit_inv, xx), '-', color='#087f8c', lw=2, label='Current linear fit')\n",
    "        for ax, xint, yint in zip(axes, [xint_tau,xint_inv], [fit_tau[1],fit_inv[1]]):\n",
    "            ax.scatter([xint,0], [0,yint], color='#c2410c', s=62, zorder=6)\n",
    "            ax.annotate(f'x = {xint:.1f}', (xint,0), xytext=(0,-20), textcoords='offset points', ha='center', fontsize=9)\n",
    "\n",
    "    for ax in axes:\n",
    "        ax.set_xlabel('Indicator buffering capacity κB')\n",
    "        ax.set_xlim(left_edge, right_edge)\n",
    "        # Display negative x values to show extrapolated x-intercepts; y remains physical (≥ 0).\n",
    "        ax.set_ylim(bottom=0)\n",
    "        ax.axvline(0, color='#111827', linewidth=1.5, zorder=2)\n",
    "        ax.axhline(0, color='#111827', linewidth=1.5, zorder=2)\n",
    "        ax.grid(alpha=0.17)\n",
    "        ax.spines[['top','right']].set_visible(False)\n",
    "        ax.legend(loc='upper left', fontsize=7.5)\n",
    "    axes[0].set_title('Slide 11 · Decay time')\n",
    "    axes[0].set_ylabel('Decay time τ (ms)')\n",
    "    axes[1].set_title('Slide 12 · Inverse amplitude')\n",
    "    axes[1].set_ylabel('1 / A (µM⁻¹)')\n",
    "    with out:\n",
    "        clear_output(wait=True)\n",
    "        display(fig)\n",
    "        plt.close(fig)\n",
    "        if show_fit.value:\n",
    "            print('PARAMETERS RECOVERED FROM THE CURRENT STRAIGHT-LINE FITS')\n",
    "            print(f'Slide 11  slope = {fit_tau[0]:.4f} ms per κB   →  γ = {1000/fit_tau[0]:.1f} s⁻¹')\n",
    "            print(f'          x-intercept = {xint_tau:.2f}   →  κS = {-xint_tau-1:.1f}')\n",
    "            print(f'          y-intercept = {fit_tau[1]:.2f} ms   →  τ₀ = {fit_tau[1]:.2f} ms')\n",
    "            print(f'Slide 12  slope = {fit_inv[0]:.4f} µM⁻¹ per κB   →  total influx = {1/fit_inv[0]:.2f} µM')\n",
    "            print(f'          x-intercept = {xint_inv:.2f}   →  κS = {-xint_inv-1:.1f}')\n",
    "            print(f'          y-intercept = {fit_inv[1]:.4f} µM⁻¹   →  A₀ = {1/fit_inv[1]:.4f} µM ({1000/fit_inv[1]:.1f} nM)')\n",
    "        else:\n",
    "            print('Noisy simulated measurements: all κB > 0.  Select “Show fitted values and intercepts” to extrapolate to κB = 0 and the negative x-intercept.')\n",
    "        if refs:\n",
    "            print('Saved simulations: ' + ' | '.join(f\"{i+1}: κS={r['ks']}, γ={r['gamma']}, influx={r['influx']} µM\" for i,r in enumerate(refs)))\n",
    "\n",
    "for control in [ks, gamma, influx, kmin, kmax, npoints, noise, show_fit]:\n",
    "    control.observe(draw, names='value')\n",
    "def on_save(_):\n",
    "    record = dict(ks=ks.value, gamma=gamma.value, influx=influx.value,\n",
    "                  kmin=kmin.value, kmax=kmax.value, npoints=npoints.value,\n",
    "                  noise=noise.value, noise_tau=noise_tau.copy(), noise_amp=noise_amp.copy())\n",
    "    if len(refs) == 2:\n",
    "        refs.pop(0)  # Keep the two most recently saved simulations.\n",
    "    refs.append(record)\n",
    "    status.value = f'<small>{len(refs)}/2 simulations saved. A third save replaces the oldest.</small>'\n",
    "    draw()\n",
    "\n",
    "def on_clear(_):\n",
    "    refs.clear()\n",
    "    status.value = '<small>No saved simulations (maximum two).</small>'\n",
    "    draw()\n",
    "\n",
    "save.on_click(on_save)\n",
    "clear_refs.on_click(on_clear)\n",
    "\n",
    "def on_reset(_):\n",
    "    refs.clear()\n",
    "    status.value = '<small>No saved simulations (maximum two).</small>'\n",
    "    for widget, value in [(ks,50),(gamma,700),(influx,20),(kmin,10),(kmax,100),(npoints,6),(noise,3.0),(show_fit,False)]:\n",
    "        widget.value=value\n",
    "reset.on_click(on_reset)\n",
    "left = widgets.VBox([widgets.HTML('<b>Physiological parameters</b>'), ks, gamma, influx,\n",
    "                     widgets.HTML('<b>Virtual experiment (κB > 0 only)</b>'), kmin, kmax, npoints, noise, show_fit, widgets.HBox([save, clear_refs]), status, reset])\n",
    "display(widgets.HBox([left, info], layout=widgets.Layout(gap='22px', flex_flow='row wrap')))\n",
    "display(out)\n",
    "draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03d1bc58",
   "metadata": {},
   "source": [
    "### Interpretation\n",
    "\n",
    "- Increasing **$\\kappa_S$** moves both lines upward and shifts their x-intercepts leftward, without changing their slopes.\n",
    "- Increasing **$\\gamma$** reduces the slope and intercept of the decay-time plot, but does not affect the inverse-amplitude plot.\n",
    "- Increasing **total calcium influx** reduces the slope and intercept of the inverse-amplitude plot, but does not affect the decay-time plot.\n",
    "- The **unperturbed** values $\\tau_0$ and $A_0$ refer to the **extrapolated** case of **no added indicator** ($\\kappa_B=0$), *not* to the absence of endogenous buffering. No data point is measured at $\\kappa_B=0$.\n",
    "\n",
    "**Model limitation:** In real experiments, added indicators may have nonlinear or slow binding, fluorescence calibration is imperfect, and measurements contain noise. The **Noise (%)** slider adds reproducible scatter to simulated measurements; the fitted line is estimated from those noisy points and therefore need not recover the exact parameters. Set noise to 0% to recover the ideal case.\n",
    "\n",
    "**Graph note:** The horizontal and vertical axes intersect at zero. Negative κB is shown solely to visualize the extrapolated x-intercept. Save up to two simulated conditions to compare them with the current condition.\n",
    "\n",
    "**How to use:** Start with 0% noise to show the ideal relationships, then increase to 3–5% to illustrate experimental variability. Save up to two conditions and change one parameter at a time.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}