diff --git a/patch_sim/__init__.py b/patch_sim/__init__.py index dcecaf86..fc54d116 100644 --- a/patch_sim/__init__.py +++ b/patch_sim/__init__.py @@ -23,6 +23,8 @@ GVPoint, HyperpolarizationAnalysisResult, ImpedanceProfile, + InactivationAnalysisResult, + InactivationPoint, IVAnalysisResult, IVPoint, PassiveProperties, @@ -48,6 +50,7 @@ boltzmann, compute_fi_point, compute_gv, + compute_inactivation, compute_iv_point, compute_sag_point, compute_sfa, @@ -78,6 +81,7 @@ CURRENT_PROTOCOLS, DEFAULT_Q10, DEFAULT_T_REF, + INACTIVATION_PROTOCOL, VOLTAGE_CLAMP, VOLTAGE_PROTOCOLS, ) @@ -124,6 +128,8 @@ "GVPoint", "HyperpolarizationAnalysisResult", "ImpedanceProfile", + "InactivationAnalysisResult", + "InactivationPoint", "IVAnalysisResult", "IVPoint", "PassiveProperties", @@ -149,6 +155,7 @@ "boltzmann", "compute_fi_point", "compute_gv", + "compute_inactivation", "compute_iv_point", "compute_sag_point", "compute_sfa", @@ -168,6 +175,7 @@ "CURRENT_PROTOCOLS", "DEFAULT_Q10", "DEFAULT_T_REF", + "INACTIVATION_PROTOCOL", "MEMBRANE_TEST_CURRENT", "MEMBRANE_TEST_POST_MS", "MEMBRANE_TEST_PRE_MS", diff --git a/patch_sim/analysis/__init__.py b/patch_sim/analysis/__init__.py index fa555502..cd4614d3 100644 --- a/patch_sim/analysis/__init__.py +++ b/patch_sim/analysis/__init__.py @@ -10,6 +10,7 @@ fi_curve: F-I curve construction from current clamp multi-sweep results. hyperpolarization: Sag and rebound analysis from hyperpolarization sweeps. impedance: Impedance-profile analysis from a chirp current-clamp run. + inactivation: Steady-state inactivation (h∞) curve from a two-pulse VC run. iv_curve: I-V curve construction from voltage clamp multi-sweep results. membrane_test: Dedicated membrane test for passive property characterisation. passive_properties: Passive membrane property extraction (R_in, τₘ, Cₘ). @@ -61,6 +62,11 @@ analyze_impedance, impedance_unavailable_reason, ) +from .inactivation import ( + InactivationAnalysisResult, + InactivationPoint, + compute_inactivation, +) from .iv_curve import IVAnalysisResult, IVPoint, analyze_iv, compute_iv_point from .membrane_test import ( MEMBRANE_TEST_CURRENT, @@ -108,6 +114,7 @@ "compute_dvdt", "compute_fi_point", "compute_gv", + "compute_inactivation", "compute_iv_point", "compute_sag_point", "compute_sfa", @@ -134,6 +141,8 @@ "FIPoint", "HyperpolarizationAnalysisResult", "ImpedanceProfile", + "InactivationAnalysisResult", + "InactivationPoint", "SagPoint", "GVAnalysisResult", "GVPoint", diff --git a/patch_sim/analysis/inactivation.py b/patch_sim/analysis/inactivation.py new file mode 100644 index 00000000..7858dd4a --- /dev/null +++ b/patch_sim/analysis/inactivation.py @@ -0,0 +1,211 @@ +"""Steady-state inactivation (h∞) curve analysis for two-pulse voltage clamp. + +Computes normalized channel availability h∞(V) = I_peak(V) / I_peak_max from an +existing :class:`~patch_sim.analysis.iv_curve.IVAnalysisResult` (whose stimulus +window is the fixed test pulse of a two-pulse protocol, so each +:class:`~patch_sim.analysis.iv_curve.IVPoint` carries the test-pulse peak inward +current at one conditioning prepulse voltage) and fits a *decreasing* +two-parameter Boltzmann sigmoid to characterise voltage-dependent inactivation. +The Boltzmann math (:func:`~patch_sim.analysis.gv_curve.boltzmann`) and the +fitted-parameter container (:class:`~patch_sim.analysis.gv_curve.BoltzmannFit`) +are shared with the activation g-V curve analysis. + +Data classes: + InactivationPoint: Per-prepulse availability record. + InactivationAnalysisResult: Aggregated h∞ analysis output. +""" + +from __future__ import annotations + +import dataclasses +import logging + +import numpy as np +from scipy.optimize import curve_fit + +from .gv_curve import BoltzmannFit, boltzmann +from .iv_curve import IVAnalysisResult + +logger = logging.getLogger(__name__) + +# Default initial guesses and bounds for the decreasing-Boltzmann curve_fit. +# The half-inactivation voltage of fast Na+ channels is more hyperpolarized than +# their half-activation voltage, hence the more negative guess and wider range. +_VHALF_GUESS = -60.0 # mV — typical Na+ half-inactivation +_K_GUESS = 7.0 # mV — typical slope factor +_VHALF_BOUNDS = (-120.0, 20.0) +_K_BOUNDS = (0.5, 40.0) + + +@dataclasses.dataclass +class InactivationPoint: + """Availability measurement for a single conditioning prepulse voltage. + + Attributes: + prepulse_voltage: Conditioning prepulse voltage held before the fixed + test pulse (mV). + peak_inward_current: Most negative (inward) current measured during the + test pulse for this prepulse (µA/cm²). + h_normalized: Availability normalized to h∞ = I_peak / I_peak_max, where + I_peak_max is the most negative peak across all prepulses + (dimensionless, clamped to 0–1). + """ + + prepulse_voltage: float + peak_inward_current: float + h_normalized: float + + +@dataclasses.dataclass +class InactivationAnalysisResult: + """Complete steady-state inactivation analysis from a two-pulse simulation. + + Points are sorted in ascending order of prepulse voltage. The convenience + properties ``prepulse_voltages``, ``peak_inward_currents``, and + ``h_normalized_values`` extract the corresponding field from each point on + demand. + + The fitted Boltzmann parameterizes a *decreasing* sigmoid + ``h∞(V) = 1 / (1 + exp((V - v_half) / k))`` with ``k > 0`` — equivalently + ``boltzmann(V, v_half, -k)`` using the shared + :func:`~patch_sim.analysis.gv_curve.boltzmann` function (which is the + *increasing* activation form). Callers reconstructing the fit curve must + therefore negate ``k``. + + Attributes: + points: Per-prepulse availability records, sorted by ascending prepulse + voltage. + boltzmann: Fitted Boltzmann parameters (may have ``converged=False``). + ``v_half`` is the half-inactivation voltage (mV) and ``k`` the + slope factor (mV, positive); the curve they describe is decreasing. + """ + + points: list[InactivationPoint] + boltzmann: BoltzmannFit + + @property + def prepulse_voltages(self) -> list[float]: + """Conditioning prepulse voltages in mV, sorted ascending.""" + return [p.prepulse_voltage for p in self.points] + + @property + def peak_inward_currents(self) -> list[float]: + """Test-pulse peak inward currents in µA/cm² at each prepulse.""" + return [p.peak_inward_current for p in self.points] + + @property + def h_normalized_values(self) -> list[float]: + """h∞ availability values (dimensionless, 0–1) at each prepulse.""" + return [p.h_normalized for p in self.points] + + +def _decreasing_boltzmann( + V: float | np.ndarray, + v_half: float, + k: float, +) -> float | np.ndarray: + """Evaluate the decreasing (inactivation) two-parameter Boltzmann sigmoid. + + Computes ``1 / (1 + exp((V - v_half) / k))`` by negating *k* and delegating + to the shared (increasing) :func:`~patch_sim.analysis.gv_curve.boltzmann` + function, so the activation and inactivation analyses share one + implementation. Used as the model passed to ``scipy.optimize.curve_fit``. + + Args: + V: Membrane voltage in mV. May be a scalar or a 1-D NumPy array. + v_half: Half-inactivation voltage (mV). + k: Slope factor (mV); positive for a decreasing sigmoid. + + Returns: + Sigmoid value in [0, 1]. Returns a NumPy array when *V* is an array, + or a float when *V* is a scalar. + """ + return boltzmann(V, v_half, -k) + + +def compute_inactivation(iv_result: IVAnalysisResult) -> InactivationAnalysisResult: + """Compute normalized availability vs. prepulse voltage from I-V results. + + Each :class:`~patch_sim.analysis.iv_curve.IVPoint` in *iv_result* is assumed + to come from one sweep of a two-pulse inactivation protocol: its + ``voltage_step`` is the conditioning prepulse voltage and its + ``peak_inward_current`` is the most negative current measured during the + fixed test pulse. Availability is:: + + h∞(V) = I_peak(V) / I_peak_max + + where ``I_peak_max`` is the most negative peak across all prepulses; the + ratio is clamped to [0, 1] (a non-negative test-pulse peak — full + inactivation, possibly net outward — maps to 0.0). + + A decreasing Boltzmann sigmoid ``h∞(V) = 1 / (1 + exp((V - v_half) / k))`` + is then fitted via ``scipy.optimize.curve_fit``. When the fit does not + converge, :attr:`BoltzmannFit.converged` is ``False`` and the parameter + values are ``0.0`` and ``1.0`` respectively. The returned ``k`` is positive + and the fit curve is reconstructed as ``boltzmann(V, v_half, -k)``. + + Args: + iv_result: Pre-computed I-V analysis result whose stimulus window is the + two-pulse protocol's fixed test pulse, with one point per + conditioning prepulse voltage. + + Returns: + An :class:`InactivationAnalysisResult` with sorted per-prepulse + availability records and Boltzmann fit parameters. ``points`` is empty + only when *iv_result* itself is empty. + """ + # Sentinel returned on failure; parameter values are never exposed to the + # user because all callers check BoltzmannFit.converged before using them. + _null_fit = BoltzmannFit(v_half=0.0, k=1.0, converged=False) + + if not iv_result.points: + return InactivationAnalysisResult(points=[], boltzmann=_null_fit) + + # iv_result.points is already sorted by voltage_step (analyze_iv sorts), but + # sort defensively so the output ordering does not depend on the caller. + ordered = sorted(iv_result.points, key=lambda p: p.voltage_step) + voltages = [p.voltage_step for p in ordered] + peaks = [p.peak_inward_current for p in ordered] + + i_max = min(peaks) # most negative ⇒ largest available inward current + if i_max >= 0.0: + # No inward current at any prepulse — pathological (e.g. no Na+ channel + # or test pulse below activation threshold). Report zero availability. + points = [ + InactivationPoint( + prepulse_voltage=v, peak_inward_current=p, h_normalized=0.0 + ) + for v, p in zip(voltages, peaks) + ] + return InactivationAnalysisResult(points=points, boltzmann=_null_fit) + + h_norm = [min(1.0, max(0.0, p / i_max)) for p in peaks] + points = [ + InactivationPoint(prepulse_voltage=v, peak_inward_current=p, h_normalized=hn) + for v, p, hn in zip(voltages, peaks, h_norm) + ] + + if len(points) < 2: + return InactivationAnalysisResult(points=points, boltzmann=_null_fit) + + v_arr = np.array(voltages, dtype=float) + h_arr = np.array(h_norm, dtype=float) + + fit = _null_fit + try: + popt, _ = curve_fit( + _decreasing_boltzmann, + v_arr, + h_arr, + p0=[_VHALF_GUESS, _K_GUESS], + bounds=( + [_VHALF_BOUNDS[0], _K_BOUNDS[0]], + [_VHALF_BOUNDS[1], _K_BOUNDS[1]], + ), + maxfev=2000, + ) + fit = BoltzmannFit(v_half=float(popt[0]), k=float(popt[1]), converged=True) + except (RuntimeError, ValueError) as exc: + logger.debug("Inactivation Boltzmann fit did not converge: %s", exc) + + return InactivationAnalysisResult(points=points, boltzmann=fit) diff --git a/patch_sim/constants.py b/patch_sim/constants.py index 28f2fd1b..0a65da00 100644 --- a/patch_sim/constants.py +++ b/patch_sim/constants.py @@ -53,6 +53,15 @@ #: the UI, so it gets its own constant rather than a bare string literal. CHIRP_PROTOCOL: str = "Chirp" +#: Protocol-type name for the two-pulse steady-state-inactivation voltage-clamp +#: protocol; a member of :data:`VOLTAGE_PROTOCOLS`. Each sweep holds at a +#: different conditioning prepulse voltage and then steps to a fixed test pulse. +#: Referenced by name in the analysis layer (the sweep executor runs the h∞ +#: inactivation analysis and skips the g-V / τ-V analyses only for this +#: protocol) and the UI, so it gets its own constant rather than a bare string +#: literal. +INACTIVATION_PROTOCOL: str = "Inactivation" + #: Current clamp protocol type names. CURRENT_PROTOCOLS: list[str] = [ "Step", @@ -68,6 +77,7 @@ "Step", "Ramp", "Pulse Train", + INACTIVATION_PROTOCOL, ] # Neuron preset names @@ -88,5 +98,6 @@ FI_CURVE: str = "F-I Curve" IV_CURVE: str = "I-V Curve" NA_CHANNEL_ACTIVATION: str = "Na+ Channel Activation" +STEADY_STATE_INACTIVATION: str = "Steady-State Inactivation" FREQUENCY_RESPONSE: str = "Frequency Response" HYPERPOLARIZATION_STEPS: str = "Hyperpolarization Steps" diff --git a/patch_sim/presets/protocols.py b/patch_sim/presets/protocols.py index 63385c09..388d71c4 100644 --- a/patch_sim/presets/protocols.py +++ b/patch_sim/presets/protocols.py @@ -17,9 +17,11 @@ FI_CURVE, FREQUENCY_RESPONSE, HYPERPOLARIZATION_STEPS, + INACTIVATION_PROTOCOL, IV_CURVE, NA_CHANNEL_ACTIVATION, REPETITIVE_FIRING, + STEADY_STATE_INACTIVATION, SUBTHRESHOLD_RESPONSE, VOLTAGE_CLAMP, ) @@ -97,6 +99,24 @@ "max_stimulus": 60.0, "stimulus_step": 10.0, }, + # Two-pulse steady-state inactivation: a long conditioning prepulse (150 ms, + # ample for the Na+ h gate to reach steady state) at each voltage from -120 + # to -20 mV in 10 mV steps, then a brief fixed test pulse to 0 mV (well past + # Na+ activation V½, near peak open probability) whose peak inward current + # measures the available (non-inactivated) channel fraction. The post-test + # tail (30 ms) is generous so the test-pulse transient sits away from the + # right edge of the trace rather than crammed against it. + STEADY_STATE_INACTIVATION: { + "clamp_mode": VOLTAGE_CLAMP, + "protocol_type": INACTIVATION_PROTOCOL, + "pre_stimulus_duration": 150.0, + "stimulus_duration": 15.0, + "post_stimulus_duration": 30.0, + "min_stimulus": -120.0, + "max_stimulus": -20.0, + "stimulus_step": 10.0, + "test_pulse_voltage": 0.0, + }, # Tuned so the chirp response stays subthreshold for the quiescent-at-rest # presets (Squid, FSI, Cortical Pyramidal, CA1) and produces a clean # impedance profile without any user tweaking — amp=0.25 µA/cm² gives a diff --git a/patch_sim/protocols/builders.py b/patch_sim/protocols/builders.py index db47a002..425b1e18 100644 --- a/patch_sim/protocols/builders.py +++ b/patch_sim/protocols/builders.py @@ -201,6 +201,7 @@ def build_voltage_protocol( pulse_amplitude: float = 20.0, pulse_width: float = 2.0, pulse_interval: float = 10.0, + test_pulse_voltage: float = 0.0, min_stimulus: float = 0.0, max_stimulus: float = 0.0, stimulus_step: float = 0.0, @@ -208,28 +209,41 @@ def build_voltage_protocol( """Build voltage clamp protocol arrays from explicit parameters. Args: - protocol_type: One of "Step", "Ramp", or "Pulse Train". + protocol_type: One of "Step", "Ramp", "Pulse Train", or "Inactivation". sampling_frequency: Sampling frequency in Hz. - pre_stimulus_duration: Duration before the stimulus in ms. - stimulus_duration: Duration of the stimulus in ms. + pre_stimulus_duration: Duration before the stimulus in ms. For the + "Inactivation" protocol this is the conditioning prepulse duration. + stimulus_duration: Duration of the stimulus in ms. For the + "Inactivation" protocol this is the fixed test-pulse duration. post_stimulus_duration: Duration after the stimulus in ms. - holding_voltage: Holding voltage in mV. + holding_voltage: Holding voltage in mV. Ignored by the "Inactivation" + protocol, which uses each sweep's conditioning prepulse as the + holding level. start_voltage: Ramp start voltage in mV. end_voltage: Ramp end voltage in mV. pulse_amplitude: Pulse amplitude in mV. pulse_width: Pulse width in ms. pulse_interval: Interval between pulse starts in ms. + test_pulse_voltage: Fixed test-pulse voltage in mV for the + "Inactivation" two-pulse protocol. The swept ``min_stimulus`` / + ``max_stimulus`` / ``stimulus_step`` range then specifies the + conditioning prepulse voltage applied as each sweep's holding level. min_stimulus: Step voltage amplitude in mV for a single-step, or - minimum voltage when running a multi-sweep range. + minimum voltage when running a multi-sweep range. For the + "Inactivation" protocol this is the most-hyperpolarized prepulse. max_stimulus: Maximum voltage for multi-sweep range in mV. Must - be >= min_stimulus. + be >= min_stimulus. For the "Inactivation" protocol this is the + most-depolarized prepulse. stimulus_step: Step size for multi-sweep range in mV. Use 0.0 (or set min_stimulus == max_stimulus) for a single-step protocol. + Must be > 0.0 for the "Inactivation" protocol, which is inherently + multi-sweep. Returns: 2-D array of shape ``(n_sweeps, n_samples)``. Single-step protocols return shape ``(1, n_samples)``. Multi-sweep Step protocols return one - row per voltage level. + row per voltage level; the "Inactivation" protocol returns one row per + conditioning prepulse voltage. Raises: ValueError: If protocol_type is unrecognized or parameters are invalid. @@ -297,6 +311,44 @@ def build_voltage_protocol( sampling_frequency=sampling_frequency, ) ] + elif protocol_type == "Inactivation": + # Two-pulse steady-state-inactivation protocol: each sweep holds at a + # different conditioning prepulse voltage (the swept min/max/step range) + # for the pre-stimulus duration, then steps to the fixed + # test_pulse_voltage for the stimulus duration, then returns to the + # prepulse level. Built from ordinary step_voltage arrays so the + # standard I-V analysis can extract the test-pulse peak current per + # sweep over the [pre, pre + stimulus] window. Unlike "Step", a + # single sweep is rejected — there is no inactivation curve to measure. + if stimulus_step < 0.0: + raise ValueError("stimulus_step must be >= 0.0") + if min_stimulus > max_stimulus: + raise ValueError( + f"min_stimulus ({min_stimulus}) must be" + f" <= max_stimulus ({max_stimulus})" + ) + if min_stimulus == max_stimulus: + raise ValueError( + "Inactivation protocol requires a conditioning prepulse range " + "(min_stimulus must differ from max_stimulus)" + ) + if stimulus_step == 0.0: + raise ValueError( + "stimulus_step must be > 0.0 for the Inactivation protocol" + ) + n_steps = round((max_stimulus - min_stimulus) / stimulus_step) + 1 + prepulse_voltages = np.linspace(min_stimulus, max_stimulus, n_steps) + sweeps = [ + step_voltage( + duration=total_duration, + voltage_amplitude=test_pulse_voltage, + step_start=pre_stimulus_duration, + step_duration=stimulus_duration, + holding_voltage=float(prepulse), + sampling_frequency=sampling_frequency, + ) + for prepulse in prepulse_voltages + ] else: raise ValueError(f"Unknown voltage protocol type: {protocol_type!r}") logger.debug( diff --git a/patch_sim_ui/components/metrics/__init__.py b/patch_sim_ui/components/metrics/__init__.py index 5b8972d3..c8db58c9 100644 --- a/patch_sim_ui/components/metrics/__init__.py +++ b/patch_sim_ui/components/metrics/__init__.py @@ -27,6 +27,7 @@ from patch_sim_ui.components.metrics.calcium_panel import _ca_transient_row from patch_sim_ui.components.metrics.fi_gv_panel import ( _gv_plot, + _inactivation_plot, _iv_curve_tab, _tau_v_plot, ) @@ -224,6 +225,7 @@ def analysis_sidebar() -> rx.Component: "_burst_row", "_ca_transient_row", "_gv_plot", + "_inactivation_plot", "_iv_curve_tab", "_spike_row", "_tau_v_plot", diff --git a/patch_sim_ui/components/metrics/fi_gv_panel.py b/patch_sim_ui/components/metrics/fi_gv_panel.py index 033e2c18..04fffb79 100644 --- a/patch_sim_ui/components/metrics/fi_gv_panel.py +++ b/patch_sim_ui/components/metrics/fi_gv_panel.py @@ -1,9 +1,10 @@ -"""F-I, g-V, τ-V, and I-V curve panels. +"""F-I, g-V, h∞, τ-V, and I-V curve panels. Embedded Plotly figures for current-clamp F-I curves and voltage-clamp I-V -analysis (with optional g-V Boltzmann fit and τ-V activation/inactivation -time constants). The :func:`_iv_curve_tab` composer also reuses the -calcium-transient summary when calcium dynamics were active. +analysis (with optional g-V activation Boltzmann fit, h∞ steady-state +inactivation Boltzmann fit, and τ-V activation/inactivation time constants). +The :func:`_iv_curve_tab` composer also reuses the calcium-transient summary +when calcium dynamics were active. """ import reflex as rx @@ -52,6 +53,28 @@ def _gv_plot() -> rx.Component: ) +def _inactivation_plot() -> rx.Component: + """Render the h∞ steady-state inactivation curve plot. + + Shown only for the two-pulse Inactivation voltage-clamp protocol. + + Returns: + A bordered flex container with the h∞ Plotly figure embedded below + the I-V curve. + """ + return rx.flex( + rx.plotly( + data=AnalysisState.inactivation_figure, + width="100%", + ), + direction="column", + width="100%", + flex_shrink="0", + border_top="1px solid var(--gray-4)", + padding="1", + ) + + def _tau_v_plot() -> rx.Component: """Render the τ-V (activation/inactivation time constants) plot. @@ -78,9 +101,12 @@ def _iv_curve_tab() -> rx.Component: Shows a Plotly I-V curve with peak inward, peak outward, and steady-state traces when voltage clamp multi-sweep data is available. When g-V data is also available, a normalized conductance plot with Boltzmann fit is shown - below the I-V curve. When calcium dynamics were active, the calcium - transient summary is shown above the I-V curve. Displays a placeholder - message when no data exists. + below the I-V curve; when h∞ steady-state inactivation data is available + (the two-pulse Inactivation protocol), a normalized availability plot with + a decreasing Boltzmann fit is shown — in that case the I-V curve above it + is the same data before normalization. When calcium dynamics were active, + the calcium transient summary is shown above the I-V curve. Displays a + placeholder message when no data exists. Returns: The full tab content as a flex column. @@ -107,6 +133,11 @@ def _iv_curve_tab() -> rx.Component: _gv_plot(), rx.box(), ), + rx.cond( + AnalysisState.has_inactivation_data, + _inactivation_plot(), + rx.box(), + ), rx.cond( AnalysisState.has_tau_v_data, _tau_v_plot(), diff --git a/patch_sim_ui/components/protocol_panel.py b/patch_sim_ui/components/protocol_panel.py index 64cdfd1b..8505b03e 100644 --- a/patch_sim_ui/components/protocol_panel.py +++ b/patch_sim_ui/components/protocol_panel.py @@ -136,6 +136,16 @@ class ParamField: ParamField("Pulse interval (ms)", "vc_pulse_interval"), ParamField("Holding voltage (mV)", "holding_voltage"), ), + # Two-pulse steady-state-inactivation protocol. The shared duration fields + # double as the conditioning-prepulse duration ("Pre-stimulus") and the + # test-pulse duration ("Stimulus"); the min/max/step trio sweeps the + # conditioning prepulse voltage and "Test pulse" is the fixed test voltage. + (VOLTAGE_CLAMP, "Inactivation"): ( + ParamField("Prepulse min (mV)", "min_stimulus"), + ParamField("Prepulse max (mV)", "max_stimulus"), + ParamField("Prepulse step (mV)", "stimulus_step"), + ParamField("Test pulse (mV)", "test_pulse_voltage"), + ), } # Single-sweep amplitude fields for the Step protocol. In single-sweep mode diff --git a/patch_sim_ui/plotting/__init__.py b/patch_sim_ui/plotting/__init__.py index b150ab14..17a8485c 100644 --- a/patch_sim_ui/plotting/__init__.py +++ b/patch_sim_ui/plotting/__init__.py @@ -9,6 +9,7 @@ from patch_sim_ui.plotting.gv_curve import build_gv_figure from patch_sim_ui.plotting.hyperpolarization import build_hyperpolarization_figure from patch_sim_ui.plotting.impedance import build_impedance_figure +from patch_sim_ui.plotting.inactivation import build_inactivation_figure from patch_sim_ui.plotting.iv_curve import build_iv_figure from patch_sim_ui.plotting.phase_plane import build_phase_plane_figure from patch_sim_ui.plotting.sfa import build_sfa_figure @@ -26,6 +27,7 @@ "build_gv_figure", "build_hyperpolarization_figure", "build_impedance_figure", + "build_inactivation_figure", "build_iv_figure", "build_phase_plane_figure", "build_sfa_figure", diff --git a/patch_sim_ui/plotting/inactivation.py b/patch_sim_ui/plotting/inactivation.py new file mode 100644 index 00000000..b831acb7 --- /dev/null +++ b/patch_sim_ui/plotting/inactivation.py @@ -0,0 +1,77 @@ +"""Steady-state inactivation (h∞) curve figure with optional Boltzmann fit.""" + +import plotly.graph_objects as go + +from patch_sim_ui.plotting._layout import ANALYSIS_FIGURE_LAYOUT + + +def build_inactivation_figure(inactivation_data: dict) -> go.Figure: + """Build a Plotly h∞ curve figure with optional Boltzmann fit overlay. + + Renders the normalized availability (h∞) data points as a scatter trace. + When the Boltzmann fit converged, a smooth decreasing fit curve is overlaid + as a dashed line and the half-inactivation voltage (V_half) and slope + factor (k) are shown as an annotation. + + Args: + inactivation_data: Dict with keys ``prepulse_voltages``, + ``h_normalized``, ``boltzmann_converged``, ``v_half``, ``k``, + ``fit_voltages``, and ``fit_h_normalized``. + + Returns: + A Plotly :class:`go.Figure` ready to be serialised and sent to the UI. + """ + prepulse_voltages = inactivation_data["prepulse_voltages"] + h_norm = inactivation_data["h_normalized"] + + fig = go.Figure() + fig.add_trace( + go.Scatter( + x=prepulse_voltages, + y=h_norm, + mode="markers", + name="h∞", + marker=dict(color="#8e44ad", size=7), + ) + ) + + if inactivation_data.get("boltzmann_converged"): + v_half = inactivation_data["v_half"] + k = inactivation_data["k"] + fig.add_trace( + go.Scatter( + x=inactivation_data["fit_voltages"], + y=inactivation_data["fit_h_normalized"], + mode="lines", + name="Boltzmann fit", + line=dict(color="#6c3483", dash="dash"), + ) + ) + fig.add_annotation( + xref="paper", + yref="paper", + x=0.02, + y=0.02, + xanchor="left", + yanchor="bottom", + text=f"V₁/₂ = {v_half:.1f} mV
k = {k:.1f} mV", + showarrow=False, + font=dict(size=10), + align="left", + ) + + fig.update_layout( + **ANALYSIS_FIGURE_LAYOUT, + xaxis_title="Prepulse voltage (mV)", + yaxis_title="h∞ (availability)", + showlegend=True, + legend=dict( + orientation="h", + yanchor="bottom", + y=1.02, + xanchor="right", + x=1, + font=dict(size=10), + ), + ) + return fig diff --git a/patch_sim_ui/state/_analysis_format.py b/patch_sim_ui/state/_analysis_format.py index 192a5a35..b744e662 100644 --- a/patch_sim_ui/state/_analysis_format.py +++ b/patch_sim_ui/state/_analysis_format.py @@ -20,8 +20,9 @@ logger = logging.getLogger(__name__) -#: Number of voltage points used for the pre-computed Boltzmann fit curve. -_GV_FIT_POINTS = 200 +#: Number of voltage points used for a pre-computed Boltzmann fit curve +#: (shared by the g-V activation curve and the h∞ inactivation curve). +_BOLTZMANN_FIT_POINTS = 200 def _fmt_optional(value: "float | None", fmt: str) -> str: @@ -839,7 +840,7 @@ def _compute_gv_data( """Compute g-V analysis data from an I-V result and a reversal potential. Calls :func:`patch_sim.compute_gv` to derive normalized conductance and fit - a Boltzmann sigmoid. A dense voltage array (:data:`_GV_FIT_POINTS` points) + a Boltzmann sigmoid. A dense voltage array (:data:`_BOLTZMANN_FIT_POINTS` points) spanning the range of included steps is pre-computed so the plotting function can draw a smooth fit curve without importing scipy. @@ -864,7 +865,7 @@ def _compute_gv_data( if fit.converged and len(gv_result.voltage_steps) >= 2: v_min = min(gv_result.voltage_steps) v_max = max(gv_result.voltage_steps) - v_arr = np.linspace(v_min, v_max, _GV_FIT_POINTS) + v_arr = np.linspace(v_min, v_max, _BOLTZMANN_FIT_POINTS) fit_voltages = v_arr.tolist() fit_gn = [float(patch_sim.boltzmann(v, fit.v_half, fit.k)) for v in v_arr] @@ -878,3 +879,52 @@ def _compute_gv_data( "fit_voltages": fit_voltages, "fit_g_normalized": fit_gn, } + + +def _compute_inactivation_data( + iv_result: "patch_sim.IVAnalysisResult", +) -> "dict[str, Any]": + """Compute steady-state inactivation (h∞) data from an I-V result. + + Calls :func:`patch_sim.compute_inactivation` to derive normalized + availability per conditioning prepulse and fit a decreasing Boltzmann + sigmoid. A dense prepulse-voltage array (:data:`_BOLTZMANN_FIT_POINTS` points) + spanning the range of measured prepulses is pre-computed so the plotting + function can draw a smooth fit curve without importing scipy. Because the + fitted ``k`` is positive but the curve is decreasing, the dense curve is + evaluated as ``boltzmann(V, v_half, -k)``. + + Args: + iv_result: Pre-computed I-V analysis result whose stimulus window is the + two-pulse protocol's fixed test pulse, with one point per + conditioning prepulse voltage. + + Returns: + A dict with keys ``prepulse_voltages``, ``h_normalized``, + ``boltzmann_converged``, ``v_half``, ``k``, ``fit_voltages``, and + ``fit_h_normalized``. Returns an empty dict when *iv_result* has no + points. + """ + result = patch_sim.compute_inactivation(iv_result) + if not result.points: + return {} + + fit = result.boltzmann + fit_voltages: list[float] = [] + fit_h: list[float] = [] + if fit.converged and len(result.prepulse_voltages) >= 2: + v_min = min(result.prepulse_voltages) + v_max = max(result.prepulse_voltages) + v_arr = np.linspace(v_min, v_max, _BOLTZMANN_FIT_POINTS) + fit_voltages = v_arr.tolist() + fit_h = [float(patch_sim.boltzmann(v, fit.v_half, -fit.k)) for v in v_arr] + + return { + "prepulse_voltages": result.prepulse_voltages, + "h_normalized": result.h_normalized_values, + "boltzmann_converged": fit.converged, + "v_half": fit.v_half, + "k": fit.k, + "fit_voltages": fit_voltages, + "fit_h_normalized": fit_h, + } diff --git a/patch_sim_ui/state/_sweep_executor.py b/patch_sim_ui/state/_sweep_executor.py index c95da34f..3de54584 100644 --- a/patch_sim_ui/state/_sweep_executor.py +++ b/patch_sim_ui/state/_sweep_executor.py @@ -15,7 +15,12 @@ import patch_sim import patch_sim.channels -from patch_sim.constants import CHIRP_PROTOCOL, CURRENT_CLAMP, VOLTAGE_CLAMP +from patch_sim.constants import ( + CHIRP_PROTOCOL, + CURRENT_CLAMP, + INACTIVATION_PROTOCOL, + VOLTAGE_CLAMP, +) from patch_sim_ui import constants from patch_sim_ui.api import traces from patch_sim_ui.plotting import TraceVisibility, build_figure @@ -26,6 +31,7 @@ _compute_cc_multi_sweep_analysis, _compute_gv_data, _compute_impedance_data, + _compute_inactivation_data, _compute_iv_data, _compute_multi_sweep_burst_data, _compute_multi_sweep_ca_transient_data, @@ -52,6 +58,7 @@ class _SimResult: iv_data: dict[str, Any] = dataclasses.field(default_factory=dict) gv_data: dict[str, Any] = dataclasses.field(default_factory=dict) tau_v_data: dict[str, Any] = dataclasses.field(default_factory=dict) + inactivation_data: dict[str, Any] = dataclasses.field(default_factory=dict) ap_metrics: list[dict[str, Any]] = dataclasses.field(default_factory=list) ap_summary: dict[str, Any] = dataclasses.field(default_factory=dict) ap_is_multi_sweep: bool = False @@ -155,7 +162,18 @@ def _compute_simulation( pre_stimulus_duration, stimulus_duration, ) - if iv_result is not None: + gv_data: dict[str, Any] = {} + tau_v_data: dict[str, Any] = {} + inactivation_data: dict[str, Any] = {} + if protocol_type == INACTIVATION_PROTOCOL: + # For the two-pulse protocol the swept command voltage is the + # conditioning prepulse, not the test voltage, so the g-V and + # τ-V analyses (which assume the swept voltage *is* the test + # voltage) would be meaningless — skip them and run the + # steady-state inactivation (h∞) analysis instead. + if iv_result is not None: + inactivation_data = _compute_inactivation_data(iv_result) + elif iv_result is not None: # Pick the first *gated* Na⁺ channel — the leak channel "NaL" # also has NernstSpec(SODIUM) but has no gating variables, so # filtering on gating_variables avoids feeding E_NaL into the @@ -171,21 +189,18 @@ def _compute_simulation( ), None, ) - gv_data = ( - _compute_gv_data(iv_result, na_channel.reversal_potential(neuron)) - if na_channel is not None - else {} + if na_channel is not None: + gv_data = _compute_gv_data( + iv_result, na_channel.reversal_potential(neuron) + ) + tau_v_data = _compute_tau_v_data( + new_sweeps, + min_stimulus, + max_stimulus, + stimulus_step, + pre_stimulus_duration, + stimulus_duration, ) - else: - gv_data = {} - tau_v_data = _compute_tau_v_data( - new_sweeps, - min_stimulus, - max_stimulus, - stimulus_step, - pre_stimulus_duration, - stimulus_duration, - ) ms_ca_metrics, ms_ca_summary = _compute_multi_sweep_ca_transient_data( new_sweeps ) @@ -195,6 +210,7 @@ def _compute_simulation( iv_data=iv_data, gv_data=gv_data, tau_v_data=tau_v_data, + inactivation_data=inactivation_data, ca_transient_metrics=ms_ca_metrics, ca_transient_summary=ms_ca_summary, ) diff --git a/patch_sim_ui/state/analysis.py b/patch_sim_ui/state/analysis.py index 429210ff..24d6d0a8 100644 --- a/patch_sim_ui/state/analysis.py +++ b/patch_sim_ui/state/analysis.py @@ -10,6 +10,7 @@ build_gv_figure, build_hyperpolarization_figure, build_impedance_figure, + build_inactivation_figure, build_iv_figure, build_phase_plane_figure, build_sfa_figure, @@ -34,6 +35,8 @@ class AnalysisState(rx.State): fi_data: dict[str, Any] = {} # Serialized FIAnalysisResult for the UI iv_data: dict[str, Any] = {} # Serialized IVAnalysisResult for the UI gv_data: dict[str, Any] = {} # Serialized GVAnalysisResult for the UI + # Serialized InactivationAnalysisResult (h∞ curve) for the UI. + inactivation_data: dict[str, Any] = {} tau_v_data: dict[str, Any] = {} # Serialized TauVAnalysisResult for the UI sfa_data: dict[str, Any] = {} # Serialized SFAAnalysisResult for the UI hyperpolarization_data: dict[str, Any] = {} # Serialized sag/rebound analysis @@ -67,6 +70,7 @@ def clear_results(self) -> None: self.fi_data = {} self.iv_data = {} self.gv_data = {} + self.inactivation_data = {} self.tau_v_data = {} self.sfa_data = {} self.hyperpolarization_data = {} @@ -129,6 +133,11 @@ def has_gv_data(self) -> bool: """Return True when g-V analysis results are available for display.""" return len(self.gv_data) > 0 + @rx.var + def has_inactivation_data(self) -> bool: + """Return True when h∞ steady-state inactivation results are available.""" + return len(self.inactivation_data) > 0 + @rx.var def has_tau_v_data(self) -> bool: """Return True when τ-V analysis results are available for display.""" @@ -164,6 +173,16 @@ def gv_figure(self) -> go.Figure: return go.Figure() return build_gv_figure(self.gv_data) + @rx.var + def inactivation_figure(self) -> go.Figure: + """Return a Plotly h∞ curve figure with decreasing Boltzmann fit overlay. + + Returns an empty figure when no inactivation data is available. + """ + if not self.inactivation_data: + return go.Figure() + return build_inactivation_figure(self.inactivation_data) + @rx.var def tau_v_figure(self) -> go.Figure: """Return a Plotly τ-V figure with activation and inactivation traces. diff --git a/patch_sim_ui/state/protocol.py b/patch_sim_ui/state/protocol.py index dfbc8363..84b808f2 100644 --- a/patch_sim_ui/state/protocol.py +++ b/patch_sim_ui/state/protocol.py @@ -8,7 +8,7 @@ import patch_sim import patch_sim.clamp_simulations -from patch_sim.constants import CURRENT_CLAMP +from patch_sim.constants import CURRENT_CLAMP, INACTIVATION_PROTOCOL from patch_sim.presets import NEURON_PROTOCOL_ADJUSTMENTS, PROTOCOL_PRESETS from patch_sim.protocols.builders import build_current_protocol, build_voltage_protocol from patch_sim_ui import constants @@ -39,6 +39,9 @@ "vc_pulse_amplitude", "vc_pulse_width", "vc_pulse_interval", + # No vc_ prefix: like ``holding_voltage`` it has no current-clamp twin, and + # protocol presets set it by this (builder-matching) name via setattr. + "test_pulse_voltage", ] logger = logging.getLogger(__name__) @@ -93,6 +96,8 @@ class ProtocolState(rx.State): vc_pulse_amplitude: float = 20.0 vc_pulse_width: float = 2.0 vc_pulse_interval: float = 10.0 + # Fixed test-pulse voltage (mV) for the two-pulse "Inactivation" protocol. + test_pulse_voltage: float = 0.0 # ------------------------------------------------------------------ # # Preset tracking # @@ -130,9 +135,13 @@ def is_step_single_sweep(self) -> bool: def can_run_continuous(self) -> bool: """True when the active protocol is compatible with continuous mode. - Multi-sweep Step configurations (min_stimulus != max_stimulus) are - excluded; all other protocols run as a single sweep and are compatible. + Continuous mode streams a single sweep, so inherently multi-sweep + protocols are excluded: multi-sweep Step configurations + (min_stimulus != max_stimulus) and the two-pulse Inactivation protocol. + Every other protocol runs as a single sweep and is compatible. """ + if self.protocol_type == INACTIVATION_PROTOCOL: + return False if self.protocol_type != "Step": return True return self.is_step_single_sweep @@ -437,10 +446,16 @@ def _build_protocols(self) -> "list[tuple[np.ndarray, str]]": pulse_amplitude=self.vc_pulse_amplitude, pulse_width=self.vc_pulse_width, pulse_interval=self.vc_pulse_interval, + test_pulse_voltage=self.test_pulse_voltage, min_stimulus=self.min_stimulus, max_stimulus=self.max_stimulus, stimulus_step=self.stimulus_step, ) if arrays.shape[0] > 1: - return self._attach_step_labels(arrays, "{:+.0f} mV") + label_fmt = ( + "{:+.0f} mV prepulse" + if self.protocol_type == INACTIVATION_PROTOCOL + else "{:+.0f} mV" + ) + return self._attach_step_labels(arrays, label_fmt) return [(arrays[0], "")] diff --git a/patch_sim_ui/state/simulation.py b/patch_sim_ui/state/simulation.py index 8cb08c69..3cc7910e 100644 --- a/patch_sim_ui/state/simulation.py +++ b/patch_sim_ui/state/simulation.py @@ -463,6 +463,7 @@ def _do_apply_simulation( analysis_st.iv_data = result.iv_data analysis_st.gv_data = result.gv_data analysis_st.tau_v_data = result.tau_v_data + analysis_st.inactivation_data = result.inactivation_data analysis_st.ap_metrics = result.ap_metrics analysis_st.ap_summary = result.ap_summary analysis_st.ap_is_multi_sweep = result.ap_is_multi_sweep diff --git a/tests/e2e/test_voltage_clamp_flow.py b/tests/e2e/test_voltage_clamp_flow.py index 8c7e6857..e2cd491f 100644 --- a/tests/e2e/test_voltage_clamp_flow.py +++ b/tests/e2e/test_voltage_clamp_flow.py @@ -8,6 +8,7 @@ IV_CURVE, NA_CHANNEL_ACTIVATION, SQUID_GIANT_AXON, + STEADY_STATE_INACTIVATION, ) from tests.e2e.conftest import StateTree, run_flow @@ -118,3 +119,52 @@ async def test_voltage_clamp_sim_token_is_set(state_tree: StateTree) -> None: ) assert state_tree.sim.sim_token != "" + + +async def test_steady_state_inactivation_preset_produces_multi_sweep( + state_tree: StateTree, +) -> None: + """The Steady-State Inactivation preset produces multiple sweeps.""" + result = await run_flow( + state_tree, + neuron_preset=SQUID_GIANT_AXON, + protocol_preset=STEADY_STATE_INACTIVATION, + ) + + assert len(result.sweeps) > 1 + + +async def test_steady_state_inactivation_preset_populates_inactivation_data( + state_tree: StateTree, +) -> None: + """The Steady-State Inactivation preset populates inactivation_data after run.""" + await run_flow( + state_tree, + neuron_preset=SQUID_GIANT_AXON, + protocol_preset=STEADY_STATE_INACTIVATION, + ) + + inact = state_tree.analysis.inactivation_data + assert inact != {} + assert "prepulse_voltages" in inact + assert "h_normalized" in inact + assert len(inact["prepulse_voltages"]) > 1 + + +async def test_steady_state_inactivation_preset_skips_gv_and_tau_v( + state_tree: StateTree, +) -> None: + """The Inactivation protocol does not populate g-V or τ-V data. + + Those analyses treat the swept command voltage as the test voltage, which + is false for the two-pulse protocol (the swept voltage is the conditioning + prepulse), so the sweep executor skips them. + """ + await run_flow( + state_tree, + neuron_preset=SQUID_GIANT_AXON, + protocol_preset=STEADY_STATE_INACTIVATION, + ) + + assert state_tree.analysis.gv_data == {} + assert state_tree.analysis.tau_v_data == {} diff --git a/tests/integration/test_inactivation_simulation.py b/tests/integration/test_inactivation_simulation.py new file mode 100644 index 00000000..0f9e92ca --- /dev/null +++ b/tests/integration/test_inactivation_simulation.py @@ -0,0 +1,68 @@ +"""Integration tests for patch_sim.analysis.inactivation against real HH sims. + +Runs a two-pulse steady-state-inactivation voltage-clamp protocol through the +full simulation pipeline and verifies that the derived h∞ curve and Boltzmann +fit are physiologically plausible. Unit tests with synthetic data live in +tests/unit/test_inactivation.py. +""" + +import numpy as np + +import patch_sim +from patch_sim.analysis.inactivation import compute_inactivation +from patch_sim.analysis.iv_curve import analyze_iv + + +def test_inactivation_curve_integration_hh_neuron(hh_model): + """h∞ from a real HH simulation is monotone-decreasing with a valid fit. + + Runs a two-pulse inactivation protocol (long conditioning prepulse swept + across voltages, then a fixed depolarizing test pulse), measures the + test-pulse peak inward current per sweep via the I-V analysis, then derives + the h∞ curve. Verifies that availability lies in [0, 1], decreases with + prepulse depolarization, and that the decreasing Boltzmann fit converges + with a physiologically plausible half-inactivation voltage (the classic HH + Na+ h∞ has V½ ≈ −60 mV). + """ + pre_ms, stim_ms = 150.0, 15.0 + min_v, max_v, step_v = -120.0, -20.0, 10.0 + + protocol_2d = patch_sim.build_voltage_protocol( + "Inactivation", + sampling_frequency=40_000, + pre_stimulus_duration=pre_ms, + stimulus_duration=stim_ms, + post_stimulus_duration=5.0, + min_stimulus=min_v, + max_stimulus=max_v, + stimulus_step=step_v, + test_pulse_voltage=0.0, + ) + n_steps = round((max_v - min_v) / step_v) + 1 + prepulses = list(np.linspace(min_v, max_v, n_steps)) + protocols = [protocol_2d[i] for i in range(protocol_2d.shape[0])] + + results = list( + patch_sim.simulate_batch(hh_model, protocols, patch_sim.simulate_voltage_clamp) + ) + time = results[0]["time"] + currents = [r["Itotal"] for r in results] + iv = analyze_iv(time, currents, prepulses, pre_ms, pre_ms + stim_ms) + result = compute_inactivation(iv) + + assert len(result.points) == n_steps + + # Availability is normalized to [0, 1] (allow a hair of float slack). + h = result.h_normalized_values + for value in h: + assert 0.0 <= value <= 1.0 + 1e-9 + + # The most-hyperpolarized prepulse leaves channels available; the + # most-depolarized inactivates them. Allow tiny numerical non-monotonicity. + assert h[0] >= h[-1] + assert np.all(np.diff(h) <= 0.05) + + # The decreasing Boltzmann fit should converge with a plausible V½ and k. + assert result.boltzmann.converged is True + assert -85.0 <= result.boltzmann.v_half <= -40.0 + assert 0.0 < result.boltzmann.k < 25.0 diff --git a/tests/ui/test_plotting.py b/tests/ui/test_plotting.py index abff04b8..a0a9c34e 100644 --- a/tests/ui/test_plotting.py +++ b/tests/ui/test_plotting.py @@ -15,6 +15,7 @@ from patch_sim_ui.plotting import ( TraceVisibility, build_figure, + build_inactivation_figure, build_tau_v_figure, compute_trace_visibility_map, ) @@ -1147,3 +1148,62 @@ def test_build_tau_v_figure_uses_log_y_axis_and_voltage_x_axis(): assert fig.layout.yaxis.type == "log" assert "Voltage" in str(fig.layout.xaxis.title.text) assert "τ" in str(fig.layout.yaxis.title.text) + + +# --------------------------------------------------------------------------- +# build_inactivation_figure +# --------------------------------------------------------------------------- + + +def test_build_inactivation_figure_points_only_when_not_converged(): + """Without a converged fit the figure has only the h∞ scatter trace.""" + data = { + "prepulse_voltages": [-100.0, -60.0, -20.0], + "h_normalized": [1.0, 0.5, 0.05], + "boltzmann_converged": False, + "v_half": 0.0, + "k": 1.0, + "fit_voltages": [], + "fit_h_normalized": [], + } + fig = build_inactivation_figure(data) + assert isinstance(fig, go.Figure) + assert len(fig.data) == 1 + assert "h" in str(fig.data[0].name) + assert len(fig.layout.annotations) == 0 + + +def test_build_inactivation_figure_includes_fit_trace_and_annotation(): + """A converged fit adds a dashed Boltzmann trace and a V½ / k annotation.""" + data = { + "prepulse_voltages": [-100.0, -60.0, -20.0], + "h_normalized": [1.0, 0.5, 0.05], + "boltzmann_converged": True, + "v_half": -62.0, + "k": 7.5, + "fit_voltages": [-100.0, -60.0, -20.0], + "fit_h_normalized": [0.99, 0.5, 0.01], + } + fig = build_inactivation_figure(data) + assert len(fig.data) == 2 + trace_names = {trace.name for trace in fig.data} + assert "Boltzmann fit" in trace_names + assert len(fig.layout.annotations) == 1 + text = fig.layout.annotations[0].text + assert "V" in text and "k" in text + + +def test_build_inactivation_figure_axis_titles(): + """The h∞ figure labels its axes Prepulse voltage (mV) and h∞.""" + data = { + "prepulse_voltages": [-100.0, -20.0], + "h_normalized": [1.0, 0.05], + "boltzmann_converged": False, + "v_half": 0.0, + "k": 1.0, + "fit_voltages": [], + "fit_h_normalized": [], + } + fig = build_inactivation_figure(data) + assert "Prepulse voltage" in str(fig.layout.xaxis.title.text) + assert "h" in str(fig.layout.yaxis.title.text) diff --git a/tests/ui/test_state.py b/tests/ui/test_state.py index de4cd4cd..474221ab 100644 --- a/tests/ui/test_state.py +++ b/tests/ui/test_state.py @@ -29,6 +29,7 @@ DOPAMINERGIC, FAST_SPIKING_INTERNEURON, HYPERPOLARIZATION_STEPS, + INACTIVATION_PROTOCOL, NA_CHANNEL_ACTIVATION, PURKINJE, REPETITIVE_FIRING, @@ -809,6 +810,13 @@ def test_can_run_continuous_true_for_ramp() -> None: assert ps.can_run_continuous is True +def test_can_run_continuous_false_for_inactivation() -> None: + """can_run_continuous is False for the multi-sweep Inactivation protocol.""" + ps = _make_protocol_state() + ps.protocol_type = INACTIVATION_PROTOCOL + assert ps.can_run_continuous is False + + # --------------------------------------------------------------------------- # Stimulus range setters — constraint logic # --------------------------------------------------------------------------- @@ -1738,6 +1746,7 @@ def _prime_analysis(an_st: AnalysisState) -> None: an_st.fi_data = {"current_steps": [0.1]} an_st.iv_data = {"voltage_steps": [-70.0]} an_st.gv_data = {"conductance": [0.5]} + an_st.inactivation_data = {"h_normalized": [1.0, 0.5]} an_st.tau_v_data = {"voltages": [-40.0, -20.0]} an_st.sfa_data = {"isi": [10.0]} an_st.hyperpolarization_data = {"sag_ratio": 0.1} @@ -1760,6 +1769,7 @@ def _assert_analysis_cleared(an_st: AnalysisState) -> None: assert an_st.fi_data == {} assert an_st.iv_data == {} assert an_st.gv_data == {} + assert an_st.inactivation_data == {} assert an_st.tau_v_data == {} assert an_st.sfa_data == {} assert an_st.hyperpolarization_data == {} diff --git a/tests/unit/test_inactivation.py b/tests/unit/test_inactivation.py new file mode 100644 index 00000000..07517d5d --- /dev/null +++ b/tests/unit/test_inactivation.py @@ -0,0 +1,141 @@ +"""Unit tests for patch_sim.analysis.inactivation. + +Covers compute_inactivation() with synthetic I-V data and edge cases. +Integration tests against real HH simulations live in +tests/integration/test_inactivation_simulation.py. +""" + +import numpy as np +import pytest + +from patch_sim.analysis.gv_curve import boltzmann +from patch_sim.analysis.inactivation import compute_inactivation +from patch_sim.analysis.iv_curve import IVAnalysisResult, IVPoint + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + + +def _make_iv_result( + prepulse_voltages: list[float], + peak_inward_currents: list[float], +) -> IVAnalysisResult: + """Build a synthetic IVAnalysisResult for testing compute_inactivation. + + Steady-state and peak-outward values are set to zero since they are not + used by compute_inactivation. + + Args: + prepulse_voltages: Conditioning prepulse voltages in mV. + peak_inward_currents: Corresponding test-pulse peak inward currents + (µA/cm²). + + Returns: + An IVAnalysisResult with one IVPoint per prepulse/current pair. + """ + points = [ + IVPoint( + voltage_step=v, + peak_inward_current=i, + peak_outward_current=0.0, + steady_state_current=0.0, + ) + for v, i in zip(prepulse_voltages, peak_inward_currents) + ] + return IVAnalysisResult(points=points) + + +# --------------------------------------------------------------------------- +# compute_inactivation() — unit tests +# --------------------------------------------------------------------------- + + +def test_compute_inactivation_normalization(): + """h∞ values are I_peak / I_peak_max (normalized by the most-negative peak).""" + prepulses = [-100.0, -80.0, -60.0, -40.0] + peaks = [-100.0, -100.0, -50.0, -10.0] + result = compute_inactivation(_make_iv_result(prepulses, peaks)) + assert result.h_normalized_values == pytest.approx([1.0, 1.0, 0.5, 0.1]) + + +def test_compute_inactivation_sorted_by_prepulse(): + """Output points are sorted by ascending prepulse voltage.""" + # Supplied out of order; -100 has the largest inward current. + prepulses = [-40.0, -100.0, -60.0] + peaks = [-10.0, -100.0, -50.0] + result = compute_inactivation(_make_iv_result(prepulses, peaks)) + assert result.prepulse_voltages == [-100.0, -60.0, -40.0] + assert result.h_normalized_values == pytest.approx([1.0, 0.5, 0.1]) + + +def test_compute_inactivation_clamps_to_unit_range(): + """A non-negative test-pulse peak (full inactivation) maps to h∞ = 0.0.""" + prepulses = [-100.0, -60.0, -20.0] + # Last sweep is net outward at the test pulse → ratio would be negative. + peaks = [-200.0, -80.0, 5.0] + result = compute_inactivation(_make_iv_result(prepulses, peaks)) + h = result.h_normalized_values + assert h[0] == pytest.approx(1.0) + assert h[1] == pytest.approx(0.4) + assert h[2] == 0.0 + for value in h: + assert 0.0 <= value <= 1.0 + + +def test_compute_inactivation_boltzmann_recovers_params(): + """The decreasing-Boltzmann fit recovers v_half and k from a known sigmoid.""" + v_half_true = -65.0 + k_true = 8.0 + i_max = -120.0 # µA/cm² at full availability + prepulses = np.linspace(-120.0, -20.0, 11) + # I_peak(V) = i_max * h∞(V), with h∞(V) = boltzmann(V, v_half, -k). + peaks = [float(i_max * boltzmann(v, v_half_true, -k_true)) for v in prepulses] + result = compute_inactivation(_make_iv_result(prepulses.tolist(), peaks)) + assert result.boltzmann.converged is True + assert result.boltzmann.v_half == pytest.approx(v_half_true, abs=1.5) + assert result.boltzmann.k == pytest.approx(k_true, abs=1.5) + + +def test_compute_inactivation_curve_is_decreasing(): + """For a monotone synthetic input the h∞ values are non-increasing.""" + prepulses = np.linspace(-120.0, -20.0, 11) + peaks = [float(-100.0 * boltzmann(v, -60.0, -7.0)) for v in prepulses] + result = compute_inactivation(_make_iv_result(prepulses.tolist(), peaks)) + h = result.h_normalized_values + assert all(h[i + 1] <= h[i] + 1e-9 for i in range(len(h) - 1)) + + +def test_compute_inactivation_empty_iv_result(): + """An empty IVAnalysisResult produces an empty result with no fit.""" + result = compute_inactivation(IVAnalysisResult(points=[])) + assert result.points == [] + assert result.boltzmann.converged is False + + +def test_compute_inactivation_single_point(): + """A single point yields one record but no Boltzmann fit (< 2 points).""" + result = compute_inactivation(_make_iv_result([-80.0], [-50.0])) + assert len(result.points) == 1 + assert result.h_normalized_values == pytest.approx([1.0]) + assert result.boltzmann.converged is False + + +def test_compute_inactivation_all_non_negative_peaks(): + """No inward current at any prepulse → all h∞ = 0.0 and no fit.""" + prepulses = [-100.0, -60.0, -20.0] + peaks = [0.0, 0.0, 1.0] + result = compute_inactivation(_make_iv_result(prepulses, peaks)) + assert result.h_normalized_values == [0.0, 0.0, 0.0] + assert result.boltzmann.converged is False + + +def test_compute_inactivation_convenience_properties_match_points(): + """Convenience list properties match the underlying points list.""" + prepulses = [-100.0, -60.0, -20.0] + peaks = [-100.0, -40.0, -10.0] + result = compute_inactivation(_make_iv_result(prepulses, peaks)) + for i, pt in enumerate(result.points): + assert result.prepulse_voltages[i] == pt.prepulse_voltage + assert result.peak_inward_currents[i] == pt.peak_inward_current + assert result.h_normalized_values[i] == pt.h_normalized diff --git a/tests/unit/test_protocol_builders.py b/tests/unit/test_protocol_builders.py index 31cc60e3..67a31c1b 100644 --- a/tests/unit/test_protocol_builders.py +++ b/tests/unit/test_protocol_builders.py @@ -61,6 +61,102 @@ def test_voltage_protocol_types(protocol_type: str) -> None: assert result.shape[0] == 1 +# --------------------------------------------------------------------------- +# build_voltage_protocol — "Inactivation" two-pulse protocol +# --------------------------------------------------------------------------- + + +def _inactivation_protocol() -> tuple[np.ndarray, float, float, float, list[float]]: + """Build a small Inactivation protocol and return it with its parameters. + + Returns: + A 5-tuple ``(arrays, pre_ms, stim_ms, post_ms, prepulses)`` where + ``arrays`` has shape ``(5, n_samples)``, the durations are in ms and + ``prepulses`` lists the five conditioning prepulse voltages in mV. + """ + pre_ms, stim_ms, post_ms = 50.0, 10.0, 5.0 + min_v, max_v, step_v = -100.0, -20.0, 20.0 + arrays = build_voltage_protocol( + protocol_type="Inactivation", + sampling_frequency=SAMPLING_FREQUENCY, + pre_stimulus_duration=pre_ms, + stimulus_duration=stim_ms, + post_stimulus_duration=post_ms, + min_stimulus=min_v, + max_stimulus=max_v, + stimulus_step=step_v, + test_pulse_voltage=0.0, + ) + n_steps = round((max_v - min_v) / step_v) + 1 + return arrays, pre_ms, stim_ms, post_ms, list(np.linspace(min_v, max_v, n_steps)) + + +def test_build_voltage_protocol_inactivation_shape() -> None: + """The Inactivation protocol returns one sweep per conditioning prepulse.""" + arrays, *_ = _inactivation_protocol() + assert _is_valid_protocol_array(arrays) + assert arrays.shape[0] == 5 + + +def test_build_voltage_protocol_inactivation_holding_varies_step_fixed() -> None: + """Each sweep holds at its prepulse and steps to the fixed test pulse.""" + arrays, pre_ms, stim_ms, _post_ms, prepulses = _inactivation_protocol() + # Sample indices comfortably inside the prepulse, test-pulse and tail. + i_prepulse = int(0.5 * pre_ms * 1e-3 * SAMPLING_FREQUENCY) + i_test = int((pre_ms + 0.5 * stim_ms) * 1e-3 * SAMPLING_FREQUENCY) + i_tail = int((pre_ms + stim_ms + 1.0) * 1e-3 * SAMPLING_FREQUENCY) + for row, prepulse in zip(arrays, prepulses): + assert row[i_prepulse] == pytest.approx(prepulse) + assert row[i_test] == pytest.approx(0.0) + assert row[i_tail] == pytest.approx(prepulse) + + +def test_build_voltage_protocol_inactivation_requires_prepulse_range() -> None: + """A missing prepulse range (min == max) is rejected.""" + with pytest.raises(ValueError, match="prepulse range"): + build_voltage_protocol( + protocol_type="Inactivation", + sampling_frequency=SAMPLING_FREQUENCY, + min_stimulus=-80.0, + max_stimulus=-80.0, + stimulus_step=10.0, + ) + + +def test_build_voltage_protocol_inactivation_requires_positive_step() -> None: + """A zero step over a prepulse range is rejected.""" + with pytest.raises(ValueError, match="stimulus_step"): + build_voltage_protocol( + protocol_type="Inactivation", + sampling_frequency=SAMPLING_FREQUENCY, + min_stimulus=-100.0, + max_stimulus=-20.0, + stimulus_step=0.0, + ) + + +def test_build_voltage_protocol_inactivation_rejects_min_gt_max() -> None: + """A prepulse range with min > max is rejected.""" + with pytest.raises(ValueError, match="min_stimulus"): + build_voltage_protocol( + protocol_type="Inactivation", + sampling_frequency=SAMPLING_FREQUENCY, + min_stimulus=-20.0, + max_stimulus=-100.0, + stimulus_step=10.0, + ) + + +def test_build_protocol_from_preset_steady_state_inactivation() -> None: + """The Steady-State Inactivation preset builds a valid multi-sweep protocol.""" + result = patch_sim.build_protocol_from_preset( + patch_sim.constants.STEADY_STATE_INACTIVATION, + sampling_frequency=SAMPLING_FREQUENCY, + ) + assert _is_valid_protocol_array(result) + assert result.shape[0] > 1 + + # --------------------------------------------------------------------------- # build_protocol_from_preset — exported from patch_sim top-level # ---------------------------------------------------------------------------