{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interactive single-compartment model of Ca²⁺ dynamics\n",
    "\n",
    "**Neurofunctional Techniques · Lesson 8**\n",
    "\n",
    "Explore how calcium influx, endogenous and indicator buffering, and extrusion shape free intracellular Ca²⁺ signals. The model assumes a well-mixed compartment, fast non-saturating buffers and linear extrusion. It shows **free Ca²⁺ elevation above baseline**, not fluorescence.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explore the model\n",
    "\n",
    "Run the cell below (▶), adjust the sliders, and compare conditions using **Save reference trace**. The equations are displayed in mathematical notation beside the controls.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1edfd7a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import display, clear_output\n",
    "\n",
    "style = {'description_width': 'initial'}\n",
    "mode = widgets.ToggleButtons(options=['Single AP', 'AP train'], value='Single AP', description='Mode:')\n",
    "ks = widgets.IntSlider(value=50, min=0, max=200, step=5, description='Endogenous κS', style=style, continuous_update=False)\n",
    "kb = widgets.IntSlider(value=20, min=0, max=200, step=5, description='Indicator κB', style=style, continuous_update=False)\n",
    "influx = widgets.IntSlider(value=20, min=1, max=50, step=1, description='Total Ca²⁺ / AP (µM)', style=style, continuous_update=False)\n",
    "gamma = widgets.IntSlider(value=700, min=100, max=1500, step=50, description='Extrusion γ (s⁻¹)', style=style, continuous_update=False)\n",
    "freq = widgets.IntSlider(value=20, min=1, max=100, step=1, description='AP frequency (Hz)', style=style, continuous_update=False)\n",
    "count = widgets.IntSlider(value=50, min=2, max=150, step=1, description='Number of APs', style=style, continuous_update=False)\n",
    "save = widgets.Button(description='Save reference trace', button_style='info', icon='save')\n",
    "clear = widgets.Button(description='Clear reference', icon='trash')\n",
    "reset = widgets.Button(description='Reset defaults', icon='refresh')\n",
    "output = widgets.Output()\n",
    "reference = {'data': None}\n",
    "\n",
    "def make_trace(settings, t):\n",
    "    m, s, b, q, g, f, number = settings\n",
    "    B = 1+s+b\n",
    "    A = q*1000/B  # µM -> nM\n",
    "    tau = B/g\n",
    "    times = np.array([0.]) if m == 'Single AP' else np.arange(number)/f\n",
    "    y = np.zeros_like(t)\n",
    "    for ti in times:\n",
    "        mask = t >= ti\n",
    "        y[mask] += A*np.exp(-(t[mask]-ti)/tau)\n",
    "    return y, times, A, tau\n",
    "\n",
    "def current():\n",
    "    return (mode.value, ks.value, kb.value, influx.value, gamma.value, freq.value, count.value)\n",
    "\n",
    "def draw(*_):\n",
    "    settings = current()\n",
    "    m,s,b,q,g,f,number = settings\n",
    "    tau = (1+s+b)/g\n",
    "    last = 0 if m == 'Single AP' else (number-1)/f\n",
    "    end = max(0.65, last+5*tau)\n",
    "    if reference['data'] is not None:\n",
    "        r = reference['data']\n",
    "        r_tau = (1+r[1]+r[2])/r[4]\n",
    "        r_last = 0 if r[0]=='Single AP' else (r[6]-1)/r[5]\n",
    "        if r[0] == m:\n",
    "            end = max(end, r_last+5*r_tau)\n",
    "    t = np.linspace(-0.03, end, 2400)\n",
    "    y, times, A, tau = make_trace(settings,t)\n",
    "    with output:\n",
    "        clear_output(wait=True)\n",
    "        fig, ax = plt.subplots(figsize=(9.5,4.5))\n",
    "        if reference['data'] is not None and reference['data'][0] == m:\n",
    "            yr,_,_,_ = make_trace(reference['data'],t)\n",
    "            ax.plot(t,yr,color='0.55',lw=2,ls='--',label='Saved reference')\n",
    "        ax.plot(t,y,color='#087f8c',lw=2.5,label='Current condition')\n",
    "        if m == 'AP train':\n",
    "            plateau = q*1000*f/g\n",
    "            ax.axhline(plateau,color='#bb6d1c',lw=1.5,ls=':',label='Long-train mean plateau')\n",
    "            for ti in times:\n",
    "                ax.axvline(ti,color='0.7',lw=.5,alpha=.28)\n",
    "        ax.set(xlabel='Time (s)',ylabel='Δ free [Ca²⁺] (nM)',title=m)\n",
    "        ax.set_xlim(-0.03,end)\n",
    "        ax.set_ylim(bottom=0)\n",
    "        ax.grid(alpha=.18)\n",
    "        ax.legend(loc='upper right')\n",
    "        plt.tight_layout()\n",
    "        plt.show()\n",
    "        print(f'A = {A:.1f} nM  |  τ = {tau*1000:.1f} ms  |  A × τ = {A*tau:.2f} nM·s')\n",
    "        if m == 'AP train':\n",
    "            print(f'Long-train time-averaged plateau P = {plateau:.1f} nM; finite-train maximum = {y.max():.1f} nM')\n",
    "\n",
    "def save_click(_):\n",
    "    reference['data']=current()\n",
    "    draw()\n",
    "def clear_click(_):\n",
    "    reference['data']=None\n",
    "    draw()\n",
    "def reset_click(_):\n",
    "    reference['data']=None\n",
    "    mode.value='Single AP';ks.value=50;kb.value=20;influx.value=20\n",
    "    gamma.value=700;freq.value=20;count.value=50\n",
    "    draw()\n",
    "save.on_click(save_click);clear.on_click(clear_click);reset.on_click(reset_click)\n",
    "for control in (mode,ks,kb,influx,gamma,freq,count):\n",
    "    control.observe(draw,names='value')\n",
    "\n",
    "def toggle_train(*_):\n",
    "    show = mode.value=='AP train'\n",
    "    freq.layout.display='' if show else 'none'\n",
    "    count.layout.display='' if show else 'none'\n",
    "mode.observe(toggle_train,names='value')\n",
    "toggle_train()\n",
    "# Render equations as an image: reliable mathematical typesetting in Colab.\n",
    "# This avoids displaying raw LaTeX commands in an HTML widget.\n",
    "import io\n",
    "fig_eq, ax_eq = plt.subplots(figsize=(5.3, 3.7), dpi=150)\n",
    "ax_eq.axis('off')\n",
    "ax_eq.text(0.03, 0.94, 'Single-compartment model', fontsize=15,\n",
    "           fontweight='bold', va='top', transform=ax_eq.transAxes)\n",
    "for y, formula in [(0.76, r'$A = \\frac{\\Delta [Ca^{2+}]_T}{1+\\kappa_S+\\kappa_B}$'),\n",
    "                   (0.53, r'$\\tau = \\frac{1+\\kappa_S+\\kappa_B}{\\gamma}$'),\n",
    "                   (0.30, r'$P = f A \\tau = \\frac{f\\,\\Delta [Ca^{2+}]_T}{\\gamma}$')]:\n",
    "    ax_eq.text(0.04, y, formula, fontsize=19, va='center', transform=ax_eq.transAxes)\n",
    "ax_eq.text(0.03, 0.09, 'A: free Ca²⁺ rise per AP   |   τ: decay time constant',\n",
    "           fontsize=10, transform=ax_eq.transAxes)\n",
    "ax_eq.text(0.03, 0.025, 'P: long-train mean plateau',\n",
    "           fontsize=10, transform=ax_eq.transAxes)\n",
    "buf = io.BytesIO()\n",
    "fig_eq.savefig(buf, format='png', bbox_inches='tight', facecolor='white')\n",
    "plt.close(fig_eq)\n",
    "equations = widgets.Image(value=buf.getvalue(), format='png',\n",
    "                          layout=widgets.Layout(width='490px'))\n",
    "controls = widgets.VBox([mode,ks,kb,influx,gamma,freq,count,\n",
    "                         widgets.HBox([save,clear,reset])],\n",
    "                        layout=widgets.Layout(width='560px'))\n",
    "display(widgets.VBox([widgets.HBox([controls,equations],\n",
    "                                    layout=widgets.Layout(align_items='flex-start',flex_flow='row wrap')),\n",
    "                      output]))\n",
    "draw()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Model scope:** The predicted plateau is the time-averaged steady-state elevation for a long regular AP train. Finite trains may not reach it. The model does not include spatial Ca²⁺ gradients, buffer saturation, or indicator fluorescence kinetics.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}