From 8e2c26354be983a1e7e948f29defe367db27a1f2 Mon Sep 17 00:00:00 2001 From: JCorson Date: Mon, 11 May 2026 15:19:20 -0500 Subject: [PATCH 1/9] feat(analysis): add impedance-profile analysis from chirp current-clamp runs MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds patch_sim.analyze_impedance / ImpedanceProfile: an FFT-based membrane impedance Z(f) = V̂(f) / Î(f) computed over the swept band of a chirp current-clamp stimulus, reporting |Z(f)| (kΩ·cm², plus absolute MΩ when the neuron declares a surface area), ∠Z(f) in degrees, the resonance frequency f_R (interior peak of |Z|), and the quality factor Q = f_R / FWHM. Closes #185 (analysis layer; UI surfacing follows). --- patch_sim/__init__.py | 4 + patch_sim/analysis/__init__.py | 4 + patch_sim/analysis/impedance.py | 310 ++++++++++++++++++ .../integration/test_impedance_simulation.py | 82 +++++ tests/unit/test_impedance.py | 237 +++++++++++++ 5 files changed, 637 insertions(+) create mode 100644 patch_sim/analysis/impedance.py create mode 100644 tests/integration/test_impedance_simulation.py create mode 100644 tests/unit/test_impedance.py diff --git a/patch_sim/__init__.py b/patch_sim/__init__.py index 1aeb9821..b1dd0253 100644 --- a/patch_sim/__init__.py +++ b/patch_sim/__init__.py @@ -22,6 +22,7 @@ GVAnalysisResult, GVPoint, HyperpolarizationAnalysisResult, + ImpedanceProfile, IVAnalysisResult, IVPoint, PassiveProperties, @@ -39,6 +40,7 @@ analyze_calcium_transients_from_result, analyze_fi, analyze_hyperpolarization, + analyze_impedance, analyze_iv, analyze_passive_properties, analyze_sfa, @@ -118,6 +120,7 @@ "GVAnalysisResult", "GVPoint", "HyperpolarizationAnalysisResult", + "ImpedanceProfile", "IVAnalysisResult", "IVPoint", "PassiveProperties", @@ -135,6 +138,7 @@ "analyze_calcium_transients_from_result", "analyze_fi", "analyze_hyperpolarization", + "analyze_impedance", "analyze_iv", "analyze_passive_properties", "analyze_sfa", diff --git a/patch_sim/analysis/__init__.py b/patch_sim/analysis/__init__.py index 7bfb21cc..9a5d81a3 100644 --- a/patch_sim/analysis/__init__.py +++ b/patch_sim/analysis/__init__.py @@ -9,6 +9,7 @@ derivatives: General-purpose derivative utilities (e.g. dV/dt). 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. 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ₘ). @@ -55,6 +56,7 @@ analyze_hyperpolarization, compute_sag_point, ) +from .impedance import ImpedanceProfile, analyze_impedance from .iv_curve import IVAnalysisResult, IVPoint, analyze_iv, compute_iv_point from .membrane_test import ( MEMBRANE_TEST_CURRENT, @@ -92,6 +94,7 @@ "analyze_calcium_transients_from_result", "analyze_fi", "analyze_hyperpolarization", + "analyze_impedance", "analyze_iv", "analyze_passive_properties", "analyze_sfa", @@ -123,6 +126,7 @@ "FIAnalysisResult", "FIPoint", "HyperpolarizationAnalysisResult", + "ImpedanceProfile", "SagPoint", "GVAnalysisResult", "GVPoint", diff --git a/patch_sim/analysis/impedance.py b/patch_sim/analysis/impedance.py new file mode 100644 index 00000000..5cb25415 --- /dev/null +++ b/patch_sim/analysis/impedance.py @@ -0,0 +1,310 @@ +"""Membrane impedance profile from a chirp current-clamp run. + +Provides :func:`analyze_impedance` for computing the FFT-based membrane +impedance ``Z(f) = V̂(f) / Î(f)`` over the swept frequency band of a chirp +(linear frequency-sweep) current-clamp stimulus. Neurons with Ih/HCN +channels show subthreshold resonance — an interior peak in ``|Z(f)|`` — +which this analysis quantifies via the resonance frequency ``f_R`` and the +quality factor ``Q``. + +Data classes: + ImpedanceProfile: Impedance magnitude/phase spectra and resonance metrics. +""" + +from __future__ import annotations + +import dataclasses +import logging + +import numpy as np + +from patch_sim.analysis.passive_properties import density_to_absolute_r_in + +logger = logging.getLogger(__name__) + +#: Minimum chirp-window duration (ms) required for a meaningful spectrum. A +#: shorter window gives a frequency resolution too coarse to resolve a +#: subthreshold resonance peak. +_MIN_WINDOW_MS: float = 50.0 + +#: Minimum number of FFT bins that must fall inside the analysis band for the +#: result to be considered meaningful. +_MIN_BAND_BINS: int = 8 + +#: Fractional margin by which the ``|Z|`` peak must exceed the low-frequency +#: band edge to count as a genuine resonance (rather than measurement ripple +#: on an otherwise monotone passive response). +_RESONANCE_REL_MARGIN: float = 0.01 + + +@dataclasses.dataclass +class ImpedanceProfile: + """Membrane impedance profile from a chirp current-clamp run. + + The per-area spectra (``magnitude``) are always populated. When + ``area_cm2`` is supplied to :func:`analyze_impedance`, the absolute + counterparts (``magnitude_mohm``, ``peak_impedance_mohm``) are filled in + too — analogous to the per-area / absolute split in + :class:`~patch_sim.analysis.passive_properties.PassiveProperties`. + + ``resonance_frequency``, ``peak_impedance`` and ``quality_factor`` are all + ``None`` when ``|Z(f)|`` has no interior peak that rises meaningfully above + the low-frequency band edge — the hallmark of a passive low-pass cell with + no resonance. ``quality_factor`` may additionally be ``None`` when an + interior peak exists but its half-power width cannot be bracketed inside + the analysis band. + + Attributes: + frequencies: Analysis-band frequency axis in Hz, ascending. + magnitude: ``|Z(f)|`` in kΩ·cm² (``|V̂ / Î|`` where V̂ is in mV and + Î is in µA/cm²), aligned to ``frequencies``. + phase: ``∠Z(f)`` in degrees, aligned to ``frequencies``; positive + values mean the voltage leads the current. + resonance_frequency: Frequency of the ``|Z|`` peak in Hz, or ``None`` + when the maximum lies at a band edge or does not rise meaningfully + above the low-frequency edge (no genuine resonance). + peak_impedance: ``|Z|`` at ``resonance_frequency`` in kΩ·cm², or + ``None`` when there is no interior peak. + quality_factor: Dimensionless ``Q = f_R / FWHM`` where FWHM is the + full width of the ``|Z|`` peak at the half-power (−3 dB) level, or + ``None`` when there is no interior peak or the half-power crossings + cannot be bracketed inside the band. + magnitude_mohm: ``|Z(f)|`` converted to absolute MΩ when ``area_cm2`` + is supplied; ``None`` otherwise. + peak_impedance_mohm: Absolute peak impedance in MΩ when ``area_cm2`` is + supplied and an interior peak exists; ``None`` otherwise. + area_cm2: Membrane surface area in cm² used for the absolute + conversion. ``None`` when only per-area values were requested. + f_start: Lower edge of the analysis band in Hz. + f_end: Upper edge of the analysis band in Hz. + """ + + frequencies: list[float] + magnitude: list[float] + phase: list[float] + resonance_frequency: float | None + peak_impedance: float | None + quality_factor: float | None + magnitude_mohm: list[float] | None = None + peak_impedance_mohm: float | None = None + area_cm2: float | None = None + f_start: float = 0.0 + f_end: float = 0.0 + + +def _half_power_crossing( + frequencies: np.ndarray, + magnitude: np.ndarray, + peak_idx: int, + half_power: float, + direction: int, +) -> float | None: + """Find the frequency where ``|Z|`` falls to the half-power level. + + Walks outward from ``peak_idx`` in ``direction`` (``-1`` toward lower + frequencies, ``+1`` toward higher) to the first bin at or below + ``half_power``, then linearly interpolates the crossing frequency between + that bin and its inward neighbour. + + Args: + frequencies: Analysis-band frequency axis in Hz, ascending. + magnitude: ``|Z(f)|`` aligned to ``frequencies``. + peak_idx: Index of the ``|Z|`` peak within the band. + half_power: Magnitude threshold (peak / sqrt(2)). + direction: ``-1`` to search toward lower bins, ``+1`` toward higher. + + Returns: + The interpolated crossing frequency in Hz, or ``None`` when the + threshold is never reached before the band edge. + """ + idx = peak_idx + n = magnitude.size + while 0 <= idx + direction < n: + nxt = idx + direction + if magnitude[nxt] <= half_power: + m0, m1 = float(magnitude[idx]), float(magnitude[nxt]) + f0, f1 = float(frequencies[idx]), float(frequencies[nxt]) + if m0 == m1: + return f1 + return f0 + (half_power - m0) * (f1 - f0) / (m1 - m0) + idx = nxt + return None + + +def _quality_factor( + frequencies: np.ndarray, + magnitude: np.ndarray, + peak_idx: int, +) -> float | None: + """Estimate the resonance quality factor ``Q = f_R / FWHM``. + + FWHM is the full width of the ``|Z|`` peak at the half-power (−3 dB) level, + i.e. where ``|Z|`` equals ``peak / sqrt(2)``. The half-power crossings on + each side of the peak are found by linear interpolation between bracketing + FFT bins. + + Args: + frequencies: Analysis-band frequency axis in Hz, ascending. + magnitude: ``|Z(f)|`` aligned to ``frequencies``. + peak_idx: Index of the ``|Z|`` peak within the band. + + Returns: + The quality factor, or ``None`` when either half-power crossing cannot + be bracketed inside the band or the bracketed width is non-positive. + """ + half_power = float(magnitude[peak_idx]) / np.sqrt(2.0) + f_low = _half_power_crossing(frequencies, magnitude, peak_idx, half_power, -1) + f_high = _half_power_crossing(frequencies, magnitude, peak_idx, half_power, +1) + if f_low is None or f_high is None: + return None + fwhm = f_high - f_low + if fwhm <= 0.0: + return None + return float(frequencies[peak_idx]) / fwhm + + +def analyze_impedance( + time: np.ndarray, + voltage: np.ndarray, + injected_current: np.ndarray, + stim_start_ms: float, + stim_end_ms: float, + f_start: float, + f_end: float, + area_cm2: float | None = None, +) -> ImpedanceProfile | None: + """Compute the membrane impedance profile from a chirp current-clamp run. + + Extracts the chirp window from the traces, removes the DC component from + both signals, takes the real FFT of each (``numpy.fft.rfft``), forms + ``Z(f) = V̂(f) / Î(f)`` over the swept band ``[f_start, f_end]``, and + reports the magnitude and phase spectra plus the resonance frequency + (argmax of ``|Z|``, only when it is an interior maximum that rises above + the low-frequency band edge) and quality factor. + + Args: + time: Time axis in ms (uniformly sampled). + voltage: Membrane voltage response in mV. + injected_current: Injected chirp current in µA/cm². + stim_start_ms: Start of the chirp stimulus window in ms. + stim_end_ms: End of the chirp stimulus window in ms. + f_start: Starting frequency of the chirp sweep in Hz. + f_end: Ending frequency of the chirp sweep in Hz. + area_cm2: Optional membrane surface area in cm². When supplied, the + absolute MΩ counterparts (``magnitude_mohm``, + ``peak_impedance_mohm``) are populated. + + Returns: + An :class:`ImpedanceProfile`, or ``None`` when the inputs are + inconsistent (mismatched lengths, non-finite values), the band is + invalid (``f_end <= f_start`` or ``f_start < 0``), the chirp window is + shorter than :data:`_MIN_WINDOW_MS`, fewer than :data:`_MIN_BAND_BINS` + FFT bins fall inside the band, or the stimulus has no spectral power in + the band. + """ + t = np.asarray(time, dtype=float) + v = np.asarray(voltage, dtype=float) + i = np.asarray(injected_current, dtype=float) + if v.shape != t.shape or i.shape != t.shape: + logger.warning("analyze_impedance: trace length mismatch; skipping.") + return None + if not (np.all(np.isfinite(v)) and np.all(np.isfinite(i))): + logger.warning("analyze_impedance: non-finite trace values; skipping.") + return None + if f_end <= f_start or f_start < 0.0: + logger.warning( + "analyze_impedance: invalid band [%s, %s] Hz; skipping.", f_start, f_end + ) + return None + if stim_end_ms - stim_start_ms < _MIN_WINDOW_MS: + logger.warning( + "analyze_impedance: chirp window %.1f ms shorter than %.1f ms; skipping.", + stim_end_ms - stim_start_ms, + _MIN_WINDOW_MS, + ) + return None + + mask = (t >= stim_start_ms) & (t < stim_end_ms) + if int(mask.sum()) < 4: + logger.warning("analyze_impedance: too few samples in chirp window; skipping.") + return None + t_win, v_win, i_win = t[mask], v[mask], i[mask] + + dt_ms = float(np.mean(np.diff(t_win))) + if dt_ms <= 0.0: + logger.warning("analyze_impedance: non-positive sampling interval; skipping.") + return None + fs_hz = 1000.0 / dt_ms + + # Remove the DC component from both signals before transforming. No + # additional taper is applied: a linear chirp already starts at zero + # phase, and its magnitude spectrum is approximately flat across the swept + # band, so a Hann (or similar) window is counter-productive here — it + # heavily attenuates the chirp's low- and high-frequency content (which + # lives in its first and last cycles) and inflates the band-edge impedance + # estimate where the windowed stimulus energy collapses toward zero. + v_d = v_win - v_win.mean() + i_d = i_win - i_win.mean() + n = v_d.size + + v_fft = np.fft.rfft(v_d) + i_fft = np.fft.rfft(i_d) + freqs = np.fft.rfftfreq(n, d=1.0 / fs_hz) + + band = (freqs >= f_start) & (freqs <= f_end) + if int(band.sum()) < _MIN_BAND_BINS: + logger.warning( + "analyze_impedance: only %d FFT bins in band; need %d; skipping.", + int(band.sum()), + _MIN_BAND_BINS, + ) + return None + fb = freqs[band] + v_band = v_fft[band] + i_band = i_fft[band] + + if float(np.max(np.abs(i_band))) <= 0.0: + logger.warning("analyze_impedance: stimulus has no power in band; skipping.") + return None + + z = v_band / i_band + mag = np.abs(z) + phase_deg = np.degrees(np.angle(z)) + + peak_idx = int(np.argmax(mag)) + low_edge = float(mag[0]) + is_resonant = 0 < peak_idx < mag.size - 1 and float(mag[peak_idx]) > low_edge * ( + 1.0 + _RESONANCE_REL_MARGIN + ) + if is_resonant: + resonance_frequency: float | None = float(fb[peak_idx]) + peak_impedance: float | None = float(mag[peak_idx]) + quality_factor: float | None = _quality_factor(fb, mag, peak_idx) + else: + resonance_frequency = None + peak_impedance = None + quality_factor = None + + if area_cm2 is not None and area_cm2 > 0.0: + magnitude_mohm: list[float] | None = [float(m) / area_cm2 / 1000.0 for m in mag] + else: + magnitude_mohm = None + peak_impedance_mohm = ( + density_to_absolute_r_in(peak_impedance, area_cm2) + if peak_impedance is not None + else None + ) + + return ImpedanceProfile( + frequencies=fb.tolist(), + magnitude=mag.tolist(), + phase=phase_deg.tolist(), + resonance_frequency=resonance_frequency, + peak_impedance=peak_impedance, + quality_factor=quality_factor, + magnitude_mohm=magnitude_mohm, + peak_impedance_mohm=peak_impedance_mohm, + area_cm2=area_cm2, + f_start=float(f_start), + f_end=float(f_end), + ) diff --git a/tests/integration/test_impedance_simulation.py b/tests/integration/test_impedance_simulation.py new file mode 100644 index 00000000..95db89f0 --- /dev/null +++ b/tests/integration/test_impedance_simulation.py @@ -0,0 +1,82 @@ +"""Integration tests for patch_sim.analysis.impedance against real simulations. + +Runs a chirp current-clamp protocol through the full simulation pipeline and +verifies that the resulting impedance profile is structurally sound and +physiologically scaled. +Unit tests with synthetic signals live in tests/unit/test_impedance.py. +""" + +import numpy as np + +import patch_sim +from patch_sim.clamp_simulations import SIM_SAMPLING_FREQ +from patch_sim.presets import make_thalamic_relay + +_DURATION_MS = 500.0 +_F_START = 1.0 +_F_END = 100.0 + + +def _run_chirp_profile(neuron, amplitude: float): + """Run a chirp current-clamp sweep and return its impedance profile + voltage. + + Args: + neuron: Neuron instance to simulate. + amplitude: Chirp amplitude in µA/cm². + + Returns: + Tuple of (ImpedanceProfile or None, voltage trace array in mV). + """ + stimulus = patch_sim.chirp_current( + duration=_DURATION_MS, + dc_offset=0.0, + amplitude=amplitude, + start_frequency=_F_START, + end_frequency=_F_END, + sampling_frequency=SIM_SAMPLING_FREQ, + ) + result = patch_sim.simulate_current_clamp(neuron, stimulus) + time = np.asarray(result["time"]) + voltage = np.asarray(result["voltage"]) + profile = patch_sim.analyze_impedance( + time, voltage, stimulus, 0.0, _DURATION_MS, _F_START, _F_END + ) + return profile, voltage + + +def test_chirp_pipeline_produces_sane_impedance_profile(hh_model) -> None: + """A subthreshold chirp on the squid axon yields a well-formed impedance profile. + + Checks the structural invariants of the returned profile: ascending + in-band frequency axis, finite positive magnitudes, bounded phase, and + consistent array lengths. + """ + profile, voltage = _run_chirp_profile(hh_model, amplitude=1.0) + + assert voltage.max() < -20.0 # stayed subthreshold + assert profile is not None + freqs = np.asarray(profile.frequencies) + mag = np.asarray(profile.magnitude) + phase = np.asarray(profile.phase) + assert len(freqs) == len(mag) == len(phase) + assert len(freqs) >= 8 + assert np.all(np.diff(freqs) > 0.0) + assert freqs.min() >= _F_START + assert freqs.max() <= _F_END + assert np.all(np.isfinite(mag)) and np.all(mag > 0.0) + assert np.all(np.abs(phase) <= 180.0) + + +def test_thalamic_relay_has_higher_low_frequency_impedance_than_squid( + hh_model, +) -> None: + """A thalamic relay neuron shows much larger low-frequency |Z| than the squid axon. + + Reflects the far higher input resistance of a small thalamic relay cell + compared with the large, low-resistance squid giant axon. + """ + relay_profile, _ = _run_chirp_profile(make_thalamic_relay(), amplitude=1.0) + squid_profile, _ = _run_chirp_profile(hh_model, amplitude=1.0) + + assert relay_profile is not None and squid_profile is not None + assert relay_profile.magnitude[0] > 5.0 * squid_profile.magnitude[0] diff --git a/tests/unit/test_impedance.py b/tests/unit/test_impedance.py new file mode 100644 index 00000000..b1c8920c --- /dev/null +++ b/tests/unit/test_impedance.py @@ -0,0 +1,237 @@ +"""Unit tests for patch_sim.analysis.impedance. + +Covers analyze_impedance with synthetic stimulus/response signals and edge +cases. Integration tests against real simulations live in +tests/integration/test_impedance_simulation.py. +""" + +import numpy as np +import pytest + +import patch_sim +from patch_sim.analysis.impedance import analyze_impedance + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +_FS_HZ = 1000.0 # Sampling rate for the synthetic signals. +_DT_MS = 1000.0 / _FS_HZ +_DURATION_MS = 4000.0 +_F_START = 1.0 +_F_END = 100.0 + + +def _chirp_and_time( + duration_ms: float = _DURATION_MS, + f_start: float = _F_START, + f_end: float = _F_END, + amplitude: float = 1.0, + dc_offset: float = 0.0, +) -> tuple[np.ndarray, np.ndarray]: + """Build a chirp stimulus and matching time axis at ``_FS_HZ``. + + Args: + duration_ms: Total duration in ms. + f_start: Starting frequency of the sweep in Hz. + f_end: Ending frequency of the sweep in Hz. + amplitude: Chirp amplitude. + dc_offset: DC offset added to the chirp. + + Returns: + Tuple of (time array in ms, chirp current array). + """ + current = patch_sim.chirp_current( + duration=duration_ms, + dc_offset=dc_offset, + amplitude=amplitude, + start_frequency=f_start, + end_frequency=f_end, + sampling_frequency=_FS_HZ, + ) + time = np.arange(current.size, dtype=float) * _DT_MS + return time, current + + +def _apply_transfer_function(current: np.ndarray, transfer: np.ndarray) -> np.ndarray: + """Apply a frequency-domain transfer function to a real signal. + + Args: + current: Real-valued input signal. + transfer: Complex transfer function evaluated at ``rfftfreq`` bins. + + Returns: + The filtered real-valued output signal, same length as ``current``. + """ + return np.fft.irfft(transfer * np.fft.rfft(current), n=current.size) + + +# --------------------------------------------------------------------------- +# Tests +# --------------------------------------------------------------------------- + + +def test_pure_resistor_flat_magnitude() -> None: + """A pure resistor (V = R·I) yields flat |Z| ≈ R and zero phase.""" + resistance = 50.0 # kΩ·cm² + time, current = _chirp_and_time() + voltage = resistance * current + + result = analyze_impedance( + time, voltage, current, 0.0, time[-1] + _DT_MS, _F_START, _F_END + ) + + assert result is not None + assert result.resonance_frequency is None + assert result.peak_impedance is None + assert result.quality_factor is None + mag = np.asarray(result.magnitude) + phase = np.asarray(result.phase) + assert mag == pytest.approx(resistance, rel=1e-6) + assert phase == pytest.approx(0.0, abs=1e-6) + + +def test_resistor_with_dc_offset_subtracted() -> None: + """DC offsets on V and I are removed before the FFT, leaving |Z| ≈ R.""" + resistance = 30.0 # kΩ·cm² + time, current = _chirp_and_time(dc_offset=5.0) + voltage = resistance * current + 12.0 + + result = analyze_impedance( + time, voltage, current, 0.0, time[-1] + _DT_MS, _F_START, _F_END + ) + + assert result is not None + assert np.asarray(result.magnitude) == pytest.approx(resistance, rel=1e-6) + assert np.asarray(result.phase) == pytest.approx(0.0, abs=1e-6) + + +def test_first_order_lowpass_no_interior_peak() -> None: + """A one-pole low-pass response has its |Z| max at the band edge, no f_R.""" + tau_s = 0.05 # 50 ms membrane time constant + time, current = _chirp_and_time() + freqs = np.fft.rfftfreq(current.size, d=1.0 / _FS_HZ) + transfer = 1.0 / (1.0 + 1j * 2.0 * np.pi * freqs * tau_s) + voltage = _apply_transfer_function(current, transfer) + + result = analyze_impedance( + time, voltage, current, 0.0, time[-1] + _DT_MS, _F_START, _F_END + ) + + assert result is not None + assert result.resonance_frequency is None + assert result.quality_factor is None + mag = np.asarray(result.magnitude) + assert int(np.argmax(mag)) == 0 + assert mag[0] > mag[-1] + + +def test_synthetic_resonance_detects_f_r() -> None: + """A second-order band-pass response yields an interior f_R and finite Q.""" + f0 = 8.0 + q0 = 2.0 + time, current = _chirp_and_time() + freqs = np.fft.rfftfreq(current.size, d=1.0 / _FS_HZ) + ratio = freqs / f0 + transfer = 1.0 / (1.0 - ratio**2 + 1j * ratio / q0) + voltage = _apply_transfer_function(current, transfer) + + result = analyze_impedance( + time, voltage, current, 0.0, time[-1] + _DT_MS, _F_START, _F_END + ) + + assert result is not None + assert result.resonance_frequency == pytest.approx(f0, abs=2.0) + assert result.peak_impedance is not None + assert result.quality_factor is not None + assert result.quality_factor > 0.3 + # The interior peak should clearly exceed the band edges. + mag = np.asarray(result.magnitude) + assert result.peak_impedance > mag[0] + assert result.peak_impedance > mag[-1] + + +def test_band_restriction() -> None: + """The analysis band is clipped to [f_start, f_end].""" + time, current = _chirp_and_time() + voltage = 10.0 * current + + result = analyze_impedance( + time, voltage, current, 0.0, time[-1] + _DT_MS, 10.0, 50.0 + ) + + assert result is not None + assert result.f_start == 10.0 + assert result.f_end == 50.0 + freqs = np.asarray(result.frequencies) + assert freqs.min() >= 10.0 + assert freqs.max() <= 50.0 + + +def test_returns_none_short_window() -> None: + """A chirp window shorter than the minimum returns None.""" + time, current = _chirp_and_time() + voltage = 10.0 * current + + result = analyze_impedance(time, voltage, current, 0.0, 10.0, _F_START, _F_END) + + assert result is None + + +def test_returns_none_invalid_band() -> None: + """A band with f_end <= f_start returns None.""" + time, current = _chirp_and_time() + voltage = 10.0 * current + + result = analyze_impedance( + time, voltage, current, 0.0, time[-1] + _DT_MS, 50.0, 50.0 + ) + + assert result is None + + +def test_returns_none_with_nans() -> None: + """Non-finite values in the voltage trace return None.""" + time, current = _chirp_and_time() + voltage = 10.0 * current + voltage[100] = np.nan + + result = analyze_impedance( + time, voltage, current, 0.0, time[-1] + _DT_MS, _F_START, _F_END + ) + + assert result is None + + +def test_returns_none_too_few_band_bins() -> None: + """A band wide enough in Hz but with too few FFT bins returns None.""" + time, current = _chirp_and_time(duration_ms=200.0) + voltage = 10.0 * current + + # 200 ms at 1 kHz → 5 Hz frequency resolution → a 1 Hz band holds < 8 bins. + result = analyze_impedance(time, voltage, current, 0.0, 200.0, 49.0, 50.0) + + assert result is None + + +def test_absolute_conversion_with_area() -> None: + """Supplying area_cm2 fills the absolute MΩ spectrum; omitting it leaves it None.""" + resistance = 40.0 # kΩ·cm² + area_cm2 = 2.0e-5 + time, current = _chirp_and_time() + voltage = resistance * current + + with_area = analyze_impedance( + time, voltage, current, 0.0, time[-1] + _DT_MS, _F_START, _F_END, area_cm2 + ) + without_area = analyze_impedance( + time, voltage, current, 0.0, time[-1] + _DT_MS, _F_START, _F_END + ) + + assert with_area is not None and without_area is not None + assert without_area.magnitude_mohm is None + assert with_area.magnitude_mohm is not None + assert len(with_area.magnitude_mohm) == len(with_area.magnitude) + expected = resistance / area_cm2 / 1000.0 + assert np.asarray(with_area.magnitude_mohm) == pytest.approx(expected, rel=1e-6) + assert with_area.area_cm2 == area_cm2 From 3843057074b3aa26f1a56a86ba3e9091f19c144d Mon Sep 17 00:00:00 2001 From: JCorson Date: Mon, 11 May 2026 15:31:54 -0500 Subject: [PATCH 2/9] feat(ui): add dedicated impedance-analysis view for the Chirp protocol MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit When a chirp (Frequency Response) current-clamp protocol is run, the analysis sidebar now shows an "Impedance Analysis" view: an f_R / Q / peak-|Z| summary plus a |Z(f)| (left axis) and ∠Z(f) in degrees (right axis) plot. The chirp impedance profile is computed in the single-sweep current-clamp branch of the sweep executor (via patch_sim.analyze_impedance), threaded through _SimResult into AnalysisState. Other protocols/clamp modes keep their existing AP-metrics / I-V views. --- patch_sim_ui/components/metrics/__init__.py | 28 ++++-- .../components/metrics/impedance_panel.py | 93 +++++++++++++++++++ patch_sim_ui/plotting/__init__.py | 2 + patch_sim_ui/plotting/impedance.py | 78 ++++++++++++++++ patch_sim_ui/state/_analysis_format.py | 63 +++++++++++++ patch_sim_ui/state/_sweep_executor.py | 20 ++++ patch_sim_ui/state/analysis.py | 20 +++- patch_sim_ui/state/simulation.py | 6 ++ tests/e2e/conftest.py | 6 ++ tests/e2e/test_current_clamp_flow.py | 40 ++++++++ 10 files changed, 347 insertions(+), 9 deletions(-) create mode 100644 patch_sim_ui/components/metrics/impedance_panel.py create mode 100644 patch_sim_ui/plotting/impedance.py diff --git a/patch_sim_ui/components/metrics/__init__.py b/patch_sim_ui/components/metrics/__init__.py index 9ffb4899..17977158 100644 --- a/patch_sim_ui/components/metrics/__init__.py +++ b/patch_sim_ui/components/metrics/__init__.py @@ -1,8 +1,8 @@ """Right-hand analysis sidebar — split into per-analysis panels. -A collapsible panel that automatically shows the relevant analysis view -based on the active clamp mode: AP metrics for current clamp, I-V curve -for voltage clamp. +A collapsible panel that automatically shows the relevant analysis view: +the impedance profile for a chirp current-clamp protocol, AP metrics for any +other current-clamp protocol, and the I-V curve in voltage clamp mode. The shell pieces (header, passive-membrane row, expand/collapse strips, ``analysis_sidebar`` entry point) live in this ``__init__.py``; the @@ -12,6 +12,7 @@ - :mod:`patch_sim_ui.components.metrics.burst_panel` — burst summary/table - :mod:`patch_sim_ui.components.metrics.calcium_panel` — calcium summary/table - :mod:`patch_sim_ui.components.metrics.fi_gv_panel` — F-I, g-V, τ-V, I-V tab +- :mod:`patch_sim_ui.components.metrics.impedance_panel` — chirp impedance profile - :mod:`patch_sim_ui.components.metrics.sfa_hyperpol_panel` — SFA/sag plots - :mod:`patch_sim_ui.components.metrics._row` — shared dict-to-row helper @@ -29,6 +30,7 @@ _iv_curve_tab, _tau_v_plot, ) +from patch_sim_ui.components.metrics.impedance_panel import impedance_tab from patch_sim_ui.state import SimulationState from patch_sim_ui.state.analysis import AnalysisState from patch_sim_ui.state.protocol import ProtocolState @@ -111,9 +113,9 @@ def _membrane_test_section() -> rx.Component: def _expanded_panel() -> rx.Component: """Render the full expanded analysis sidebar. - The analysis view shown depends on the active clamp mode: AP metrics - are shown in current clamp mode, and the I-V curve is shown in voltage - clamp mode. + The analysis view shown depends on the active clamp mode and protocol: + the impedance profile for a chirp current-clamp protocol, AP metrics for + any other current-clamp protocol, and the I-V curve in voltage clamp mode. Returns: A fixed-width flex column with a header and the mode-appropriate @@ -124,7 +126,12 @@ def _expanded_panel() -> rx.Component: rx.icon("chart-line", size=14), rx.cond( ProtocolState.clamp_mode == CURRENT_CLAMP, - rx.text("AP Analysis", size="4", weight="bold"), + rx.cond( + # "Chirp" is the protocol-type name from CURRENT_PROTOCOLS. + ProtocolState.protocol_type == "Chirp", + rx.text("Impedance Analysis", size="4", weight="bold"), + rx.text("AP Analysis", size="4", weight="bold"), + ), rx.text("I-V Analysis", size="4", weight="bold"), ), rx.spacer(), @@ -145,7 +152,12 @@ def _expanded_panel() -> rx.Component: rx.box( rx.cond( ProtocolState.clamp_mode == CURRENT_CLAMP, - _ap_metrics_tab(), + rx.cond( + # "Chirp" is the protocol-type name from CURRENT_PROTOCOLS. + ProtocolState.protocol_type == "Chirp", + impedance_tab(), + _ap_metrics_tab(), + ), _iv_curve_tab(), ), display="flex", diff --git a/patch_sim_ui/components/metrics/impedance_panel.py b/patch_sim_ui/components/metrics/impedance_panel.py new file mode 100644 index 00000000..d564a834 --- /dev/null +++ b/patch_sim_ui/components/metrics/impedance_panel.py @@ -0,0 +1,93 @@ +"""Impedance-analysis panel: resonance summary and |Z| / phase plot. + +Shown in the analysis sidebar in place of the AP-metrics panel when the +active current-clamp protocol is a chirp (Frequency Response). +""" + +import reflex as rx + +from patch_sim_ui.state.analysis import AnalysisState + + +def _impedance_summary() -> rx.Component: + """Render the resonance-summary section (resonance frequency, Q, peak |Z|). + + Returns: + A compact grid of labeled impedance metrics drawn from + ``AnalysisState.impedance_data``. + """ + d = AnalysisState.impedance_data + return rx.box( + rx.grid( + rx.text("Resonance f_R", size="1", color="gray"), + rx.text(d["resonance_frequency"].to(str) + " Hz", size="1", align="left"), + rx.text("Quality Q", size="1", color="gray"), + rx.text(d["quality_factor"].to(str), size="1", align="left"), + rx.text("Peak |Z|", size="1", color="gray"), + rx.text( + d["peak_impedance"].to(str) + " " + d["units"].to(str), + size="1", + align="left", + ), + columns="2", + spacing="2", + width="100%", + ), + padding="3", + border_bottom="1px solid var(--gray-4)", + width="100%", + ) + + +def _impedance_plot() -> rx.Component: + """Render the embedded impedance-profile plot (|Z| and phase vs frequency). + + Returns: + A flex container holding the impedance Plotly figure. + """ + return rx.flex( + rx.plotly( + data=AnalysisState.impedance_figure, + width="100%", + ), + direction="column", + width="100%", + flex_shrink="0", + border_top="1px solid var(--gray-4)", + padding="1", + ) + + +def impedance_tab() -> rx.Component: + """Render the impedance-analysis tab for the chirp current-clamp pane. + + Shows the resonance summary and the |Z| / phase plot when impedance data + is available, and a placeholder otherwise. + + Returns: + The impedance-analysis content as a scrollable flex column. + """ + return rx.cond( + AnalysisState.has_impedance_data, + rx.scroll_area( + rx.flex( + _impedance_summary(), + _impedance_plot(), + direction="column", + width="100%", + ), + height="100%", + width="100%", + ), + rx.flex( + rx.text( + "Run a Frequency Response (chirp) current clamp simulation to " + "see the impedance profile.", + size="1", + color="gray", + text_align="center", + ), + padding="4", + justify="center", + ), + ) diff --git a/patch_sim_ui/plotting/__init__.py b/patch_sim_ui/plotting/__init__.py index 51ba2bfd..b150ab14 100644 --- a/patch_sim_ui/plotting/__init__.py +++ b/patch_sim_ui/plotting/__init__.py @@ -8,6 +8,7 @@ from patch_sim_ui.plotting.fi_curve import build_fi_figure 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.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 @@ -24,6 +25,7 @@ "build_fi_figure", "build_gv_figure", "build_hyperpolarization_figure", + "build_impedance_figure", "build_iv_figure", "build_phase_plane_figure", "build_sfa_figure", diff --git a/patch_sim_ui/plotting/impedance.py b/patch_sim_ui/plotting/impedance.py new file mode 100644 index 00000000..8cf6087b --- /dev/null +++ b/patch_sim_ui/plotting/impedance.py @@ -0,0 +1,78 @@ +"""Impedance-profile figure (|Z| magnitude and phase vs frequency).""" + +import plotly.graph_objects as go + +from patch_sim_ui.plotting._layout import ANALYSIS_FIGURE_LAYOUT + + +def build_impedance_figure(imp_data: dict) -> go.Figure: + """Build a Plotly impedance-profile figure from serialised chirp results. + + Renders the impedance magnitude (left y-axis, teal) and phase (right + y-axis, orange, dashed) plotted against frequency (Hz). When a resonance + frequency was detected, a dotted vertical marker is drawn at that + frequency. + + Args: + imp_data: Dict with keys ``frequencies``, ``magnitude``, ``phase`` + (lists aligned by index), ``units`` (magnitude unit string), and + ``resonance_frequency`` (pre-formatted string, ``"—"`` when none). + + Returns: + A Plotly :class:`go.Figure` ready to be serialised and sent to the UI. + """ + freqs = imp_data["frequencies"] + magnitude = imp_data["magnitude"] + phase = imp_data["phase"] + units = imp_data.get("units", "kΩ·cm²") + + fig = go.Figure() + fig.add_trace( + go.Scatter( + x=freqs, + y=magnitude, + mode="lines+markers", + name="|Z|", + line=dict(color="#16a085"), + marker=dict(size=4), + yaxis="y1", + ) + ) + fig.add_trace( + go.Scatter( + x=freqs, + y=phase, + mode="lines", + name="Phase", + line=dict(color="#e67e22", dash="dash"), + yaxis="y2", + ) + ) + fig.update_layout( + **ANALYSIS_FIGURE_LAYOUT, + xaxis_title="Frequency (Hz)", + yaxis=dict(title=f"|Z| ({units})", color="#16a085"), + yaxis2=dict( + title="Phase (°)", + overlaying="y", + side="right", + color="#e67e22", + ), + showlegend=True, + legend=dict( + orientation="h", + yanchor="bottom", + y=1.02, + xanchor="right", + x=1, + font=dict(size=10), + ), + ) + + f_r = imp_data.get("resonance_frequency") + if f_r is not None: + try: + fig.add_vline(x=float(f_r), line_dash="dot", line_color="gray") + except (TypeError, ValueError): + pass + return fig diff --git a/patch_sim_ui/state/_analysis_format.py b/patch_sim_ui/state/_analysis_format.py index ce5b7a79..1dde6dd4 100644 --- a/patch_sim_ui/state/_analysis_format.py +++ b/patch_sim_ui/state/_analysis_format.py @@ -464,6 +464,69 @@ def _build_phase_plane_data(sweeps: "list[Sweep]") -> dict[str, Any]: return {"sweeps": eligible} if eligible else {} +def _compute_impedance_data( + sweep: "Sweep", + start_frequency: float, + end_frequency: float, + area_cm2: "float | None", +) -> dict[str, Any]: + """Serialise a chirp current-clamp sweep into impedance data for AnalysisState. + + Computes the FFT-based membrane impedance over the chirp's swept band via + :func:`patch_sim.analyze_impedance`. The whole recording is treated as the + chirp window (the Frequency Response preset has no pre/post-stimulus + padding). Absolute MΩ values are used when the neuron declares a surface + area, otherwise the per-area kΩ·cm² spectrum is reported. + + Args: + sweep: The single current-clamp sweep produced by a chirp protocol; + its ``stimulus`` field holds the injected current (µA/cm²). + start_frequency: Chirp sweep start frequency in Hz. + end_frequency: Chirp sweep end frequency in Hz. + area_cm2: Membrane surface area in cm², or ``None`` when undeclared. + + Returns: + A dict with keys ``frequencies``, ``magnitude``, ``phase`` (lists), + ``resonance_frequency``, ``quality_factor``, ``peak_impedance`` + (pre-formatted strings), and ``units``. Returns an empty dict when the + impedance analysis is not applicable (see :func:`patch_sim.analyze_impedance`). + """ + time_arr = np.asarray(sweep.time, dtype=float) + if len(time_arr) < 2: + return {} + voltage_arr = np.asarray(sweep.voltage, dtype=float) + current_arr = np.asarray(sweep.stimulus, dtype=float) + profile = patch_sim.analyze_impedance( + time_arr, + voltage_arr, + current_arr, + float(time_arr[0]), + float(time_arr[-1]), + start_frequency, + end_frequency, + area_cm2=area_cm2, + ) + if profile is None: + return {} + if profile.magnitude_mohm is not None: + magnitude = profile.magnitude_mohm + peak = profile.peak_impedance_mohm + units = "MΩ" + else: + magnitude = profile.magnitude + peak = profile.peak_impedance + units = "kΩ·cm²" + return { + "frequencies": profile.frequencies, + "magnitude": magnitude, + "phase": profile.phase, + "resonance_frequency": _fmt_optional(profile.resonance_frequency, ".2f"), + "quality_factor": _fmt_optional(profile.quality_factor, ".2f"), + "peak_impedance": _fmt_optional(peak, ".3g"), + "units": units, + } + + def _compute_cc_multi_sweep_analysis( sweeps: "list[Sweep]", min_stimulus: float, diff --git a/patch_sim_ui/state/_sweep_executor.py b/patch_sim_ui/state/_sweep_executor.py index f0dae094..bf2db74b 100644 --- a/patch_sim_ui/state/_sweep_executor.py +++ b/patch_sim_ui/state/_sweep_executor.py @@ -25,6 +25,7 @@ _compute_ca_transient_data, _compute_cc_multi_sweep_analysis, _compute_gv_data, + _compute_impedance_data, _compute_iv_data, _compute_multi_sweep_burst_data, _compute_multi_sweep_ca_transient_data, @@ -62,6 +63,7 @@ class _SimResult: sfa_data: dict[str, Any] = dataclasses.field(default_factory=dict) hyperpolarization_data: dict[str, Any] = dataclasses.field(default_factory=dict) phase_plane_data: dict[str, Any] = dataclasses.field(default_factory=dict) + impedance_data: dict[str, Any] = dataclasses.field(default_factory=dict) def _compute_simulation( @@ -75,6 +77,9 @@ def _compute_simulation( stimulus_step: float, pre_stimulus_duration: float, stimulus_duration: float, + protocol_type: str = "", + start_frequency: float = 0.0, + end_frequency: float = 0.0, ) -> _SimResult: """Run the simulation synchronously and compute all analysis. @@ -95,6 +100,12 @@ def _compute_simulation( stimulus_step: Stimulus step size for analysis range. pre_stimulus_duration: Pre-stimulus duration (ms) for analysis windows. stimulus_duration: Stimulus duration (ms) for analysis windows. + protocol_type: Protocol-type name (e.g. ``"Step"``, ``"Chirp"``). Used + to decide which protocol-specific analyses to run. + start_frequency: Chirp sweep start frequency (Hz); only used when + ``protocol_type == "Chirp"``. + end_frequency: Chirp sweep end frequency (Hz); only used when + ``protocol_type == "Chirp"``. Returns: A :class:`_SimResult` containing sweeps, figure token, and all @@ -239,6 +250,14 @@ def _compute_simulation( burst_metrics, burst_summary = _compute_burst_data( ap_result, np.asarray(result["time"]) ) + # "Chirp" is the protocol-type name from constants.CURRENT_PROTOCOLS. + impedance_data = ( + _compute_impedance_data( + sweep, start_frequency, end_frequency, neuron.area_cm2 + ) + if protocol_type == "Chirp" + else {} + ) return _SimResult( sweeps=[sweep], sim_token=sim_token, @@ -273,6 +292,7 @@ def _compute_simulation( else {} ), phase_plane_data=_build_phase_plane_data([sweep]), + impedance_data=impedance_data, ) # Single VC sweep — IV/GV require multi-sweep, but calcium transients can diff --git a/patch_sim_ui/state/analysis.py b/patch_sim_ui/state/analysis.py index 17f12990..325b6cc5 100644 --- a/patch_sim_ui/state/analysis.py +++ b/patch_sim_ui/state/analysis.py @@ -9,6 +9,7 @@ build_fi_figure, build_gv_figure, build_hyperpolarization_figure, + build_impedance_figure, build_iv_figure, build_phase_plane_figure, build_sfa_figure, @@ -17,7 +18,7 @@ class AnalysisState(rx.State): - """State for AP, F-I, I-V, g-V, SFA, and phase-plane analysis results.""" + """State for AP, F-I, I-V, g-V, SFA, phase-plane, and impedance results.""" ap_metrics: list[dict[str, Any]] = [] # Per-spike metrics (serialized) ap_summary: dict[str, Any] = {} # Aggregate summary statistics @@ -37,6 +38,7 @@ class AnalysisState(rx.State): sfa_data: dict[str, Any] = {} # Serialized SFAAnalysisResult for the UI hyperpolarization_data: dict[str, Any] = {} # Serialized sag/rebound analysis phase_plane_data: dict[str, Any] = {} # Serialized V vs dV/dt sweep data + impedance_data: dict[str, Any] = {} # Serialized chirp impedance profile # Membrane test results — persisted across protocol/simulation changes. # Only invalidated when neuron parameters change (neuron_fingerprint mismatch). @@ -69,6 +71,7 @@ def clear_results(self) -> None: self.sfa_data = {} self.hyperpolarization_data = {} self.phase_plane_data = {} + self.impedance_data = {} @rx.var def has_membrane_test(self) -> bool: @@ -215,3 +218,18 @@ def phase_plane_figure(self) -> go.Figure: if not self.phase_plane_data: return go.Figure() return build_phase_plane_figure(self.phase_plane_data) + + @rx.var + def has_impedance_data(self) -> bool: + """Return True when chirp impedance-profile results are available.""" + return len(self.impedance_data) > 0 + + @rx.var + def impedance_figure(self) -> go.Figure: + """Return a Plotly impedance-profile figure (|Z| and phase vs frequency). + + Returns an empty figure when no impedance data is available. + """ + if not self.impedance_data: + return go.Figure() + return build_impedance_figure(self.impedance_data) diff --git a/patch_sim_ui/state/simulation.py b/patch_sim_ui/state/simulation.py index 9f1b6ad6..8cb08c69 100644 --- a/patch_sim_ui/state/simulation.py +++ b/patch_sim_ui/state/simulation.py @@ -474,6 +474,7 @@ def _do_apply_simulation( analysis_st.sfa_data = result.sfa_data analysis_st.hyperpolarization_data = result.hyperpolarization_data analysis_st.phase_plane_data = result.phase_plane_data + analysis_st.impedance_data = result.impedance_data # ------------------------------------------------------------------ # # Continuous simulation mode # @@ -699,6 +700,8 @@ async def run_simulation(self) -> AsyncGenerator[Any, None]: stimulus_step = proto_st.stimulus_step pre_stimulus_duration = proto_st.pre_stimulus_duration stimulus_duration = proto_st.stimulus_duration + start_frequency = proto_st.start_frequency + end_frequency = proto_st.end_frequency logger.info("Simulation started: mode=%s, protocol=%s", mode, ptype) @@ -719,6 +722,9 @@ async def run_simulation(self) -> AsyncGenerator[Any, None]: stimulus_step, pre_stimulus_duration, stimulus_duration, + ptype, + start_frequency, + end_frequency, ) except ValueError as exc: logger.exception("Simulation error: %s", exc) diff --git a/tests/e2e/conftest.py b/tests/e2e/conftest.py index 64b60abf..769d7328 100644 --- a/tests/e2e/conftest.py +++ b/tests/e2e/conftest.py @@ -217,6 +217,9 @@ async def run_flow( stimulus_step=tree.protocol.stimulus_step, pre_stimulus_duration=tree.protocol.pre_stimulus_duration, stimulus_duration=tree.protocol.stimulus_duration, + protocol_type=tree.protocol.protocol_type, + start_frequency=tree.protocol.start_frequency, + end_frequency=tree.protocol.end_frequency, ) tree.sim._do_apply_simulation(result, tree.analysis) return result @@ -249,6 +252,9 @@ def simulate_and_apply(tree: StateTree) -> _SimResult: stimulus_step=tree.protocol.stimulus_step, pre_stimulus_duration=tree.protocol.pre_stimulus_duration, stimulus_duration=tree.protocol.stimulus_duration, + protocol_type=tree.protocol.protocol_type, + start_frequency=tree.protocol.start_frequency, + end_frequency=tree.protocol.end_frequency, ) tree.sim._do_apply_simulation(result, tree.analysis) return result diff --git a/tests/e2e/test_current_clamp_flow.py b/tests/e2e/test_current_clamp_flow.py index 5a80a7b0..97bdc786 100644 --- a/tests/e2e/test_current_clamp_flow.py +++ b/tests/e2e/test_current_clamp_flow.py @@ -12,6 +12,7 @@ CORTICAL_PYRAMIDAL, DOPAMINERGIC, FAST_SPIKING_INTERNEURON, + FREQUENCY_RESPONSE, PURKINJE, REPETITIVE_FIRING, SQUID_GIANT_AXON, @@ -111,6 +112,45 @@ async def test_repetitive_firing_populates_phase_plane(state_tree: StateTree) -> assert state_tree.analysis.phase_plane_data != {} +async def test_frequency_response_preset_populates_impedance( + state_tree: StateTree, +) -> None: + """The Frequency Response (chirp) preset populates impedance_data after a run.""" + await run_flow( + state_tree, + neuron_preset=SQUID_GIANT_AXON, + protocol_preset=FREQUENCY_RESPONSE, + ) + + impedance = state_tree.analysis.impedance_data + assert impedance != {} + expected_keys = { + "frequencies", + "magnitude", + "phase", + "resonance_frequency", + "quality_factor", + "peak_impedance", + "units", + } + assert expected_keys.issubset(impedance.keys()) + assert len(impedance["frequencies"]) == len(impedance["magnitude"]) + assert len(impedance["frequencies"]) > 0 + + +async def test_action_potential_preset_leaves_impedance_empty( + state_tree: StateTree, +) -> None: + """A non-chirp protocol does not populate impedance_data.""" + await run_flow( + state_tree, + neuron_preset=SQUID_GIANT_AXON, + protocol_preset=ACTION_POTENTIAL, + ) + + assert state_tree.analysis.impedance_data == {} + + async def test_apply_simulation_clears_previous_analysis( state_tree: StateTree, ) -> None: From 4d55f586df3f46953d314a3d410ab41ca22d2030 Mon Sep 17 00:00:00 2001 From: JCorson Date: Mon, 11 May 2026 15:55:01 -0500 Subject: [PATCH 3/9] fix(analysis): only report chirp impedance for subthreshold responses MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Addresses code review of #365: - analyze_impedance now returns None when the windowed response is suprathreshold — membrane impedance is a linear small-signal quantity and spike transients would otherwise dominate the spectrum and the f_R/Q metrics. The UI shows the placeholder (with a hint to reduce the chirp amplitude) for cells that spike. - Trim contiguous band-edge FFT bins where the chirp has negligible spectral power, so dividing by a near-zero stimulus can't inflate the band-edge estimate or steer the resonance-peak detection. - Restrict the analysis to the actual chirp window (pre/stimulus durations) rather than the whole recording, so protocol padding is honored. - Use a small low-frequency bin average as the resonance reference (robust to first-bin spectral leakage). - Document the FREQUENCY_RESPONSE preset's suprathreshold default amplitude. - Tests: subthreshold synthetic signals; a suprathreshold→None unit test; a CA1 Ih-resonance integration test; a suprathreshold→None integration test; e2e coverage of the populated, empty, and non-chirp paths. --- patch_sim/analysis/impedance.py | 69 +++++++++++++---- patch_sim/presets/protocols.py | 5 ++ .../components/metrics/impedance_panel.py | 3 +- patch_sim_ui/state/_analysis_format.py | 24 ++++-- patch_sim_ui/state/_sweep_executor.py | 7 +- tests/e2e/test_current_clamp_flow.py | 43 ++++++++--- .../integration/test_impedance_simulation.py | 43 ++++++++--- tests/unit/test_impedance.py | 75 +++++++++++++------ 8 files changed, 203 insertions(+), 66 deletions(-) diff --git a/patch_sim/analysis/impedance.py b/patch_sim/analysis/impedance.py index 5cb25415..009caf74 100644 --- a/patch_sim/analysis/impedance.py +++ b/patch_sim/analysis/impedance.py @@ -18,7 +18,10 @@ import numpy as np -from patch_sim.analysis.passive_properties import density_to_absolute_r_in +from patch_sim.analysis.passive_properties import ( + density_to_absolute_r_in, + is_subthreshold, +) logger = logging.getLogger(__name__) @@ -32,10 +35,16 @@ _MIN_BAND_BINS: int = 8 #: Fractional margin by which the ``|Z|`` peak must exceed the low-frequency -#: band edge to count as a genuine resonance (rather than measurement ripple -#: on an otherwise monotone passive response). +#: reference to count as a genuine resonance (rather than measurement ripple on +#: an otherwise monotone passive response). _RESONANCE_REL_MARGIN: float = 0.01 +#: A linear chirp's spectral power tapers toward the band edges. Contiguous +#: leading/trailing FFT bins whose stimulus magnitude falls below this fraction +#: of the in-band maximum are trimmed before forming ``Z(f)``, so that dividing +#: by a near-zero stimulus does not inflate the band-edge impedance estimate. +_MIN_STIM_FRAC: float = 0.05 + @dataclasses.dataclass class ImpedanceProfile: @@ -49,7 +58,7 @@ class ImpedanceProfile: ``resonance_frequency``, ``peak_impedance`` and ``quality_factor`` are all ``None`` when ``|Z(f)|`` has no interior peak that rises meaningfully above - the low-frequency band edge — the hallmark of a passive low-pass cell with + the low-frequency reference — the hallmark of a passive low-pass cell with no resonance. ``quality_factor`` may additionally be ``None`` when an interior peak exists but its half-power width cannot be bracketed inside the analysis band. @@ -177,10 +186,15 @@ def analyze_impedance( Extracts the chirp window from the traces, removes the DC component from both signals, takes the real FFT of each (``numpy.fft.rfft``), forms - ``Z(f) = V̂(f) / Î(f)`` over the swept band ``[f_start, f_end]``, and - reports the magnitude and phase spectra plus the resonance frequency - (argmax of ``|Z|``, only when it is an interior maximum that rises above - the low-frequency band edge) and quality factor. + ``Z(f) = V̂(f) / Î(f)`` over the swept band ``[f_start, f_end]`` (trimming + band-edge bins where the chirp has negligible spectral power), and reports + the magnitude and phase spectra plus the resonance frequency (argmax of + ``|Z|``, only when it is an interior maximum that rises above the + low-frequency reference) and quality factor. + + Membrane impedance is a linear, small-signal quantity: the analysis bails + out (returns ``None``) when the windowed response is suprathreshold, since + spike transients dominate the spectrum and the result would be meaningless. Args: time: Time axis in ms (uniformly sampled). @@ -198,9 +212,9 @@ def analyze_impedance( An :class:`ImpedanceProfile`, or ``None`` when the inputs are inconsistent (mismatched lengths, non-finite values), the band is invalid (``f_end <= f_start`` or ``f_start < 0``), the chirp window is - shorter than :data:`_MIN_WINDOW_MS`, fewer than :data:`_MIN_BAND_BINS` - FFT bins fall inside the band, or the stimulus has no spectral power in - the band. + shorter than :data:`_MIN_WINDOW_MS`, the windowed response is + suprathreshold, the stimulus has no spectral power in the band, or + fewer than :data:`_MIN_BAND_BINS` FFT bins remain after trimming. """ t = np.asarray(time, dtype=float) v = np.asarray(voltage, dtype=float) @@ -230,6 +244,13 @@ def analyze_impedance( return None t_win, v_win, i_win = t[mask], v[mask], i[mask] + if not is_subthreshold(t_win, v_win): + logger.warning( + "analyze_impedance: response is suprathreshold; impedance is only " + "defined in the subthreshold regime; skipping." + ) + return None + dt_ms = float(np.mean(np.diff(t_win))) if dt_ms <= 0.0: logger.warning("analyze_impedance: non-positive sampling interval; skipping.") @@ -263,17 +284,36 @@ def analyze_impedance( v_band = v_fft[band] i_band = i_fft[band] - if float(np.max(np.abs(i_band))) <= 0.0: + i_mag = np.abs(i_band) + i_max = float(np.max(i_mag)) + if i_max <= 0.0: logger.warning("analyze_impedance: stimulus has no power in band; skipping.") return None + # Trim contiguous band-edge bins where the chirp has negligible power so + # that dividing by a near-zero stimulus does not inflate the estimate. + keep = np.flatnonzero(i_mag >= _MIN_STIM_FRAC * i_max) + lo, hi = int(keep[0]), int(keep[-1]) + 1 + if hi - lo < _MIN_BAND_BINS: + logger.warning( + "analyze_impedance: only %d FFT bins with usable stimulus power; " + "need %d; skipping.", + hi - lo, + _MIN_BAND_BINS, + ) + return None + fb = fb[lo:hi] + v_band = v_band[lo:hi] + i_band = i_band[lo:hi] z = v_band / i_band mag = np.abs(z) phase_deg = np.degrees(np.angle(z)) peak_idx = int(np.argmax(mag)) - low_edge = float(mag[0]) - is_resonant = 0 < peak_idx < mag.size - 1 and float(mag[peak_idx]) > low_edge * ( + # Reference the low-frequency end with a small average to be robust to + # spectral leakage in the first bin. + low_ref = float(np.mean(mag[: min(3, mag.size)])) + is_resonant = 0 < peak_idx < mag.size - 1 and float(mag[peak_idx]) > low_ref * ( 1.0 + _RESONANCE_REL_MARGIN ) if is_resonant: @@ -286,6 +326,7 @@ def analyze_impedance( quality_factor = None if area_cm2 is not None and area_cm2 > 0.0: + # Per the density_to_absolute_r_in derivation: kΩ·cm² / area / 1000 = MΩ. magnitude_mohm: list[float] | None = [float(m) / area_cm2 / 1000.0 for m in mag] else: magnitude_mohm = None diff --git a/patch_sim/presets/protocols.py b/patch_sim/presets/protocols.py index a1c4d948..2da38237 100644 --- a/patch_sim/presets/protocols.py +++ b/patch_sim/presets/protocols.py @@ -97,6 +97,11 @@ "max_stimulus": 60.0, "stimulus_step": 10.0, }, + # NOTE: this DC offset / amplitude drives most presets above spike + # threshold, so the FFT-based impedance analysis (valid only in the linear + # subthreshold regime) reports nothing for those cells. To measure a clean + # impedance / subthreshold-resonance profile, reduce the amplitude (and zero + # the DC offset) so the chirp response stays below threshold. FREQUENCY_RESPONSE: { "clamp_mode": CURRENT_CLAMP, "protocol_type": "Chirp", diff --git a/patch_sim_ui/components/metrics/impedance_panel.py b/patch_sim_ui/components/metrics/impedance_panel.py index d564a834..736f4c90 100644 --- a/patch_sim_ui/components/metrics/impedance_panel.py +++ b/patch_sim_ui/components/metrics/impedance_panel.py @@ -82,7 +82,8 @@ def impedance_tab() -> rx.Component: rx.flex( rx.text( "Run a Frequency Response (chirp) current clamp simulation to " - "see the impedance profile.", + "see the impedance profile. The chirp must keep the cell " + "subthreshold — reduce the amplitude if the response spikes.", size="1", color="gray", text_align="center", diff --git a/patch_sim_ui/state/_analysis_format.py b/patch_sim_ui/state/_analysis_format.py index 1dde6dd4..bfee1a83 100644 --- a/patch_sim_ui/state/_analysis_format.py +++ b/patch_sim_ui/state/_analysis_format.py @@ -466,6 +466,8 @@ def _build_phase_plane_data(sweeps: "list[Sweep]") -> dict[str, Any]: def _compute_impedance_data( sweep: "Sweep", + pre_stimulus_duration: float, + stimulus_duration: float, start_frequency: float, end_frequency: float, area_cm2: "float | None", @@ -473,14 +475,18 @@ def _compute_impedance_data( """Serialise a chirp current-clamp sweep into impedance data for AnalysisState. Computes the FFT-based membrane impedance over the chirp's swept band via - :func:`patch_sim.analyze_impedance`. The whole recording is treated as the - chirp window (the Frequency Response preset has no pre/post-stimulus - padding). Absolute MΩ values are used when the neuron declares a surface - area, otherwise the per-area kΩ·cm² spectrum is reported. + :func:`patch_sim.analyze_impedance`, restricting the analysis to the chirp + window ``[pre_stimulus_duration, pre_stimulus_duration + stimulus_duration]`` + (or to the whole recording when ``stimulus_duration`` is zero, matching the + chirp builder). Absolute MΩ values are used when the neuron declares a + surface area, otherwise the per-area kΩ·cm² spectrum is reported. Args: sweep: The single current-clamp sweep produced by a chirp protocol; its ``stimulus`` field holds the injected current (µA/cm²). + pre_stimulus_duration: Pre-stimulus padding before the chirp (ms). + stimulus_duration: Duration of the chirp stimulus (ms); zero means the + chirp fills the recording from ``pre_stimulus_duration`` onward. start_frequency: Chirp sweep start frequency in Hz. end_frequency: Chirp sweep end frequency in Hz. area_cm2: Membrane surface area in cm², or ``None`` when undeclared. @@ -496,12 +502,18 @@ def _compute_impedance_data( return {} voltage_arr = np.asarray(sweep.voltage, dtype=float) current_arr = np.asarray(sweep.stimulus, dtype=float) + stim_start_ms = float(pre_stimulus_duration) + stim_end_ms = ( + float(pre_stimulus_duration + stimulus_duration) + if stimulus_duration > 0.0 + else float(time_arr[-1]) + ) profile = patch_sim.analyze_impedance( time_arr, voltage_arr, current_arr, - float(time_arr[0]), - float(time_arr[-1]), + stim_start_ms, + stim_end_ms, start_frequency, end_frequency, area_cm2=area_cm2, diff --git a/patch_sim_ui/state/_sweep_executor.py b/patch_sim_ui/state/_sweep_executor.py index bf2db74b..825f9c30 100644 --- a/patch_sim_ui/state/_sweep_executor.py +++ b/patch_sim_ui/state/_sweep_executor.py @@ -253,7 +253,12 @@ def _compute_simulation( # "Chirp" is the protocol-type name from constants.CURRENT_PROTOCOLS. impedance_data = ( _compute_impedance_data( - sweep, start_frequency, end_frequency, neuron.area_cm2 + sweep, + pre_stimulus_duration, + stimulus_duration, + start_frequency, + end_frequency, + neuron.area_cm2, ) if protocol_type == "Chirp" else {} diff --git a/tests/e2e/test_current_clamp_flow.py b/tests/e2e/test_current_clamp_flow.py index 97bdc786..ed3a713b 100644 --- a/tests/e2e/test_current_clamp_flow.py +++ b/tests/e2e/test_current_clamp_flow.py @@ -18,7 +18,12 @@ SQUID_GIANT_AXON, THALAMIC_RELAY, ) -from tests.e2e.conftest import StateTree, run_flow +from tests.e2e.conftest import ( + StateTree, + patch_get_state, + run_flow, + simulate_and_apply, +) @pytest.mark.parametrize( @@ -112,15 +117,20 @@ async def test_repetitive_firing_populates_phase_plane(state_tree: StateTree) -> assert state_tree.analysis.phase_plane_data != {} -async def test_frequency_response_preset_populates_impedance( - state_tree: StateTree, -) -> None: - """The Frequency Response (chirp) preset populates impedance_data after a run.""" - await run_flow( - state_tree, - neuron_preset=SQUID_GIANT_AXON, - protocol_preset=FREQUENCY_RESPONSE, - ) +async def test_subthreshold_chirp_populates_impedance(state_tree: StateTree) -> None: + """A subthreshold chirp through the full pipeline populates impedance_data. + + Loads the Frequency Response (chirp) preset on the squid axon, then dials + the chirp down to a subthreshold amplitude so ``analyze_impedance`` has a + valid linear response to work with. + """ + async with patch_get_state(state_tree): + [_ async for _ in state_tree.neuron.load_neuron_preset(SQUID_GIANT_AXON)] + [_ async for _ in state_tree.protocol.load_protocol_preset(FREQUENCY_RESPONSE)] + state_tree.protocol.dc_offset = 0.0 + state_tree.protocol.amplitude = 1.0 + + simulate_and_apply(state_tree) impedance = state_tree.analysis.impedance_data assert impedance != {} @@ -138,6 +148,19 @@ async def test_frequency_response_preset_populates_impedance( assert len(impedance["frequencies"]) > 0 +async def test_suprathreshold_chirp_leaves_impedance_empty( + state_tree: StateTree, +) -> None: + """The default (large) Frequency Response chirp drives spiking → no impedance.""" + await run_flow( + state_tree, + neuron_preset=SQUID_GIANT_AXON, + protocol_preset=FREQUENCY_RESPONSE, + ) + + assert state_tree.analysis.impedance_data == {} + + async def test_action_potential_preset_leaves_impedance_empty( state_tree: StateTree, ) -> None: diff --git a/tests/integration/test_impedance_simulation.py b/tests/integration/test_impedance_simulation.py index 95db89f0..5ef61b2a 100644 --- a/tests/integration/test_impedance_simulation.py +++ b/tests/integration/test_impedance_simulation.py @@ -1,8 +1,8 @@ """Integration tests for patch_sim.analysis.impedance against real simulations. Runs a chirp current-clamp protocol through the full simulation pipeline and -verifies that the resulting impedance profile is structurally sound and -physiologically scaled. +verifies that the resulting impedance profile is structurally sound, +physiologically scaled, and only produced for subthreshold responses. Unit tests with synthetic signals live in tests/unit/test_impedance.py. """ @@ -10,7 +10,7 @@ import patch_sim from patch_sim.clamp_simulations import SIM_SAMPLING_FREQ -from patch_sim.presets import make_thalamic_relay +from patch_sim.presets import make_ca1_pyramidal, make_thalamic_relay _DURATION_MS = 500.0 _F_START = 1.0 @@ -67,16 +67,35 @@ def test_chirp_pipeline_produces_sane_impedance_profile(hh_model) -> None: assert np.all(np.abs(phase) <= 180.0) -def test_thalamic_relay_has_higher_low_frequency_impedance_than_squid( - hh_model, -) -> None: - """A thalamic relay neuron shows much larger low-frequency |Z| than the squid axon. +def test_ca1_chirp_shows_resonance_and_high_impedance(hh_model) -> None: + """A hippocampal CA1 cell (with Ih) shows a low-frequency |Z| resonance. - Reflects the far higher input resistance of a small thalamic relay cell - compared with the large, low-resistance squid giant axon. + The CA1 pyramidal preset has an Ih (HCN) conductance, so a subthreshold + chirp produces an interior peak in |Z(f)| in the few-Hz range, and its + low-frequency impedance is much larger than the low-resistance squid axon. """ - relay_profile, _ = _run_chirp_profile(make_thalamic_relay(), amplitude=1.0) + ca1_profile, ca1_voltage = _run_chirp_profile(make_ca1_pyramidal(), amplitude=1.0) squid_profile, _ = _run_chirp_profile(hh_model, amplitude=1.0) - assert relay_profile is not None and squid_profile is not None - assert relay_profile.magnitude[0] > 5.0 * squid_profile.magnitude[0] + assert ca1_voltage.max() < -20.0 # stayed subthreshold + assert ca1_profile is not None and squid_profile is not None + assert ca1_profile.resonance_frequency is not None + assert 1.0 < ca1_profile.resonance_frequency < 30.0 + assert ca1_profile.peak_impedance is not None + mag = np.asarray(ca1_profile.magnitude) + assert ca1_profile.peak_impedance > mag[0] + assert ca1_profile.peak_impedance > mag[-1] + assert ca1_profile.magnitude[0] > 5.0 * squid_profile.magnitude[0] + + +def test_suprathreshold_chirp_returns_none() -> None: + """A chirp that drives spiking yields no impedance profile (linear regime only). + + A thalamic relay neuron fires low-threshold Ca²⁺ spikes in response to even + a small chirp around rest, so ``analyze_impedance`` declines to report a + profile. + """ + profile, voltage = _run_chirp_profile(make_thalamic_relay(), amplitude=1.0) + + assert voltage.max() > -20.0 # the cell spiked + assert profile is None diff --git a/tests/unit/test_impedance.py b/tests/unit/test_impedance.py index b1c8920c..8c21f35f 100644 --- a/tests/unit/test_impedance.py +++ b/tests/unit/test_impedance.py @@ -20,13 +20,15 @@ _DURATION_MS = 4000.0 _F_START = 1.0 _F_END = 100.0 +_AMP = 0.05 # Small chirp amplitude so |Z|·amplitude stays subthreshold. +_BASELINE_MV = -65.0 # Physiological resting potential for the response trace. def _chirp_and_time( duration_ms: float = _DURATION_MS, f_start: float = _F_START, f_end: float = _F_END, - amplitude: float = 1.0, + amplitude: float = _AMP, dc_offset: float = 0.0, ) -> tuple[np.ndarray, np.ndarray]: """Build a chirp stimulus and matching time axis at ``_FS_HZ``. @@ -66,19 +68,32 @@ def _apply_transfer_function(current: np.ndarray, transfer: np.ndarray) -> np.nd return np.fft.irfft(transfer * np.fft.rfft(current), n=current.size) +def _window_end(time: np.ndarray) -> float: + """Return a window-end time that includes the final sample of ``time``. + + Args: + time: Ascending time axis in ms. + + Returns: + A value slightly past ``time[-1]`` so the ``t < stim_end`` mask keeps + every sample. + """ + return float(time[-1]) + _DT_MS + + # --------------------------------------------------------------------------- # Tests # --------------------------------------------------------------------------- def test_pure_resistor_flat_magnitude() -> None: - """A pure resistor (V = R·I) yields flat |Z| ≈ R and zero phase.""" + """A pure resistor (V = R·I + V_rest) yields flat |Z| ≈ R and zero phase.""" resistance = 50.0 # kΩ·cm² time, current = _chirp_and_time() - voltage = resistance * current + voltage = _BASELINE_MV + resistance * current result = analyze_impedance( - time, voltage, current, 0.0, time[-1] + _DT_MS, _F_START, _F_END + time, voltage, current, 0.0, _window_end(time), _F_START, _F_END ) assert result is not None @@ -94,11 +109,12 @@ def test_pure_resistor_flat_magnitude() -> None: def test_resistor_with_dc_offset_subtracted() -> None: """DC offsets on V and I are removed before the FFT, leaving |Z| ≈ R.""" resistance = 30.0 # kΩ·cm² - time, current = _chirp_and_time(dc_offset=5.0) - voltage = resistance * current + 12.0 + time, chirp_ac = _chirp_and_time() + current = chirp_ac + 5.0 # I carries a DC offset of 5.0. + voltage = _BASELINE_MV + resistance * chirp_ac # V carries a different offset. result = analyze_impedance( - time, voltage, current, 0.0, time[-1] + _DT_MS, _F_START, _F_END + time, voltage, current, 0.0, _window_end(time), _F_START, _F_END ) assert result is not None @@ -112,10 +128,10 @@ def test_first_order_lowpass_no_interior_peak() -> None: time, current = _chirp_and_time() freqs = np.fft.rfftfreq(current.size, d=1.0 / _FS_HZ) transfer = 1.0 / (1.0 + 1j * 2.0 * np.pi * freqs * tau_s) - voltage = _apply_transfer_function(current, transfer) + voltage = _BASELINE_MV + _apply_transfer_function(current, transfer) result = analyze_impedance( - time, voltage, current, 0.0, time[-1] + _DT_MS, _F_START, _F_END + time, voltage, current, 0.0, _window_end(time), _F_START, _F_END ) assert result is not None @@ -134,10 +150,10 @@ def test_synthetic_resonance_detects_f_r() -> None: freqs = np.fft.rfftfreq(current.size, d=1.0 / _FS_HZ) ratio = freqs / f0 transfer = 1.0 / (1.0 - ratio**2 + 1j * ratio / q0) - voltage = _apply_transfer_function(current, transfer) + voltage = _BASELINE_MV + _apply_transfer_function(current, transfer) result = analyze_impedance( - time, voltage, current, 0.0, time[-1] + _DT_MS, _F_START, _F_END + time, voltage, current, 0.0, _window_end(time), _F_START, _F_END ) assert result is not None @@ -154,10 +170,10 @@ def test_synthetic_resonance_detects_f_r() -> None: def test_band_restriction() -> None: """The analysis band is clipped to [f_start, f_end].""" time, current = _chirp_and_time() - voltage = 10.0 * current + voltage = _BASELINE_MV + 10.0 * current result = analyze_impedance( - time, voltage, current, 0.0, time[-1] + _DT_MS, 10.0, 50.0 + time, voltage, current, 0.0, _window_end(time), 10.0, 50.0 ) assert result is not None @@ -171,7 +187,7 @@ def test_band_restriction() -> None: def test_returns_none_short_window() -> None: """A chirp window shorter than the minimum returns None.""" time, current = _chirp_and_time() - voltage = 10.0 * current + voltage = _BASELINE_MV + 10.0 * current result = analyze_impedance(time, voltage, current, 0.0, 10.0, _F_START, _F_END) @@ -181,10 +197,10 @@ def test_returns_none_short_window() -> None: def test_returns_none_invalid_band() -> None: """A band with f_end <= f_start returns None.""" time, current = _chirp_and_time() - voltage = 10.0 * current + voltage = _BASELINE_MV + 10.0 * current result = analyze_impedance( - time, voltage, current, 0.0, time[-1] + _DT_MS, 50.0, 50.0 + time, voltage, current, 0.0, _window_end(time), 50.0, 50.0 ) assert result is None @@ -193,11 +209,11 @@ def test_returns_none_invalid_band() -> None: def test_returns_none_with_nans() -> None: """Non-finite values in the voltage trace return None.""" time, current = _chirp_and_time() - voltage = 10.0 * current + voltage = _BASELINE_MV + 10.0 * current voltage[100] = np.nan result = analyze_impedance( - time, voltage, current, 0.0, time[-1] + _DT_MS, _F_START, _F_END + time, voltage, current, 0.0, _window_end(time), _F_START, _F_END ) assert result is None @@ -206,7 +222,7 @@ def test_returns_none_with_nans() -> None: def test_returns_none_too_few_band_bins() -> None: """A band wide enough in Hz but with too few FFT bins returns None.""" time, current = _chirp_and_time(duration_ms=200.0) - voltage = 10.0 * current + voltage = _BASELINE_MV + 10.0 * current # 200 ms at 1 kHz → 5 Hz frequency resolution → a 1 Hz band holds < 8 bins. result = analyze_impedance(time, voltage, current, 0.0, 200.0, 49.0, 50.0) @@ -214,18 +230,33 @@ def test_returns_none_too_few_band_bins() -> None: assert result is None +def test_returns_none_when_suprathreshold() -> None: + """A response that crosses spike threshold returns None (impedance undefined).""" + time, current = _chirp_and_time(amplitude=1.0) + # A short +30 mV depolarization with a sharp upstroke reads as a spike. + voltage = np.full_like(time, _BASELINE_MV) + spike_window = (time >= 1000.0) & (time < 1005.0) + voltage[spike_window] = 30.0 + + result = analyze_impedance( + time, voltage, current, 0.0, _window_end(time), _F_START, _F_END + ) + + assert result is None + + def test_absolute_conversion_with_area() -> None: """Supplying area_cm2 fills the absolute MΩ spectrum; omitting it leaves it None.""" resistance = 40.0 # kΩ·cm² area_cm2 = 2.0e-5 time, current = _chirp_and_time() - voltage = resistance * current + voltage = _BASELINE_MV + resistance * current with_area = analyze_impedance( - time, voltage, current, 0.0, time[-1] + _DT_MS, _F_START, _F_END, area_cm2 + time, voltage, current, 0.0, _window_end(time), _F_START, _F_END, area_cm2 ) without_area = analyze_impedance( - time, voltage, current, 0.0, time[-1] + _DT_MS, _F_START, _F_END + time, voltage, current, 0.0, _window_end(time), _F_START, _F_END ) assert with_area is not None and without_area is not None From 61ef32748f251606f9f50f59e68834e85f5dd51c Mon Sep 17 00:00:00 2001 From: JCorson Date: Mon, 11 May 2026 17:19:30 -0500 Subject: [PATCH 4/9] feat(analysis): recover impedance from spike-free chirp sub-windows MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit When the chirp window contains spikes, analyze_impedance now falls back to the longest contiguous spike-free segment instead of giving up. A linear chirp's instantaneous frequency rises with time, so the recovered profile covers only a sub-band of the requested [f_start, f_end] — the existing band-edge trimming (_MIN_STIM_FRAC) and minimum-bin check (_MIN_BAND_BINS) already discard low-power bins, so almost no new validation logic is needed. The covered band is reported via frequencies[0]/frequencies[-1] and the sub-window duration via the new ImpedanceProfile.analyzed_window_ms field. Spike-free-segment lookup lives in a new patch_sim.analysis.passive_properties.longest_subthreshold_run helper next to is_subthreshold; it runs analyze_aps and excises a guard-padded window around each detected spike (_SPIKE_GUARD_PRE_MS/_SPIKE_GUARD_POST_MS) to clear the after-hyperpolarization, then returns the widest remaining gap. The cheap pre-FFT guard checks are refactored into a private _impedance_unavailable_reason so the *condition + message* is defined once, and a public wrapper impedance_unavailable_reason is exposed for the UI layer to surface a specific reason when analysis bails out (rather than the previous generic placeholder). --- patch_sim/__init__.py | 4 + patch_sim/analysis/__init__.py | 9 +- patch_sim/analysis/impedance.py | 199 +++++++++++++++++++---- patch_sim/analysis/passive_properties.py | 90 ++++++++++ 4 files changed, 265 insertions(+), 37 deletions(-) diff --git a/patch_sim/__init__.py b/patch_sim/__init__.py index b1dd0253..decc9250 100644 --- a/patch_sim/__init__.py +++ b/patch_sim/__init__.py @@ -56,7 +56,9 @@ density_to_absolute_r_in, double_exp_decay, estimate_rheobase, + impedance_unavailable_reason, is_subthreshold, + longest_subthreshold_run, run_membrane_test, single_exp_decay, single_exp_rise, @@ -154,7 +156,9 @@ "density_to_absolute_r_in", "double_exp_decay", "estimate_rheobase", + "impedance_unavailable_reason", "is_subthreshold", + "longest_subthreshold_run", "run_membrane_test", "single_exp_decay", "single_exp_rise", diff --git a/patch_sim/analysis/__init__.py b/patch_sim/analysis/__init__.py index 9a5d81a3..fa555502 100644 --- a/patch_sim/analysis/__init__.py +++ b/patch_sim/analysis/__init__.py @@ -56,7 +56,11 @@ analyze_hyperpolarization, compute_sag_point, ) -from .impedance import ImpedanceProfile, analyze_impedance +from .impedance import ( + ImpedanceProfile, + analyze_impedance, + impedance_unavailable_reason, +) from .iv_curve import IVAnalysisResult, IVPoint, analyze_iv, compute_iv_point from .membrane_test import ( MEMBRANE_TEST_CURRENT, @@ -71,6 +75,7 @@ density_to_absolute_c_m, density_to_absolute_r_in, is_subthreshold, + longest_subthreshold_run, ) from .sfa import SFAAnalysisResult, SFACurve, analyze_sfa, compute_sfa from .tau_v import ( @@ -111,7 +116,9 @@ "density_to_absolute_r_in", "double_exp_decay", "estimate_rheobase", + "impedance_unavailable_reason", "is_subthreshold", + "longest_subthreshold_run", "run_membrane_test", "single_exp_decay", "single_exp_rise", diff --git a/patch_sim/analysis/impedance.py b/patch_sim/analysis/impedance.py index 009caf74..56eee8bd 100644 --- a/patch_sim/analysis/impedance.py +++ b/patch_sim/analysis/impedance.py @@ -21,6 +21,7 @@ from patch_sim.analysis.passive_properties import ( density_to_absolute_r_in, is_subthreshold, + longest_subthreshold_run, ) logger = logging.getLogger(__name__) @@ -84,8 +85,15 @@ class ImpedanceProfile: supplied and an interior peak exists; ``None`` otherwise. area_cm2: Membrane surface area in cm² used for the absolute conversion. ``None`` when only per-area values were requested. - f_start: Lower edge of the analysis band in Hz. - f_end: Upper edge of the analysis band in Hz. + f_start: Lower edge of the requested analysis band in Hz. + f_end: Upper edge of the requested analysis band in Hz. + analyzed_window_ms: ``None`` when the full chirp window was analyzed. + Otherwise the duration in ms of the spike-free sub-window that was + analyzed instead (the chirp contained spikes; impedance was + recovered from the longest contiguous spike-free segment, which + covers only a sub-band of ``[f_start, f_end]`` — see + ``frequencies[0]`` / ``frequencies[-1]`` for the actually-covered + band). """ frequencies: list[float] @@ -99,6 +107,7 @@ class ImpedanceProfile: area_cm2: float | None = None f_start: float = 0.0 f_end: float = 0.0 + analyzed_window_ms: float | None = None def _half_power_crossing( @@ -172,6 +181,120 @@ def _quality_factor( return float(frequencies[peak_idx]) / fwhm +def _impedance_unavailable_reason( + time: np.ndarray, + voltage: np.ndarray, + injected_current: np.ndarray, + stim_start_ms: float, + stim_end_ms: float, + f_start: float, + f_end: float, +) -> str | None: + """Return a human-readable reason why impedance analysis can't proceed. + + Covers the cheap pre-FFT failure modes — shape / NaN / band sanity, a + too-short chirp window, and a suprathreshold response with no spike-free + fallback segment long enough for the FFT. The post-FFT failures (no + stimulus power in the band, too few in-band bins after trimming) are not + detected here and must be diagnosed by running :func:`analyze_impedance` + itself. + + Args: + time: Time axis in ms. + voltage: Membrane voltage response in mV. + injected_current: Injected chirp current in µA/cm². + stim_start_ms: Start of the chirp stimulus window in ms. + stim_end_ms: End of the chirp stimulus window in ms. + f_start: Starting frequency of the chirp sweep in Hz. + f_end: Ending frequency of the chirp sweep in Hz. + + Returns: + ``None`` when the preamble checks pass and the analysis can proceed; a + single-sentence reason string otherwise. + """ + t = np.asarray(time, dtype=float) + v = np.asarray(voltage, dtype=float) + i = np.asarray(injected_current, dtype=float) + if v.shape != t.shape or i.shape != t.shape: + return "internal: trace length mismatch." + if not (np.all(np.isfinite(v)) and np.all(np.isfinite(i))): + return "internal: trace contains non-finite values." + if f_end <= f_start or f_start < 0.0: + return f"invalid frequency band [{f_start}, {f_end}] Hz." + window_ms = stim_end_ms - stim_start_ms + if window_ms < _MIN_WINDOW_MS: + return ( + f"the chirp window ({window_ms:.1f} ms) is shorter than the " + f"minimum {_MIN_WINDOW_MS:.1f} ms needed for a meaningful spectrum." + ) + mask = (t >= stim_start_ms) & (t < stim_end_ms) + if int(mask.sum()) < 4: + return "too few samples in the chirp window." + t_win, v_win = t[mask], v[mask] + if is_subthreshold(t_win, v_win): + return None + run = longest_subthreshold_run(t_win, v_win) + if run is None: + return ( + "the cell fired throughout the chirp — reduce the amplitude or " + "apply a hyperpolarizing holding current." + ) + lo, hi = run + span = float(t_win[hi - 1] - t_win[lo]) + if span < _MIN_WINDOW_MS: + return ( + f"the longest spike-free segment ({span:.1f} ms) is shorter than the " + f"minimum {_MIN_WINDOW_MS:.1f} ms — reduce the amplitude or apply a " + "hyperpolarizing holding current." + ) + return None + + +def impedance_unavailable_reason( + time: np.ndarray, + voltage: np.ndarray, + injected_current: np.ndarray, + stim_start_ms: float, + stim_end_ms: float, + f_start: float, + f_end: float, +) -> str: + """Public wrapper around :func:`_impedance_unavailable_reason`. + + Returns the empty string when the preamble checks pass (and the analysis + would proceed past the windowing / suprathreshold guards), so callers can + use the result directly as UI copy. Only the cheap pre-FFT failure modes + are reported — when this returns ``""`` and :func:`analyze_impedance` still + returns ``None``, the caller should fall back to a generic + "too little usable signal in the band" message. + + Args: + time: Time axis in ms. + voltage: Membrane voltage response in mV. + injected_current: Injected chirp current in µA/cm². + stim_start_ms: Start of the chirp stimulus window in ms. + stim_end_ms: End of the chirp stimulus window in ms. + f_start: Starting frequency of the chirp sweep in Hz. + f_end: Ending frequency of the chirp sweep in Hz. + + Returns: + The reason string from :func:`_impedance_unavailable_reason`, or ``""`` + when that function returns ``None``. + """ + return ( + _impedance_unavailable_reason( + time, + voltage, + injected_current, + stim_start_ms, + stim_end_ms, + f_start, + f_end, + ) + or "" + ) + + def analyze_impedance( time: np.ndarray, voltage: np.ndarray, @@ -192,9 +315,14 @@ def analyze_impedance( ``|Z|``, only when it is an interior maximum that rises above the low-frequency reference) and quality factor. - Membrane impedance is a linear, small-signal quantity: the analysis bails - out (returns ``None``) when the windowed response is suprathreshold, since - spike transients dominate the spectrum and the result would be meaningless. + Membrane impedance is a linear, small-signal quantity, so spike transients + are excluded: when the windowed response contains spikes the analysis falls + back to the longest contiguous spike-free sub-window (a linear chirp's + frequency is time-dependent, so the recovered profile then covers only a + sub-band of ``[f_start, f_end]`` — surfaced via + :attr:`ImpedanceProfile.analyzed_window_ms`). The analysis bails out + (returns ``None``) when no such spike-free segment is long enough — see + :func:`impedance_unavailable_reason` for the user-facing reason. Args: time: Time axis in ms (uniformly sampled). @@ -212,44 +340,42 @@ def analyze_impedance( An :class:`ImpedanceProfile`, or ``None`` when the inputs are inconsistent (mismatched lengths, non-finite values), the band is invalid (``f_end <= f_start`` or ``f_start < 0``), the chirp window is - shorter than :data:`_MIN_WINDOW_MS`, the windowed response is - suprathreshold, the stimulus has no spectral power in the band, or - fewer than :data:`_MIN_BAND_BINS` FFT bins remain after trimming. + shorter than :data:`_MIN_WINDOW_MS`, the response fires throughout the + chirp (no spike-free segment of at least :data:`_MIN_WINDOW_MS`), the + stimulus has no spectral power in the band, or fewer than + :data:`_MIN_BAND_BINS` FFT bins remain after band-edge trimming. """ + reason = _impedance_unavailable_reason( + time, + voltage, + injected_current, + stim_start_ms, + stim_end_ms, + f_start, + f_end, + ) + if reason is not None: + logger.warning("analyze_impedance: %s", reason) + return None + t = np.asarray(time, dtype=float) v = np.asarray(voltage, dtype=float) i = np.asarray(injected_current, dtype=float) - if v.shape != t.shape or i.shape != t.shape: - logger.warning("analyze_impedance: trace length mismatch; skipping.") - return None - if not (np.all(np.isfinite(v)) and np.all(np.isfinite(i))): - logger.warning("analyze_impedance: non-finite trace values; skipping.") - return None - if f_end <= f_start or f_start < 0.0: - logger.warning( - "analyze_impedance: invalid band [%s, %s] Hz; skipping.", f_start, f_end - ) - return None - if stim_end_ms - stim_start_ms < _MIN_WINDOW_MS: - logger.warning( - "analyze_impedance: chirp window %.1f ms shorter than %.1f ms; skipping.", - stim_end_ms - stim_start_ms, - _MIN_WINDOW_MS, - ) - return None - mask = (t >= stim_start_ms) & (t < stim_end_ms) - if int(mask.sum()) < 4: - logger.warning("analyze_impedance: too few samples in chirp window; skipping.") - return None t_win, v_win, i_win = t[mask], v[mask], i[mask] - if not is_subthreshold(t_win, v_win): - logger.warning( - "analyze_impedance: response is suprathreshold; impedance is only " - "defined in the subthreshold regime; skipping." - ) - return None + analyzed_window_ms: float | None + if is_subthreshold(t_win, v_win): + analyzed_window_ms = None + else: + # _impedance_unavailable_reason already verified that a long-enough + # spike-free run exists; fall back to it. + run = longest_subthreshold_run(t_win, v_win) + if run is None: # defensive — guarded above + return None + lo, hi = run + t_win, v_win, i_win = t_win[lo:hi], v_win[lo:hi], i_win[lo:hi] + analyzed_window_ms = float(t_win[-1] - t_win[0]) dt_ms = float(np.mean(np.diff(t_win))) if dt_ms <= 0.0: @@ -348,4 +474,5 @@ def analyze_impedance( area_cm2=area_cm2, f_start=float(f_start), f_end=float(f_end), + analyzed_window_ms=analyzed_window_ms, ) diff --git a/patch_sim/analysis/passive_properties.py b/patch_sim/analysis/passive_properties.py index 302a9ae6..2804a5fa 100644 --- a/patch_sim/analysis/passive_properties.py +++ b/patch_sim/analysis/passive_properties.py @@ -37,6 +37,16 @@ #: Default initial guess for the membrane time constant (ms). _TAU_INITIAL_GUESS_MS: float = 5.0 +#: Guard padding (ms) trimmed *before* a detected spike's threshold crossing +#: when carving out the spike-free portion of a trace. Spike onset is sharp, +#: so a small pad suffices. +_SPIKE_GUARD_PRE_MS: float = 2.0 + +#: Guard padding (ms) trimmed *after* a detected spike's peak. Must be long +#: enough to clear the after-hyperpolarization so the recovered segment carries +#: only the linear subthreshold response. +_SPIKE_GUARD_POST_MS: float = 20.0 + @dataclasses.dataclass class PassiveProperties: @@ -150,6 +160,86 @@ def is_subthreshold( return result.spike_count == 0 +def longest_subthreshold_run( + time: np.ndarray, + voltage: np.ndarray, + *, + dvdt_threshold: float = 20.0, + min_spike_height: float = -20.0, +) -> tuple[int, int] | None: + """Find the longest contiguous spike-free span of a voltage trace. + + Runs :func:`analyze_aps` to locate spikes, excises a guard-padded window + around each one — from :data:`_SPIKE_GUARD_PRE_MS` before its threshold + crossing to :data:`_SPIKE_GUARD_POST_MS` after its peak, so the + after-hyperpolarization transient is removed too — and returns the widest + remaining gap. When the trace contains no spikes the whole trace is + returned. + + Args: + time: Time array in ms, shape ``(N,)``, ascending. + voltage: Membrane voltage in mV, shape ``(N,)``, same length as + ``time``. + dvdt_threshold: dV/dt value (mV/ms) used to identify spike onset. + min_spike_height: Minimum peak voltage (mV) for a candidate event to be + classified as a spike. + + Returns: + Half-open index bounds ``(start_idx, stop_idx)`` into ``time`` / + ``voltage`` of the longest spike-free span (``(0, N)`` when there are no + spikes), or ``None`` when no spike-free span exists — every sample falls + inside a guard window. + """ + time = np.asarray(time, dtype=float) + voltage = np.asarray(voltage, dtype=float) + n = time.size + if n == 0: + return None + + result = analyze_aps( + time, + voltage, + dvdt_threshold=dvdt_threshold, + min_spike_height=min_spike_height, + ) + if result.spike_count == 0: + return (0, n) + + # Merge the guard-padded "bad" time intervals around each spike. + intervals = sorted( + (s.threshold_time - _SPIKE_GUARD_PRE_MS, s.peak_time + _SPIKE_GUARD_POST_MS) + for s in result.spikes + ) + merged: list[list[float]] = [list(intervals[0])] + for lo, hi in intervals[1:]: + if lo <= merged[-1][1]: + merged[-1][1] = max(merged[-1][1], hi) + else: + merged.append([lo, hi]) + + # Spike-free time gaps: before the first interval, between consecutive + # intervals, and after the last interval. + gaps: list[tuple[float, float]] = [(float(time[0]), merged[0][0])] + for idx in range(len(merged) - 1): + gaps.append((merged[idx][1], merged[idx + 1][0])) + gaps.append((merged[-1][1], float(time[-1]))) + + best: tuple[int, int] | None = None + best_span = 0.0 + for gap_start, gap_end in gaps: + if gap_end <= gap_start: + continue + start_idx = int(np.searchsorted(time, gap_start, side="left")) + stop_idx = int(np.searchsorted(time, gap_end, side="right")) + if stop_idx - start_idx < 2: + continue + span = float(time[stop_idx - 1] - time[start_idx]) + if span > best_span: + best_span = span + best = (start_idx, stop_idx) + return best + + def _exponential_model(t: np.ndarray, v_ss: float, a: float, tau: float) -> np.ndarray: """Single-exponential model for the voltage step response. From 82eedc835cd7ea2f1a99cb302361f6381f81adc3 Mon Sep 17 00:00:00 2001 From: JCorson Date: Mon, 11 May 2026 17:19:42 -0500 Subject: [PATCH 5/9] feat(ui): surface impedance unavailability reason and sub-window caption MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The analysis sidebar's impedance pane now distinguishes the failure modes of analyze_impedance: - When the analysis recovers a profile from a spike-free *sub*-window (rather than the full chirp), the pane shows a small caption above the plot — Computed from a NNN ms spike-free segment of the chirp (X.X–Y.Y Hz). — so the user can see the band was narrowed. - When the analysis bails out entirely, the pane shows the specific reason returned by patch_sim.impedance_unavailable_reason (e.g. "the cell fired throughout the chirp — reduce the amplitude or apply a hyperpolarizing holding current") instead of the generic "reduce the amplitude" placeholder. The reason + caption ride inside the existing impedance_data dict (unavailable_reason / caption keys), so no _SimResult or sweep-executor changes are needed. AnalysisState.has_impedance_data now keys off the presence of "frequencies" so a reason-only dict doesn't try to render a figure. --- .../components/metrics/impedance_panel.py | 57 +++++++++++++++---- patch_sim_ui/state/_analysis_format.py | 38 +++++++++++-- patch_sim_ui/state/analysis.py | 30 +++++++++- 3 files changed, 104 insertions(+), 21 deletions(-) diff --git a/patch_sim_ui/components/metrics/impedance_panel.py b/patch_sim_ui/components/metrics/impedance_panel.py index 736f4c90..86865813 100644 --- a/patch_sim_ui/components/metrics/impedance_panel.py +++ b/patch_sim_ui/components/metrics/impedance_panel.py @@ -42,10 +42,26 @@ def _impedance_summary() -> rx.Component: def _impedance_plot() -> rx.Component: """Render the embedded impedance-profile plot (|Z| and phase vs frequency). + Includes a sub-window caption when ``analyze_impedance`` recovered the + profile from a spike-free segment rather than the full chirp window. + Returns: - A flex container holding the impedance Plotly figure. + A flex container holding the optional caption and the impedance Plotly + figure. """ return rx.flex( + rx.cond( + AnalysisState.impedance_caption != "", + rx.text( + AnalysisState.impedance_caption, + size="1", + color="gray", + text_align="center", + padding_x="3", + padding_top="2", + ), + rx.fragment(), + ), rx.plotly( data=AnalysisState.impedance_figure, width="100%", @@ -62,7 +78,8 @@ def impedance_tab() -> rx.Component: """Render the impedance-analysis tab for the chirp current-clamp pane. Shows the resonance summary and the |Z| / phase plot when impedance data - is available, and a placeholder otherwise. + is available, the specific unavailability reason when the analysis bailed + out, and a generic placeholder otherwise. Returns: The impedance-analysis content as a scrollable flex column. @@ -79,16 +96,32 @@ def impedance_tab() -> rx.Component: height="100%", width="100%", ), - rx.flex( - rx.text( - "Run a Frequency Response (chirp) current clamp simulation to " - "see the impedance profile. The chirp must keep the cell " - "subthreshold — reduce the amplitude if the response spikes.", - size="1", - color="gray", - text_align="center", + rx.cond( + AnalysisState.impedance_unavailable_reason != "", + rx.flex( + rx.text( + "Impedance not computed: " + + AnalysisState.impedance_unavailable_reason, + size="1", + color="gray", + text_align="center", + ), + padding="4", + justify="center", + ), + rx.flex( + rx.text( + "Run a Frequency Response (chirp) current clamp simulation " + "to see the impedance profile. The chirp must keep the " + "cell subthreshold — reduce the amplitude and, for cells " + "that fire spontaneously, apply a hyperpolarizing DC " + "offset (holding current).", + size="1", + color="gray", + text_align="center", + ), + padding="4", + justify="center", ), - padding="4", - justify="center", ), ) diff --git a/patch_sim_ui/state/_analysis_format.py b/patch_sim_ui/state/_analysis_format.py index bfee1a83..c8abf14e 100644 --- a/patch_sim_ui/state/_analysis_format.py +++ b/patch_sim_ui/state/_analysis_format.py @@ -492,10 +492,14 @@ def _compute_impedance_data( area_cm2: Membrane surface area in cm², or ``None`` when undeclared. Returns: - A dict with keys ``frequencies``, ``magnitude``, ``phase`` (lists), - ``resonance_frequency``, ``quality_factor``, ``peak_impedance`` - (pre-formatted strings), and ``units``. Returns an empty dict when the - impedance analysis is not applicable (see :func:`patch_sim.analyze_impedance`). + On success, a dict with keys ``frequencies``, ``magnitude``, ``phase`` + (lists), ``resonance_frequency``, ``quality_factor``, ``peak_impedance`` + (pre-formatted strings), and ``units`` — plus a ``caption`` string when + the profile was recovered from a spike-free sub-window rather than the + full chirp window. On failure, a dict with a single + ``unavailable_reason`` key carrying a human-readable sentence. Returns + an empty dict only when the sweep itself is too short to attempt + analysis. """ time_arr = np.asarray(sweep.time, dtype=float) if len(time_arr) < 2: @@ -519,7 +523,22 @@ def _compute_impedance_data( area_cm2=area_cm2, ) if profile is None: - return {} + # impedance_unavailable_reason covers the cheap pre-FFT failures; if it + # is silent, the FFT itself rejected the band (no stimulus power, too + # few in-band bins after trimming) — fall back to a generic message. + reason = ( + patch_sim.impedance_unavailable_reason( + time_arr, + voltage_arr, + current_arr, + stim_start_ms, + stim_end_ms, + start_frequency, + end_frequency, + ) + or "the chirp produced too little usable signal in the swept band." + ) + return {"unavailable_reason": reason} if profile.magnitude_mohm is not None: magnitude = profile.magnitude_mohm peak = profile.peak_impedance_mohm @@ -528,7 +547,7 @@ def _compute_impedance_data( magnitude = profile.magnitude peak = profile.peak_impedance units = "kΩ·cm²" - return { + result: dict[str, Any] = { "frequencies": profile.frequencies, "magnitude": magnitude, "phase": profile.phase, @@ -537,6 +556,13 @@ def _compute_impedance_data( "peak_impedance": _fmt_optional(peak, ".3g"), "units": units, } + if profile.analyzed_window_ms is not None and profile.frequencies: + result["caption"] = ( + f"Computed from a {profile.analyzed_window_ms:.0f} ms spike-free " + f"segment of the chirp ({profile.frequencies[0]:.1f}–" + f"{profile.frequencies[-1]:.1f} Hz)." + ) + return result def _compute_cc_multi_sweep_analysis( diff --git a/patch_sim_ui/state/analysis.py b/patch_sim_ui/state/analysis.py index 325b6cc5..429210ff 100644 --- a/patch_sim_ui/state/analysis.py +++ b/patch_sim_ui/state/analysis.py @@ -221,8 +221,13 @@ def phase_plane_figure(self) -> go.Figure: @rx.var def has_impedance_data(self) -> bool: - """Return True when chirp impedance-profile results are available.""" - return len(self.impedance_data) > 0 + """Return True when a renderable chirp impedance-profile is available. + + A dict carrying only an ``unavailable_reason`` (the analysis bailed out) + is intentionally treated as *no* impedance data so the panel falls + through to the reason-aware placeholder. + """ + return "frequencies" in self.impedance_data @rx.var def impedance_figure(self) -> go.Figure: @@ -230,6 +235,25 @@ def impedance_figure(self) -> go.Figure: Returns an empty figure when no impedance data is available. """ - if not self.impedance_data: + if "frequencies" not in self.impedance_data: return go.Figure() return build_impedance_figure(self.impedance_data) + + @rx.var + def impedance_unavailable_reason(self) -> str: + """Return the reason the chirp impedance analysis didn't run, or ``""``. + + Set by :func:`patch_sim_ui.state._analysis_format._compute_impedance_data` + when :func:`patch_sim.analyze_impedance` returns ``None``; the impedance + panel renders this in place of the generic placeholder. + """ + return self.impedance_data.get("unavailable_reason", "") + + @rx.var + def impedance_caption(self) -> str: + """Return the sub-window caption for the impedance plot, or ``""``. + + Populated when ``analyze_impedance`` recovered the profile from a + spike-free sub-window rather than the full chirp window. + """ + return self.impedance_data.get("caption", "") From a93ac367e3c1fd8422ea97156008723f29c436ab Mon Sep 17 00:00:00 2001 From: JCorson Date: Mon, 11 May 2026 17:19:52 -0500 Subject: [PATCH 6/9] feat(presets): tune FREQUENCY_RESPONSE so the chirp runs out of the box MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The previous default (`dc_offset=8`, `amplitude=4`) drove almost every preset above spike threshold, so the FFT-based impedance analysis returned nothing for the exact presets the PR test plan recommends (Thalamic Relay / STN / CA1). Soften the global default to `dc_offset=0`, `amplitude=0.25` µA/cm², and extend `stimulus_duration` to 1000 ms — small enough to stay in the linear regime for every quiescent-at-rest preset (Squid, FSI, Cortical, Thalamic Relay¹, CA1), long enough to resolve the few-Hz Ih resonance. For autonomous pacemakers (Purkinje, SNc Dopaminergic, STN, TRN) the softened default still triggers tonic firing, so each preset's `PROTOCOL_ADJUSTMENTS[FREQUENCY_RESPONSE]` now supplies a small hyperpolarizing holding current (`dc_offset` of −0.5 to −5 µA/cm²) that silences spontaneous firing during the chirp — mirroring how a real ZAP/impedance experiment is run on a pacemaker cell. Values chosen empirically via `scratch/tune_frequency_response.py` so every preset's chirp response stays subthreshold across the full sweep and `analyze_impedance` returns a full-window profile. ¹ Thalamic Relay also needs a small hyperpolarizing hold because its slow-inactivating ICaT gives it a very low LTS rheobase that the chirp would otherwise trip on each depolarizing swing. --- patch_sim/presets/dopaminergic.py | 10 ++++++++++ patch_sim/presets/protocols.py | 20 ++++++++++++-------- patch_sim/presets/purkinje.py | 11 +++++++++++ patch_sim/presets/stn.py | 10 ++++++++++ patch_sim/presets/thalamocortical.py | 12 ++++++++++++ patch_sim/presets/trn.py | 11 +++++++++++ 6 files changed, 66 insertions(+), 8 deletions(-) diff --git a/patch_sim/presets/dopaminergic.py b/patch_sim/presets/dopaminergic.py index 3c1f554c..b591661e 100644 --- a/patch_sim/presets/dopaminergic.py +++ b/patch_sim/presets/dopaminergic.py @@ -16,6 +16,7 @@ from patch_sim.constants import ( ACTION_POTENTIAL, FI_CURVE, + FREQUENCY_RESPONSE, HYPERPOLARIZATION_STEPS, REPETITIVE_FIRING, SUBTHRESHOLD_RESPONSE, @@ -74,6 +75,15 @@ "max_stimulus": -5.0, "stimulus_step": 5.0, }, + # Autonomous pacemaker: a hyperpolarizing holding current is required + # to silence tonic firing during the chirp (mirrors a real ZAP + # experiment under current-clamp hold). −2 µA/cm² holds the cell at + # ≈ −72 mV (below the INaP/Cav1.3 oscillator window); amp=0.25 keeps + # the response in the linear regime. + FREQUENCY_RESPONSE: { + "dc_offset": -2.0, + "amplitude": 0.25, + }, } diff --git a/patch_sim/presets/protocols.py b/patch_sim/presets/protocols.py index 2da38237..63385c09 100644 --- a/patch_sim/presets/protocols.py +++ b/patch_sim/presets/protocols.py @@ -97,19 +97,23 @@ "max_stimulus": 60.0, "stimulus_step": 10.0, }, - # NOTE: this DC offset / amplitude drives most presets above spike - # threshold, so the FFT-based impedance analysis (valid only in the linear - # subthreshold regime) reports nothing for those cells. To measure a clean - # impedance / subthreshold-resonance profile, reduce the amplitude (and zero - # the DC offset) so the chirp response stays below threshold. + # 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 + # few-mV peak-to-peak swing, small enough for the linear regime and large + # enough for usable SNR; 1000 ms is long enough to resolve the few-Hz Ih + # resonance band. Autonomous pacemakers (Purkinje, SNc DA, STN, TRN) and + # the burst-mode Thalamic Relay still need a hyperpolarizing holding + # current; their per-neuron PROTOCOL_ADJUSTMENTS[FREQUENCY_RESPONSE] + # overrides supply it. FREQUENCY_RESPONSE: { "clamp_mode": CURRENT_CLAMP, "protocol_type": "Chirp", "pre_stimulus_duration": 0.0, - "stimulus_duration": 500.0, + "stimulus_duration": 1000.0, "post_stimulus_duration": 0.0, - "dc_offset": 8.0, - "amplitude": 4.0, + "dc_offset": 0.0, + "amplitude": 0.25, "start_frequency": 1.0, "end_frequency": 100.0, }, diff --git a/patch_sim/presets/purkinje.py b/patch_sim/presets/purkinje.py index fc91900d..7c0a33d5 100644 --- a/patch_sim/presets/purkinje.py +++ b/patch_sim/presets/purkinje.py @@ -16,6 +16,7 @@ ) from patch_sim.constants import ( ACTION_POTENTIAL, + FREQUENCY_RESPONSE, HYPERPOLARIZATION_STEPS, REPETITIVE_FIRING, SUBTHRESHOLD_RESPONSE, @@ -55,6 +56,16 @@ "max_stimulus": -0.5, "stimulus_step": 0.5, }, + # Autonomous pacemaker: a hyperpolarizing holding current is required + # to silence tonic firing during the chirp (mirrors a real ZAP + # experiment under current-clamp hold). −5 µA/cm² holds the cell at + # ≈ −80 mV (below the INaP-driven oscillator window); amp=0.25 keeps + # the response in the linear regime so the FFT-based impedance + # analysis reports a clean profile. + FREQUENCY_RESPONSE: { + "dc_offset": -5.0, + "amplitude": 0.25, + }, } diff --git a/patch_sim/presets/stn.py b/patch_sim/presets/stn.py index 2ac43e5b..0f2f380a 100644 --- a/patch_sim/presets/stn.py +++ b/patch_sim/presets/stn.py @@ -18,6 +18,7 @@ from patch_sim.constants import ( ACTION_POTENTIAL, FI_CURVE, + FREQUENCY_RESPONSE, REPETITIVE_FIRING, SUBTHRESHOLD_RESPONSE, ) @@ -56,6 +57,15 @@ # Very high ICaT conductance (g=5.0 mS/cm²) produces a prominent # post-inhibitory rebound burst on step release; Ih (g=1.0 mS/cm²) # adds a depolarizing overshoot that can trigger additional spikes. + # Autonomous pacemaker: a hyperpolarizing holding current silences + # tonic firing during the chirp (mirrors a real ZAP experiment under + # current-clamp hold). −2 µA/cm² holds the cell at ≈ −69 mV (below + # the INaP-driven oscillator window); amp=0.25 keeps the response in + # the linear regime. + FREQUENCY_RESPONSE: { + "dc_offset": -2.0, + "amplitude": 0.25, + }, } diff --git a/patch_sim/presets/thalamocortical.py b/patch_sim/presets/thalamocortical.py index f23b7b25..2e94645d 100644 --- a/patch_sim/presets/thalamocortical.py +++ b/patch_sim/presets/thalamocortical.py @@ -13,6 +13,7 @@ from patch_sim.constants import ( ACTION_POTENTIAL, FI_CURVE, + FREQUENCY_RESPONSE, REPETITIVE_FIRING, SUBTHRESHOLD_RESPONSE, ) @@ -62,6 +63,17 @@ # (g=2.5 mS/cm²); on release a textbook post-inhibitory LTS burst # fires (issue #287; McCormick & Huguenard 1992). Ih (g=1.0 mS/cm²) # also contributes via sag and post-step overshoot. + # Very low LTS rheobase (≈ 0.05 µA/cm² with slow-inactivating ICaT, + # see SUBTHRESHOLD_RESPONSE) makes the global chirp default fire an + # LTS even at amp=0.1, so a small hyperpolarizing hold is applied. + # −1 µA/cm² parks the cell near −71 mV (further de-inactivating ICaT + # without crossing the LTS threshold during the chirp swing); amp=0.25 + # then yields a ≈ 5 mV peak-to-peak subthreshold response over the + # full chirp window. + FREQUENCY_RESPONSE: { + "dc_offset": -1.0, + "amplitude": 0.25, + }, } diff --git a/patch_sim/presets/trn.py b/patch_sim/presets/trn.py index bdc0e327..97bc0de7 100644 --- a/patch_sim/presets/trn.py +++ b/patch_sim/presets/trn.py @@ -15,6 +15,7 @@ from patch_sim.constants import ( ACTION_POTENTIAL, FI_CURVE, + FREQUENCY_RESPONSE, HYPERPOLARIZATION_STEPS, REPETITIVE_FIRING, SUBTHRESHOLD_RESPONSE, @@ -68,6 +69,16 @@ "max_stimulus": -1.0, "stimulus_step": 1.0, }, + # Autonomous pacemaker: a small hyperpolarizing holding current + # silences the tonic train during the chirp (mirrors a real ZAP + # experiment under current-clamp hold). −0.5 µA/cm² parks the cell + # near −81 mV, slightly below v_rest = −80 mV — far enough to suppress + # firing while still keeping ICaT partially de-inactivated; amp=0.25 + # keeps every chirp swing below the LTS threshold. + FREQUENCY_RESPONSE: { + "dc_offset": -0.5, + "amplitude": 0.25, + }, } From 58242e5f16ff08eb50b2a9b8c906527b7a188180 Mon Sep 17 00:00:00 2001 From: JCorson Date: Mon, 11 May 2026 17:20:08 -0500 Subject: [PATCH 7/9] test: cover spike-free sub-window fallback and tuned FREQUENCY_RESPONSE MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Unit: - tests/unit/test_passive_properties.py: longest_subthreshold_run tests — full trace returned when no spikes, trailing gap returned when one early spike, None when guard windows cover the whole trace. - tests/unit/test_impedance.py: * Renamed test_returns_none_when_suprathreshold → test_returns_none_when_spikes_throughout and rewrote it to tile spikes every 50 ms so every excised guard window butts against the next (a single isolated spike is now tolerated via the sub-window fallback). * test_recovers_profile_from_spike_free_segment: synthetic resonator with one early spike → profile recovered, reported analyzed_window_ms ≈ trace duration minus the spike + guard. * test_impedance_unavailable_reason_messages: public reason wrapper returns "" for OK inputs and the expected sentences for a short window / spikes-throughout trace. Integration: - tests/integration/test_impedance_simulation.py: new parametrized test_frequency_response_preset_yields_profile_for_every_neuron — for every shipped preset, building the chirp via build_protocol_from_preset(FREQUENCY_RESPONSE, neuron_preset=name) and running analyze_impedance yields a non-None profile. Replaces the old test_suprathreshold_chirp_returns_none (whose intent — "a spiking trace returns None" — is now better covered by the unit test, since whether the sub-window fallback kicks in depends on preset-specific kinetics that make a stable integration assertion brittle). - tests/integration/test_preset_protocols.py: lower the "voltage is responding to the chirp" sanity thresholds (std > 0.05 mV, range > 0.2 mV) — the tuned chirp is small by design, so low-R_in cells like the squid axon swing by only a few tenths of a mV. The new parametrized impedance test gives a stronger positive-direction check anyway. E2e: - tests/e2e/test_current_clamp_flow.py: replace test_suprathreshold_chirp_leaves_impedance_empty with test_oversized_chirp_surfaces_unavailable_reason — load the FREQUENCY_RESPONSE preset, crank dc_offset/amplitude high enough to drive sustained spiking, and assert impedance_data carries an unavailable_reason but no frequencies (so the panel falls through to the reason-aware placeholder). --- tests/e2e/test_current_clamp_flow.py | 28 +++++-- .../integration/test_impedance_simulation.py | 50 +++++++++-- tests/integration/test_preset_protocols.py | 20 +++-- tests/unit/test_impedance.py | 84 +++++++++++++++++-- tests/unit/test_passive_properties.py | 60 +++++++++++++ 5 files changed, 212 insertions(+), 30 deletions(-) diff --git a/tests/e2e/test_current_clamp_flow.py b/tests/e2e/test_current_clamp_flow.py index ed3a713b..3fecfbb4 100644 --- a/tests/e2e/test_current_clamp_flow.py +++ b/tests/e2e/test_current_clamp_flow.py @@ -148,17 +148,29 @@ async def test_subthreshold_chirp_populates_impedance(state_tree: StateTree) -> assert len(impedance["frequencies"]) > 0 -async def test_suprathreshold_chirp_leaves_impedance_empty( +async def test_oversized_chirp_surfaces_unavailable_reason( state_tree: StateTree, ) -> None: - """The default (large) Frequency Response chirp drives spiking → no impedance.""" - await run_flow( - state_tree, - neuron_preset=SQUID_GIANT_AXON, - protocol_preset=FREQUENCY_RESPONSE, - ) + """A chirp that drives sustained spiking populates an ``unavailable_reason``. - assert state_tree.analysis.impedance_data == {} + With the tuned global default the squid axon now stays subthreshold, so we + crank the chirp up explicitly to drive a suprathreshold-throughout + response. ``impedance_data`` carries the reason string but no + ``frequencies`` key, so the panel falls through to the reason-aware + placeholder. + """ + async with patch_get_state(state_tree): + [_ async for _ in state_tree.neuron.load_neuron_preset(SQUID_GIANT_AXON)] + [_ async for _ in state_tree.protocol.load_protocol_preset(FREQUENCY_RESPONSE)] + state_tree.protocol.dc_offset = 30.0 + state_tree.protocol.amplitude = 10.0 + + simulate_and_apply(state_tree) + + impedance = state_tree.analysis.impedance_data + assert "frequencies" not in impedance + assert "unavailable_reason" in impedance + assert impedance["unavailable_reason"] # non-empty async def test_action_potential_preset_leaves_impedance_empty( diff --git a/tests/integration/test_impedance_simulation.py b/tests/integration/test_impedance_simulation.py index 5ef61b2a..badc237e 100644 --- a/tests/integration/test_impedance_simulation.py +++ b/tests/integration/test_impedance_simulation.py @@ -7,10 +7,16 @@ """ import numpy as np +import pytest import patch_sim from patch_sim.clamp_simulations import SIM_SAMPLING_FREQ -from patch_sim.presets import make_ca1_pyramidal, make_thalamic_relay +from patch_sim.constants import FREQUENCY_RESPONSE +from patch_sim.presets import ( + NEURON_PRESETS, + make_ca1_pyramidal, +) +from patch_sim.presets.protocols import build_protocol_from_preset _DURATION_MS = 500.0 _F_START = 1.0 @@ -88,14 +94,40 @@ def test_ca1_chirp_shows_resonance_and_high_impedance(hh_model) -> None: assert ca1_profile.magnitude[0] > 5.0 * squid_profile.magnitude[0] -def test_suprathreshold_chirp_returns_none() -> None: - """A chirp that drives spiking yields no impedance profile (linear regime only). +# Integration-level coverage for the "no usable spike-free segment → None" +# branch lives at the unit level (tests/unit/test_impedance.py +# ``test_returns_none_when_spikes_throughout``) because the size/duration +# combination needed to wipe out every 50 ms gap is sensitive to preset +# excitability and gating kinetics, making a stable integration test brittle. +# The parametrized test below gives the cleaner positive-direction coverage — +# every shipped preset yields a profile out of the box. - A thalamic relay neuron fires low-threshold Ca²⁺ spikes in response to even - a small chirp around rest, so ``analyze_impedance`` declines to report a - profile. + +@pytest.mark.parametrize("name", list(NEURON_PRESETS.keys())) +def test_frequency_response_preset_yields_profile_for_every_neuron(name: str) -> None: + """``FREQUENCY_RESPONSE`` + per-neuron overrides recover a profile for every preset. + + Regression guard for the tuned global default plus the per-pacemaker + holding-current overrides: building the chirp via + ``build_protocol_from_preset`` (same code path the UI uses) and running it + through ``analyze_impedance`` must yield a non-``None`` profile with finite, + positive magnitudes for every preset. """ - profile, voltage = _run_chirp_profile(make_thalamic_relay(), amplitude=1.0) + neuron = NEURON_PRESETS[name]() + stim = build_protocol_from_preset( + FREQUENCY_RESPONSE, + neuron_preset=name, + sampling_frequency=SIM_SAMPLING_FREQ, + )[0] + result = patch_sim.simulate_current_clamp(neuron, stim) + time = np.asarray(result["time"]) + voltage = np.asarray(result["voltage"]) + duration = float(time[-1]) + profile = patch_sim.analyze_impedance( + time, voltage, stim, 0.0, duration, _F_START, _F_END, area_cm2=neuron.area_cm2 + ) - assert voltage.max() > -20.0 # the cell spiked - assert profile is None + assert profile is not None, f"{name}: analyze_impedance returned None" + mag = np.asarray(profile.magnitude) + assert mag.size >= 8 + assert np.all(np.isfinite(mag)) and np.all(mag > 0.0) diff --git a/tests/integration/test_preset_protocols.py b/tests/integration/test_preset_protocols.py index b14d651c..515dada8 100644 --- a/tests/integration/test_preset_protocols.py +++ b/tests/integration/test_preset_protocols.py @@ -289,10 +289,14 @@ def test_fi_curve_preset(preset_name: str) -> None: def test_frequency_response_preset(preset_name: str) -> None: """Frequency Response protocol produces a measurable voltage response. - The chirp stimulus sweeps from 1 to 100 Hz and elicits a mixture of - subthreshold and suprathreshold responses. The test asserts that the - voltage trace is not flat: voltage range must exceed 2 mV and standard - deviation must exceed 0.5 mV. + The tuned chirp deliberately uses a small amplitude (0.25 µA/cm² by + default) plus per-pacemaker hyperpolarizing-holding-current overrides so + the response stays in the linear subthreshold regime — small-R_in cells + like the squid axon and FSI swing by only a few tenths of a mV, so the + threshold for "responding" is lenient here. The companion test + ``test_frequency_response_preset_yields_profile_for_every_neuron`` + (tests/integration/test_impedance_simulation.py) provides the stronger + positive-direction check that ``analyze_impedance`` recovers a profile. Args: preset_name: Name of the neuron preset under test. @@ -308,12 +312,12 @@ def test_frequency_response_preset(preset_name: str) -> None: v_std = float(np.std(result["voltage"])) v_range = float(result["voltage"].max() - result["voltage"].min()) - assert v_std > 0.5, ( - f"{preset_name}: Frequency Response voltage std {v_std:.3f} mV ≤ 0.5 mV " + assert v_std > 0.05, ( + f"{preset_name}: Frequency Response voltage std {v_std:.3f} mV ≤ 0.05 mV " "(neuron may not be responding to the chirp stimulus)" ) - assert v_range > 2.0, ( - f"{preset_name}: Frequency Response voltage range {v_range:.3f} mV ≤ 2.0 mV " + assert v_range > 0.2, ( + f"{preset_name}: Frequency Response voltage range {v_range:.3f} mV ≤ 0.2 mV " "(neuron may not be responding to the chirp stimulus)" ) diff --git a/tests/unit/test_impedance.py b/tests/unit/test_impedance.py index 8c21f35f..2a9ddbf0 100644 --- a/tests/unit/test_impedance.py +++ b/tests/unit/test_impedance.py @@ -230,13 +230,21 @@ def test_returns_none_too_few_band_bins() -> None: assert result is None -def test_returns_none_when_suprathreshold() -> None: - """A response that crosses spike threshold returns None (impedance undefined).""" +def test_returns_none_when_spikes_throughout() -> None: + """Spikes covering the entire window leave no spike-free segment → None. + + A single isolated spike is now tolerated via the spike-free sub-window + fallback (covered by ``test_recovers_profile_from_spike_free_segment``); to + exercise the genuine "no usable segment" path we tile spikes every 50 ms + so every excised guard window butts against the next. + """ time, current = _chirp_and_time(amplitude=1.0) - # A short +30 mV depolarization with a sharp upstroke reads as a spike. voltage = np.full_like(time, _BASELINE_MV) - spike_window = (time >= 1000.0) & (time < 1005.0) - voltage[spike_window] = 30.0 + # Spike every 50 ms (= _REFRACTORY_SAMPLES at this 1 kHz sample rate, so + # each event is detected) across the entire trace. + for centre_ms in np.arange(50.0, _DURATION_MS, 50.0): + spike_window = (time >= centre_ms) & (time < centre_ms + 5.0) + voltage[spike_window] = 30.0 result = analyze_impedance( time, voltage, current, 0.0, _window_end(time), _F_START, _F_END @@ -245,6 +253,72 @@ def test_returns_none_when_suprathreshold() -> None: assert result is None +def test_recovers_profile_from_spike_free_segment() -> None: + """A brief spike + a quiet resonator recovers via the sub-window fallback. + + The full chirp window contains a spike, but the spike-free remainder is long + enough (≫ the 50 ms minimum) to carry an FFT, so ``analyze_impedance`` + returns a profile and reports the analyzed sub-window duration. The + reported band is narrower than the requested band because the recovered + segment covers only the late portion of the linear frequency sweep. + """ + f0 = 8.0 + q0 = 2.0 + time, current = _chirp_and_time() + freqs = np.fft.rfftfreq(current.size, d=1.0 / _FS_HZ) + ratio = freqs / f0 + transfer = 1.0 / (1.0 - ratio**2 + 1j * ratio / q0) + voltage = _BASELINE_MV + _apply_transfer_function(current, transfer) + # Inject one synthetic spike near the start of the trace. + spike_window = (time >= 100.0) & (time < 105.0) + voltage[spike_window] = 30.0 + + result = analyze_impedance( + time, voltage, current, 0.0, _window_end(time), _F_START, _F_END + ) + + assert result is not None + assert result.analyzed_window_ms is not None + # The recovered window covers everything past the spike + guard padding, + # which on a 4000 ms trace is well over 3500 ms. + assert result.analyzed_window_ms > _DURATION_MS - 200.0 + # ``f_start`` / ``f_end`` report the requested band; the actually-covered + # band lives in ``frequencies``. + assert result.f_start == _F_START + assert result.f_end == _F_END + assert min(result.frequencies) >= _F_START + + +def test_impedance_unavailable_reason_messages() -> None: + """The public reason wrapper returns "" for OK inputs and a sentence otherwise.""" + time, current = _chirp_and_time() + clean_voltage = _BASELINE_MV + 10.0 * current + + # OK case → empty string. + assert ( + patch_sim.impedance_unavailable_reason( + time, clean_voltage, current, 0.0, _window_end(time), _F_START, _F_END + ) + == "" + ) + + # Window shorter than _MIN_WINDOW_MS → "shorter than the minimum". + short_reason = patch_sim.impedance_unavailable_reason( + time, clean_voltage, current, 0.0, 10.0, _F_START, _F_END + ) + assert "shorter than the minimum" in short_reason + + # Spikes throughout → "fired throughout" message. + busy_voltage = np.full_like(time, _BASELINE_MV) + for centre_ms in np.arange(50.0, _DURATION_MS, 50.0): + sp = (time >= centre_ms) & (time < centre_ms + 5.0) + busy_voltage[sp] = 30.0 + busy_reason = patch_sim.impedance_unavailable_reason( + time, busy_voltage, current, 0.0, _window_end(time), _F_START, _F_END + ) + assert "fired throughout" in busy_reason or "spike-free segment" in busy_reason + + def test_absolute_conversion_with_area() -> None: """Supplying area_cm2 fills the absolute MΩ spectrum; omitting it leaves it None.""" resistance = 40.0 # kΩ·cm² diff --git a/tests/unit/test_passive_properties.py b/tests/unit/test_passive_properties.py index f0e65e15..ba62aa99 100644 --- a/tests/unit/test_passive_properties.py +++ b/tests/unit/test_passive_properties.py @@ -9,10 +9,12 @@ import pytest from patch_sim.analysis.passive_properties import ( + _SPIKE_GUARD_POST_MS, analyze_passive_properties, density_to_absolute_c_m, density_to_absolute_r_in, is_subthreshold, + longest_subthreshold_run, ) # --------------------------------------------------------------------------- @@ -105,6 +107,64 @@ def test_is_subthreshold_flat_trace() -> None: assert is_subthreshold(time, voltage) is True +# --------------------------------------------------------------------------- +# longest_subthreshold_run +# --------------------------------------------------------------------------- + + +def _inject_spike(voltage: np.ndarray, time: np.ndarray, centre_ms: float) -> None: + """Inject a brief synthetic spike around ``centre_ms`` into a flat voltage trace. + + The spike is a 1 ms-wide +30 mV pulse — sharp enough to read as a + threshold-crossing event for ``analyze_aps`` at the default + ``dvdt_threshold=20`` mV/ms. + + Args: + voltage: Voltage trace to modify in place. + time: Time axis aligned with ``voltage``. + centre_ms: Centre of the spike in ms. + """ + mask = (time >= centre_ms) & (time < centre_ms + 1.0) + voltage[mask] = 30.0 + + +def test_longest_subthreshold_run_no_spikes_returns_whole_trace() -> None: + """A spike-free trace returns the full ``(0, N)`` index range.""" + time, voltage = _flat_trace(duration_ms=200.0) + run = longest_subthreshold_run(time, voltage) + assert run == (0, time.size) + + +def test_longest_subthreshold_run_one_spike_returns_trailing_gap() -> None: + """One early spike leaves a long trailing spike-free segment past the AHP guard.""" + time, voltage = _flat_trace(duration_ms=500.0) + _inject_spike(voltage, time, centre_ms=50.0) + + run = longest_subthreshold_run(time, voltage) + + assert run is not None + start_idx, stop_idx = run + # The recovered segment must start well after the spike's peak time + + # the post-guard window. + assert float(time[start_idx]) >= 50.0 + _SPIKE_GUARD_POST_MS + assert stop_idx == time.size + # And it should cover most of the remaining trace duration. + assert float(time[stop_idx - 1] - time[start_idx]) > 400.0 + + +def test_longest_subthreshold_run_spikes_throughout_returns_none() -> None: + """Tightly-packed spikes leave no spike-free span; the function returns ``None``.""" + time, voltage = _flat_trace(duration_ms=300.0) + # Refractory at 40 kHz is 1 ms, so 5 ms spacing is comfortably resolved. + # Spaced 5 ms apart, guard windows (~22 ms wide) fully overlap into one + # merged excision interval — and we start at 1 ms and extend past the end + # of the trace so the leading and trailing gaps both collapse to empty. + for centre in np.arange(1.0, 305.0, 5.0): + _inject_spike(voltage, time, centre_ms=float(centre)) + + assert longest_subthreshold_run(time, voltage) is None + + # --------------------------------------------------------------------------- # analyze_passive_properties — guard conditions # --------------------------------------------------------------------------- From 9c6c27ca161274a58f860c44e2d5395bacd65191 Mon Sep 17 00:00:00 2001 From: JCorson Date: Tue, 12 May 2026 07:45:35 -0500 Subject: [PATCH 8/9] fix(analysis): report impedance resonance only when f_R/Q/peak are all defined MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `analyze_impedance` used to set `resonance_frequency` for any interior `|Z|` maximum more than 1% above the low-frequency reference, while `quality_factor` (and `peak_impedance`, defined as "|Z| at f_R") required the −3 dB half-power crossings to be bracketable inside the band. The result was confusing rows like "f_R = 66 Hz / Q = — / Peak |Z| = —" for cells with no measurable resonance (squid axon, fast-spiking interneuron, the hyperpolarized pacemakers). Now a resonance is reported only when the interior peak both clears a 5% prominence margin (raised from 1%) and has a bracketable −3 dB width; `resonance_frequency`, `peak_impedance` and `quality_factor` populate or blank together. Of the shipped presets only Hippocampal CA1 (which has an Ih/HCN current) reports a resonance (f_R ≈ 9 Hz, Q ≈ 1.06); every other preset shows "—" for all three resonance fields, while the full |Z(f)| / phase plot is unchanged. No UI changes needed — the panel already maps None → "—" and only draws the f_R marker when one is reported. Tests: - tests/unit/test_impedance.py: new `test_shallow_interior_bump_reports_no_resonance` — a barely-underdamped 2nd-order low-pass has a genuine but sub-percent interior |Z| bump, so f_R / peak_impedance / Q are all None. - tests/integration/test_impedance_simulation.py: lengthen the test chirp to 1000 ms (1 Hz resolution — fine enough to bracket CA1's few-Hz Ih resonance, which a 500 ms / 2 Hz-resolution chirp put too close to the 1 Hz band edge); add `quality_factor is not None` to the CA1 resonance test. --- patch_sim/analysis/impedance.py | 53 +++++++++++-------- .../integration/test_impedance_simulation.py | 16 ++++-- tests/unit/test_impedance.py | 34 ++++++++++++ 3 files changed, 78 insertions(+), 25 deletions(-) diff --git a/patch_sim/analysis/impedance.py b/patch_sim/analysis/impedance.py index 56eee8bd..e15dcf2f 100644 --- a/patch_sim/analysis/impedance.py +++ b/patch_sim/analysis/impedance.py @@ -36,9 +36,11 @@ _MIN_BAND_BINS: int = 8 #: Fractional margin by which the ``|Z|`` peak must exceed the low-frequency -#: reference to count as a genuine resonance (rather than measurement ripple on -#: an otherwise monotone passive response). -_RESONANCE_REL_MARGIN: float = 0.01 +#: reference to be considered for a genuine resonance (rather than measurement +#: ripple on an otherwise monotone passive response). Clearing this margin is +#: necessary but not sufficient — a resonance is only reported when the peak +#: also has a bracketable -3 dB half-power width (see :func:`_quality_factor`). +_RESONANCE_REL_MARGIN: float = 0.05 #: A linear chirp's spectral power tapers toward the band edges. Contiguous #: leading/trailing FFT bins whose stimulus magnitude falls below this fraction @@ -58,11 +60,11 @@ class ImpedanceProfile: :class:`~patch_sim.analysis.passive_properties.PassiveProperties`. ``resonance_frequency``, ``peak_impedance`` and ``quality_factor`` are all - ``None`` when ``|Z(f)|`` has no interior peak that rises meaningfully above - the low-frequency reference — the hallmark of a passive low-pass cell with - no resonance. ``quality_factor`` may additionally be ``None`` when an - interior peak exists but its half-power width cannot be bracketed inside - the analysis band. + ``None`` together unless ``|Z(f)|`` has a genuine resonance — an interior + peak that both rises meaningfully above the low-frequency reference *and* + has a half-power (−3 dB) width that can be bracketed inside the analysis + band. A passive low-pass cell, or one whose mild interior bump has no + bracketable bandwidth, reports ``None`` for all three. Attributes: frequencies: Analysis-band frequency axis in Hz, ascending. @@ -70,19 +72,19 @@ class ImpedanceProfile: Î is in µA/cm²), aligned to ``frequencies``. phase: ``∠Z(f)`` in degrees, aligned to ``frequencies``; positive values mean the voltage leads the current. - resonance_frequency: Frequency of the ``|Z|`` peak in Hz, or ``None`` - when the maximum lies at a band edge or does not rise meaningfully - above the low-frequency edge (no genuine resonance). + resonance_frequency: Frequency of the ``|Z|`` resonance peak in Hz, or + ``None`` when there is no genuine resonance (the maximum lies at a + band edge, does not clear the prominence margin, or has no + bracketable −3 dB width). peak_impedance: ``|Z|`` at ``resonance_frequency`` in kΩ·cm², or - ``None`` when there is no interior peak. + ``None`` when there is no genuine resonance. quality_factor: Dimensionless ``Q = f_R / FWHM`` where FWHM is the full width of the ``|Z|`` peak at the half-power (−3 dB) level, or - ``None`` when there is no interior peak or the half-power crossings - cannot be bracketed inside the band. + ``None`` when there is no genuine resonance. magnitude_mohm: ``|Z(f)|`` converted to absolute MΩ when ``area_cm2`` is supplied; ``None`` otherwise. peak_impedance_mohm: Absolute peak impedance in MΩ when ``area_cm2`` is - supplied and an interior peak exists; ``None`` otherwise. + supplied and a genuine resonance exists; ``None`` otherwise. area_cm2: Membrane surface area in cm² used for the absolute conversion. ``None`` when only per-area values were requested. f_start: Lower edge of the requested analysis band in Hz. @@ -311,9 +313,12 @@ def analyze_impedance( both signals, takes the real FFT of each (``numpy.fft.rfft``), forms ``Z(f) = V̂(f) / Î(f)`` over the swept band ``[f_start, f_end]`` (trimming band-edge bins where the chirp has negligible spectral power), and reports - the magnitude and phase spectra plus the resonance frequency (argmax of - ``|Z|``, only when it is an interior maximum that rises above the - low-frequency reference) and quality factor. + the magnitude and phase spectra plus, when a genuine resonance is present, + the resonance frequency and quality factor. A genuine resonance requires + an interior ``|Z|`` maximum that both clears the prominence margin over the + low-frequency reference and has a half-power (−3 dB) width bracketable + inside the band; ``resonance_frequency``, ``peak_impedance`` and + ``quality_factor`` are otherwise all ``None``. Membrane impedance is a linear, small-signal quantity, so spike transients are excluded: when the windowed response contains spikes the analysis falls @@ -439,13 +444,19 @@ def analyze_impedance( # Reference the low-frequency end with a small average to be robust to # spectral leakage in the first bin. low_ref = float(np.mean(mag[: min(3, mag.size)])) - is_resonant = 0 < peak_idx < mag.size - 1 and float(mag[peak_idx]) > low_ref * ( + # A resonance is reported only when the interior |Z| peak both clears the + # prominence margin and has a measurable -3 dB width, so that f_R, the peak + # impedance, and Q always populate (or blank) together — a "resonance" with + # no bracketable half-power bandwidth is just measurement ripple on an + # otherwise monotone passive response. + peak_prominent = 0 < peak_idx < mag.size - 1 and float(mag[peak_idx]) > low_ref * ( 1.0 + _RESONANCE_REL_MARGIN ) - if is_resonant: + q = _quality_factor(fb, mag, peak_idx) if peak_prominent else None + if q is not None: resonance_frequency: float | None = float(fb[peak_idx]) peak_impedance: float | None = float(mag[peak_idx]) - quality_factor: float | None = _quality_factor(fb, mag, peak_idx) + quality_factor: float | None = q else: resonance_frequency = None peak_impedance = None diff --git a/tests/integration/test_impedance_simulation.py b/tests/integration/test_impedance_simulation.py index badc237e..97c9c267 100644 --- a/tests/integration/test_impedance_simulation.py +++ b/tests/integration/test_impedance_simulation.py @@ -18,7 +18,11 @@ ) from patch_sim.presets.protocols import build_protocol_from_preset -_DURATION_MS = 500.0 +# 1000 ms matches the tuned FREQUENCY_RESPONSE preset and gives 1 Hz frequency +# resolution — fine enough to bracket the −3 dB width of CA1's few-Hz Ih +# resonance (a 500 ms / 2 Hz-resolution chirp puts that peak too close to the +# 1 Hz band edge to bracket its low-side half-power crossing). +_DURATION_MS = 1000.0 _F_START = 1.0 _F_END = 100.0 @@ -74,11 +78,13 @@ def test_chirp_pipeline_produces_sane_impedance_profile(hh_model) -> None: def test_ca1_chirp_shows_resonance_and_high_impedance(hh_model) -> None: - """A hippocampal CA1 cell (with Ih) shows a low-frequency |Z| resonance. + """A hippocampal CA1 cell (with Ih) shows a genuine low-frequency |Z| resonance. The CA1 pyramidal preset has an Ih (HCN) conductance, so a subthreshold - chirp produces an interior peak in |Z(f)| in the few-Hz range, and its - low-frequency impedance is much larger than the low-resistance squid axon. + chirp produces an interior peak in |Z(f)| in the few-Hz range that is both + prominent and sharp enough to bracket a −3 dB width — so ``resonance_frequency``, + ``peak_impedance`` *and* ``quality_factor`` are all populated. Its + low-frequency impedance is also much larger than the low-resistance squid axon. """ ca1_profile, ca1_voltage = _run_chirp_profile(make_ca1_pyramidal(), amplitude=1.0) squid_profile, _ = _run_chirp_profile(hh_model, amplitude=1.0) @@ -88,6 +94,8 @@ def test_ca1_chirp_shows_resonance_and_high_impedance(hh_model) -> None: assert ca1_profile.resonance_frequency is not None assert 1.0 < ca1_profile.resonance_frequency < 30.0 assert ca1_profile.peak_impedance is not None + assert ca1_profile.quality_factor is not None + assert ca1_profile.quality_factor > 0.0 mag = np.asarray(ca1_profile.magnitude) assert ca1_profile.peak_impedance > mag[0] assert ca1_profile.peak_impedance > mag[-1] diff --git a/tests/unit/test_impedance.py b/tests/unit/test_impedance.py index 2a9ddbf0..e4c02419 100644 --- a/tests/unit/test_impedance.py +++ b/tests/unit/test_impedance.py @@ -167,6 +167,40 @@ def test_synthetic_resonance_detects_f_r() -> None: assert result.peak_impedance > mag[-1] +def test_shallow_interior_bump_reports_no_resonance() -> None: + """A barely-underdamped 2nd-order low-pass has an interior |Z| max but no resonance. + + With ``Q0`` just above the ``1/√2`` threshold the response has a genuine + interior maximum, but it rises only a fraction of a percent above the + low-frequency reference — well below the prominence margin — so + ``resonance_frequency``, ``peak_impedance`` and ``quality_factor`` are all + ``None`` together (a "resonance" with no measurable −3 dB width is just + ripple on a passive response). The magnitude/phase spectra are still + populated. + """ + f0 = 30.0 + q0 = 0.72 # just above 1/sqrt(2) ≈ 0.707 → a tiny interior peak exists + time, current = _chirp_and_time() + freqs = np.fft.rfftfreq(current.size, d=1.0 / _FS_HZ) + ratio = freqs / f0 + transfer = 1.0 / (1.0 - ratio**2 + 1j * ratio / q0) + voltage = _BASELINE_MV + _apply_transfer_function(current, transfer) + + result = analyze_impedance( + time, voltage, current, 0.0, _window_end(time), _F_START, _F_END + ) + + assert result is not None + mag = np.asarray(result.magnitude) + # There IS an interior maximum... + assert 0 < int(np.argmax(mag)) < mag.size - 1 + # ...but it is too shallow to be reported as a resonance. + assert result.resonance_frequency is None + assert result.peak_impedance is None + assert result.quality_factor is None + assert len(result.magnitude) == len(result.phase) > 0 + + def test_band_restriction() -> None: """The analysis band is clipped to [f_start, f_end].""" time, current = _chirp_and_time() From c32197975e14b7523964564195379abfc6ad150d Mon Sep 17 00:00:00 2001 From: JCorson Date: Tue, 12 May 2026 08:26:43 -0500 Subject: [PATCH 9/9] refactor(analysis): tighten impedance analysis per code review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - analyze_impedance now reports a band-maximum impedance (max_impedance / max_impedance_mohm) unconditionally, so the UI summary shows a useful scalar (≈ input impedance) even for passive cells with no resonance — previously the whole "Impedance Analysis" summary row was em-dashes for the squid axon, FSI, and the hyperpolarized pacemakers. peak_impedance (|Z| at f_R) still blanks together with f_R and Q. - Band-edge trimming keeps the longest *contiguous* run of FFT bins with usable stimulus power (was: the span from first-good to last-good bin), so an interior near-zero stimulus bin is dropped instead of dividing V̂/Î through it and spiking the estimate; _quality_factor also rejects a peak narrower than the FFT frequency resolution. New unit test covers the interior-dropout case. - Collapse the input-validation / windowing preamble into one _prepare_impedance_window helper so analyze_impedance and the public impedance_unavailable_reason wrapper each run analyze_aps once instead of 2–3 times per call. - Add patch_sim.constants.CHIRP_PROTOCOL and use it in place of the bare "Chirp" string literal in the sweep executor and the metrics sidebar. - US spelling: neighbor (impedance.py docstring), serialized (plotting/impedance.py), center_ms / center (test helpers). --- patch_sim/__init__.py | 2 + patch_sim/analysis/impedance.py | 240 +++++++++++--------- patch_sim/constants.py | 8 +- patch_sim_ui/components/metrics/__init__.py | 8 +- patch_sim_ui/plotting/impedance.py | 4 +- patch_sim_ui/state/_analysis_format.py | 20 +- patch_sim_ui/state/_sweep_executor.py | 20 +- tests/unit/test_impedance.py | 45 +++- tests/unit/test_passive_properties.py | 14 +- 9 files changed, 222 insertions(+), 139 deletions(-) diff --git a/patch_sim/__init__.py b/patch_sim/__init__.py index decc9250..dcecaf86 100644 --- a/patch_sim/__init__.py +++ b/patch_sim/__init__.py @@ -73,6 +73,7 @@ simulate_voltage_clamp_from_state, ) from .constants import ( + CHIRP_PROTOCOL, CURRENT_CLAMP, CURRENT_PROTOCOLS, DEFAULT_Q10, @@ -162,6 +163,7 @@ "run_membrane_test", "single_exp_decay", "single_exp_rise", + "CHIRP_PROTOCOL", "CURRENT_CLAMP", "CURRENT_PROTOCOLS", "DEFAULT_Q10", diff --git a/patch_sim/analysis/impedance.py b/patch_sim/analysis/impedance.py index e15dcf2f..bc787561 100644 --- a/patch_sim/analysis/impedance.py +++ b/patch_sim/analysis/impedance.py @@ -20,7 +20,6 @@ from patch_sim.analysis.passive_properties import ( density_to_absolute_r_in, - is_subthreshold, longest_subthreshold_run, ) @@ -42,29 +41,39 @@ #: also has a bracketable -3 dB half-power width (see :func:`_quality_factor`). _RESONANCE_REL_MARGIN: float = 0.05 -#: A linear chirp's spectral power tapers toward the band edges. Contiguous -#: leading/trailing FFT bins whose stimulus magnitude falls below this fraction -#: of the in-band maximum are trimmed before forming ``Z(f)``, so that dividing -#: by a near-zero stimulus does not inflate the band-edge impedance estimate. +#: A linear chirp's spectral power tapers toward the band edges (and can have +#: occasional interior dropouts from spectral leakage). After forming the +#: in-band spectrum the analysis keeps only the longest *contiguous* run of FFT +#: bins whose stimulus magnitude is at least this fraction of the in-band +#: maximum, so that dividing by a near-zero stimulus cannot inflate the +#: estimate at the band edges or spike it at an interior dropout. _MIN_STIM_FRAC: float = 0.05 +#: Windowed traces returned by :func:`_prepare_impedance_window` on success: +#: ``(t_win, v_win, i_win, analyzed_window_ms)`` where ``analyzed_window_ms`` is +#: ``None`` for the full chirp window or the sub-window duration in ms. +_PreparedWindow = tuple[np.ndarray, np.ndarray, np.ndarray, float | None] + @dataclasses.dataclass class ImpedanceProfile: """Membrane impedance profile from a chirp current-clamp run. - The per-area spectra (``magnitude``) are always populated. When - ``area_cm2`` is supplied to :func:`analyze_impedance`, the absolute - counterparts (``magnitude_mohm``, ``peak_impedance_mohm``) are filled in - too — analogous to the per-area / absolute split in + The per-area spectra (``magnitude``) and the band-maximum impedance + (``max_impedance``) are always populated. When ``area_cm2`` is supplied to + :func:`analyze_impedance`, the absolute counterparts (``magnitude_mohm``, + ``max_impedance_mohm``, ``peak_impedance_mohm``) are filled in too — + analogous to the per-area / absolute split in :class:`~patch_sim.analysis.passive_properties.PassiveProperties`. ``resonance_frequency``, ``peak_impedance`` and ``quality_factor`` are all ``None`` together unless ``|Z(f)|`` has a genuine resonance — an interior peak that both rises meaningfully above the low-frequency reference *and* has a half-power (−3 dB) width that can be bracketed inside the analysis - band. A passive low-pass cell, or one whose mild interior bump has no - bracketable bandwidth, reports ``None`` for all three. + band and is at least one FFT bin wide. A passive low-pass cell, or one + whose mild interior bump has no bracketable bandwidth, reports ``None`` for + all three; its ``max_impedance`` (the band-edge value, ≈ the input + impedance) is still available. Attributes: frequencies: Analysis-band frequency axis in Hz, ascending. @@ -72,17 +81,26 @@ class ImpedanceProfile: Î is in µA/cm²), aligned to ``frequencies``. phase: ``∠Z(f)`` in degrees, aligned to ``frequencies``; positive values mean the voltage leads the current. + max_impedance: The maximum of ``|Z(f)|`` over the analysis band in + kΩ·cm². For a cell with a genuine resonance this equals + ``peak_impedance``; for a passive low-pass cell it is the + lowest-frequency (≈ DC) value, i.e. roughly the input impedance. + Always populated. resonance_frequency: Frequency of the ``|Z|`` resonance peak in Hz, or ``None`` when there is no genuine resonance (the maximum lies at a band edge, does not clear the prominence margin, or has no bracketable −3 dB width). peak_impedance: ``|Z|`` at ``resonance_frequency`` in kΩ·cm², or - ``None`` when there is no genuine resonance. + ``None`` when there is no genuine resonance. When non-``None`` it + is identical to ``max_impedance`` (the resonance peak *is* the + band maximum). quality_factor: Dimensionless ``Q = f_R / FWHM`` where FWHM is the full width of the ``|Z|`` peak at the half-power (−3 dB) level, or ``None`` when there is no genuine resonance. magnitude_mohm: ``|Z(f)|`` converted to absolute MΩ when ``area_cm2`` is supplied; ``None`` otherwise. + max_impedance_mohm: ``max_impedance`` in MΩ when ``area_cm2`` is + supplied; ``None`` otherwise. peak_impedance_mohm: Absolute peak impedance in MΩ when ``area_cm2`` is supplied and a genuine resonance exists; ``None`` otherwise. area_cm2: Membrane surface area in cm² used for the absolute @@ -101,10 +119,12 @@ class ImpedanceProfile: frequencies: list[float] magnitude: list[float] phase: list[float] + max_impedance: float resonance_frequency: float | None peak_impedance: float | None quality_factor: float | None magnitude_mohm: list[float] | None = None + max_impedance_mohm: float | None = None peak_impedance_mohm: float | None = None area_cm2: float | None = None f_start: float = 0.0 @@ -124,7 +144,7 @@ def _half_power_crossing( Walks outward from ``peak_idx`` in ``direction`` (``-1`` toward lower frequencies, ``+1`` toward higher) to the first bin at or below ``half_power``, then linearly interpolates the crossing frequency between - that bin and its inward neighbour. + that bin and its inward neighbor. Args: frequencies: Analysis-band frequency axis in Hz, ascending. @@ -161,16 +181,20 @@ def _quality_factor( FWHM is the full width of the ``|Z|`` peak at the half-power (−3 dB) level, i.e. where ``|Z|`` equals ``peak / sqrt(2)``. The half-power crossings on each side of the peak are found by linear interpolation between bracketing - FFT bins. + FFT bins. A peak narrower than the FFT frequency resolution is not + resolved and is rejected (this also screens out lone single-bin spikes, + e.g. from an interior stimulus dropout). Args: - frequencies: Analysis-band frequency axis in Hz, ascending. + frequencies: Analysis-band frequency axis in Hz, ascending and + uniformly spaced. magnitude: ``|Z(f)|`` aligned to ``frequencies``. peak_idx: Index of the ``|Z|`` peak within the band. Returns: The quality factor, or ``None`` when either half-power crossing cannot - be bracketed inside the band or the bracketed width is non-positive. + be bracketed inside the band, the bracketed width is non-positive, or + the width is narrower than the frequency resolution. """ half_power = float(magnitude[peak_idx]) / np.sqrt(2.0) f_low = _half_power_crossing(frequencies, magnitude, peak_idx, half_power, -1) @@ -178,12 +202,36 @@ def _quality_factor( if f_low is None or f_high is None: return None fwhm = f_high - f_low - if fwhm <= 0.0: + df = float(frequencies[1] - frequencies[0]) + if fwhm < df: return None return float(frequencies[peak_idx]) / fwhm -def _impedance_unavailable_reason( +def _longest_usable_run(values: np.ndarray, threshold: float) -> tuple[int, int]: + """Return half-open bounds of the longest contiguous ``values >= threshold`` run. + + Args: + values: 1-D array (here: per-bin stimulus magnitudes). + threshold: Inclusive lower bound a value must meet to count as usable. + + Returns: + ``(lo, hi)`` such that ``values[lo:hi]`` is the longest contiguous + block of indices all meeting ``threshold``. ``(0, 0)`` when no value + meets the threshold. + """ + idx = np.flatnonzero(values >= threshold) + if idx.size == 0: + return 0, 0 + # Group boundaries: positions in `idx` where consecutive indices jump by >1. + breaks = np.flatnonzero(np.diff(idx) > 1) + 1 + starts = np.concatenate(([0], breaks)) + stops = np.concatenate((breaks, [idx.size])) + k = int(np.argmax(stops - starts)) + return int(idx[starts[k]]), int(idx[stops[k] - 1]) + 1 + + +def _prepare_impedance_window( time: np.ndarray, voltage: np.ndarray, injected_current: np.ndarray, @@ -191,15 +239,16 @@ def _impedance_unavailable_reason( stim_end_ms: float, f_start: float, f_end: float, -) -> str | None: - """Return a human-readable reason why impedance analysis can't proceed. +) -> tuple[str, None] | tuple[None, _PreparedWindow]: + """Validate the chirp inputs and extract the analyzable window. - Covers the cheap pre-FFT failure modes — shape / NaN / band sanity, a - too-short chirp window, and a suprathreshold response with no spike-free - fallback segment long enough for the FFT. The post-FFT failures (no - stimulus power in the band, too few in-band bins after trimming) are not - detected here and must be diagnosed by running :func:`analyze_impedance` - itself. + Performs the cheap pre-FFT checks (shape / NaN / band sanity, a too-short + chirp window, and a suprathreshold response that leaves no long-enough + spike-free segment), and on success returns the windowed traces — the full + chirp window when the response is subthreshold throughout, or the longest + contiguous spike-free sub-window otherwise. The post-FFT failures (no + stimulus power in the band, too few in-band bins after trimming) are *not* + detected here. Args: time: Time axis in ms. @@ -211,45 +260,52 @@ def _impedance_unavailable_reason( f_end: Ending frequency of the chirp sweep in Hz. Returns: - ``None`` when the preamble checks pass and the analysis can proceed; a - single-sentence reason string otherwise. + ``(reason, None)`` with a single-sentence failure reason, or + ``(None, (t_win, v_win, i_win, analyzed_window_ms))`` on success where + ``analyzed_window_ms`` is ``None`` for the full chirp window or the + sub-window duration (ms) when a spike-free segment was used. Exactly + one of the two tuple elements is non-``None``. """ t = np.asarray(time, dtype=float) v = np.asarray(voltage, dtype=float) i = np.asarray(injected_current, dtype=float) if v.shape != t.shape or i.shape != t.shape: - return "internal: trace length mismatch." + return "internal: trace length mismatch.", None if not (np.all(np.isfinite(v)) and np.all(np.isfinite(i))): - return "internal: trace contains non-finite values." + return "internal: trace contains non-finite values.", None if f_end <= f_start or f_start < 0.0: - return f"invalid frequency band [{f_start}, {f_end}] Hz." + return f"invalid frequency band [{f_start}, {f_end}] Hz.", None window_ms = stim_end_ms - stim_start_ms if window_ms < _MIN_WINDOW_MS: return ( f"the chirp window ({window_ms:.1f} ms) is shorter than the " f"minimum {_MIN_WINDOW_MS:.1f} ms needed for a meaningful spectrum." - ) + ), None mask = (t >= stim_start_ms) & (t < stim_end_ms) if int(mask.sum()) < 4: - return "too few samples in the chirp window." - t_win, v_win = t[mask], v[mask] - if is_subthreshold(t_win, v_win): - return None + return "too few samples in the chirp window.", None + t_win, v_win, i_win = t[mask], v[mask], i[mask] + + # longest_subthreshold_run returns (0, N) when the window has no spikes, so + # a single call covers both "subthreshold throughout" and "fall back to the + # widest spike-free segment" — no separate is_subthreshold check needed. run = longest_subthreshold_run(t_win, v_win) if run is None: return ( "the cell fired throughout the chirp — reduce the amplitude or " "apply a hyperpolarizing holding current." - ) + ), None lo, hi = run + if (lo, hi) == (0, t_win.size): + return None, (t_win, v_win, i_win, None) span = float(t_win[hi - 1] - t_win[lo]) if span < _MIN_WINDOW_MS: return ( f"the longest spike-free segment ({span:.1f} ms) is shorter than the " f"minimum {_MIN_WINDOW_MS:.1f} ms — reduce the amplitude or apply a " "hyperpolarizing holding current." - ) - return None + ), None + return None, (t_win[lo:hi], v_win[lo:hi], i_win[lo:hi], span) def impedance_unavailable_reason( @@ -261,14 +317,14 @@ def impedance_unavailable_reason( f_start: float, f_end: float, ) -> str: - """Public wrapper around :func:`_impedance_unavailable_reason`. + """Return a human-readable reason impedance analysis can't proceed, or ``""``. - Returns the empty string when the preamble checks pass (and the analysis - would proceed past the windowing / suprathreshold guards), so callers can - use the result directly as UI copy. Only the cheap pre-FFT failure modes - are reported — when this returns ``""`` and :func:`analyze_impedance` still - returns ``None``, the caller should fall back to a generic - "too little usable signal in the band" message. + Returns the empty string when the cheap pre-FFT checks pass (so the + analysis would proceed past the windowing / suprathreshold guards), letting + callers use the result directly as UI copy. Only those pre-FFT failure + modes are reported — when this returns ``""`` and :func:`analyze_impedance` + still returns ``None``, the caller should fall back to a generic "too little + usable signal in the band" message. Args: time: Time axis in ms. @@ -280,21 +336,13 @@ def impedance_unavailable_reason( f_end: Ending frequency of the chirp sweep in Hz. Returns: - The reason string from :func:`_impedance_unavailable_reason`, or ``""`` - when that function returns ``None``. + A single-sentence reason string, or ``""`` when the preamble checks + pass. """ - return ( - _impedance_unavailable_reason( - time, - voltage, - injected_current, - stim_start_ms, - stim_end_ms, - f_start, - f_end, - ) - or "" + reason, _ = _prepare_impedance_window( + time, voltage, injected_current, stim_start_ms, stim_end_ms, f_start, f_end ) + return reason or "" def analyze_impedance( @@ -311,14 +359,16 @@ def analyze_impedance( Extracts the chirp window from the traces, removes the DC component from both signals, takes the real FFT of each (``numpy.fft.rfft``), forms - ``Z(f) = V̂(f) / Î(f)`` over the swept band ``[f_start, f_end]`` (trimming - band-edge bins where the chirp has negligible spectral power), and reports - the magnitude and phase spectra plus, when a genuine resonance is present, - the resonance frequency and quality factor. A genuine resonance requires - an interior ``|Z|`` maximum that both clears the prominence margin over the - low-frequency reference and has a half-power (−3 dB) width bracketable - inside the band; ``resonance_frequency``, ``peak_impedance`` and - ``quality_factor`` are otherwise all ``None``. + ``Z(f) = V̂(f) / Î(f)`` over the swept band ``[f_start, f_end]`` (keeping + only the longest contiguous run of FFT bins where the chirp has usable + spectral power), and reports the magnitude and phase spectra, the + band-maximum impedance, plus — when a genuine resonance is present — the + resonance frequency and quality factor. A genuine resonance requires an + interior ``|Z|`` maximum that clears the prominence margin over the + low-frequency reference and has a half-power (−3 dB) width that is both + bracketable inside the band and at least one FFT bin wide; + ``resonance_frequency``, ``peak_impedance`` and ``quality_factor`` are + otherwise all ``None``. Membrane impedance is a linear, small-signal quantity, so spike transients are excluded: when the windowed response contains spikes the analysis falls @@ -330,7 +380,8 @@ def analyze_impedance( :func:`impedance_unavailable_reason` for the user-facing reason. Args: - time: Time axis in ms (uniformly sampled). + time: Time axis in ms, assumed uniformly sampled (the FFT frequency + axis is derived from the mean sample interval). voltage: Membrane voltage response in mV. injected_current: Injected chirp current in µA/cm². stim_start_ms: Start of the chirp stimulus window in ms. @@ -338,7 +389,7 @@ def analyze_impedance( f_start: Starting frequency of the chirp sweep in Hz. f_end: Ending frequency of the chirp sweep in Hz. area_cm2: Optional membrane surface area in cm². When supplied, the - absolute MΩ counterparts (``magnitude_mohm``, + absolute MΩ counterparts (``magnitude_mohm``, ``max_impedance_mohm``, ``peak_impedance_mohm``) are populated. Returns: @@ -350,37 +401,15 @@ def analyze_impedance( stimulus has no spectral power in the band, or fewer than :data:`_MIN_BAND_BINS` FFT bins remain after band-edge trimming. """ - reason = _impedance_unavailable_reason( - time, - voltage, - injected_current, - stim_start_ms, - stim_end_ms, - f_start, - f_end, + # _prepare_impedance_window returns exactly one of (reason, None) / + # (None, prepared), so a None `prepared` always carries a non-empty reason. + reason, prepared = _prepare_impedance_window( + time, voltage, injected_current, stim_start_ms, stim_end_ms, f_start, f_end ) - if reason is not None: + if prepared is None: logger.warning("analyze_impedance: %s", reason) return None - - t = np.asarray(time, dtype=float) - v = np.asarray(voltage, dtype=float) - i = np.asarray(injected_current, dtype=float) - mask = (t >= stim_start_ms) & (t < stim_end_ms) - t_win, v_win, i_win = t[mask], v[mask], i[mask] - - analyzed_window_ms: float | None - if is_subthreshold(t_win, v_win): - analyzed_window_ms = None - else: - # _impedance_unavailable_reason already verified that a long-enough - # spike-free run exists; fall back to it. - run = longest_subthreshold_run(t_win, v_win) - if run is None: # defensive — guarded above - return None - lo, hi = run - t_win, v_win, i_win = t_win[lo:hi], v_win[lo:hi], i_win[lo:hi] - analyzed_window_ms = float(t_win[-1] - t_win[0]) + t_win, v_win, i_win, analyzed_window_ms = prepared dt_ms = float(np.mean(np.diff(t_win))) if dt_ms <= 0.0: @@ -420,10 +449,10 @@ def analyze_impedance( if i_max <= 0.0: logger.warning("analyze_impedance: stimulus has no power in band; skipping.") return None - # Trim contiguous band-edge bins where the chirp has negligible power so - # that dividing by a near-zero stimulus does not inflate the estimate. - keep = np.flatnonzero(i_mag >= _MIN_STIM_FRAC * i_max) - lo, hi = int(keep[0]), int(keep[-1]) + 1 + # Keep only the longest contiguous run of bins with usable stimulus power so + # that dividing by a near-zero stimulus cannot inflate the band-edge + # estimate or spike it at an interior dropout. + lo, hi = _longest_usable_run(i_mag, _MIN_STIM_FRAC * i_max) if hi - lo < _MIN_BAND_BINS: logger.warning( "analyze_impedance: only %d FFT bins with usable stimulus power; " @@ -440,6 +469,7 @@ def analyze_impedance( mag = np.abs(z) phase_deg = np.degrees(np.angle(z)) + max_impedance = float(np.max(mag)) peak_idx = int(np.argmax(mag)) # Reference the low-frequency end with a small average to be robust to # spectral leakage in the first bin. @@ -454,8 +484,11 @@ def analyze_impedance( ) q = _quality_factor(fb, mag, peak_idx) if peak_prominent else None if q is not None: + # The resonance peak is the band maximum (peak_idx == argmax(mag)), so + # peak_impedance == max_impedance here; keep both for callers that want + # the resonance-specific value to be None when there is no resonance. resonance_frequency: float | None = float(fb[peak_idx]) - peak_impedance: float | None = float(mag[peak_idx]) + peak_impedance: float | None = max_impedance quality_factor: float | None = q else: resonance_frequency = None @@ -467,6 +500,7 @@ def analyze_impedance( magnitude_mohm: list[float] | None = [float(m) / area_cm2 / 1000.0 for m in mag] else: magnitude_mohm = None + max_impedance_mohm = density_to_absolute_r_in(max_impedance, area_cm2) peak_impedance_mohm = ( density_to_absolute_r_in(peak_impedance, area_cm2) if peak_impedance is not None @@ -477,10 +511,12 @@ def analyze_impedance( frequencies=fb.tolist(), magnitude=mag.tolist(), phase=phase_deg.tolist(), + max_impedance=max_impedance, resonance_frequency=resonance_frequency, peak_impedance=peak_impedance, quality_factor=quality_factor, magnitude_mohm=magnitude_mohm, + max_impedance_mohm=max_impedance_mohm, peak_impedance_mohm=peak_impedance_mohm, area_cm2=area_cm2, f_start=float(f_start), diff --git a/patch_sim/constants.py b/patch_sim/constants.py index 92abf6bd..28f2fd1b 100644 --- a/patch_sim/constants.py +++ b/patch_sim/constants.py @@ -47,13 +47,19 @@ CURRENT_CLAMP: str = "Current Clamp" VOLTAGE_CLAMP: str = "Voltage Clamp" +#: Protocol-type name for the chirp (linear frequency-sweep) current-clamp +#: protocol; a member of :data:`CURRENT_PROTOCOLS`. Referenced by name in the +#: analysis layer (which runs impedance analysis only for this protocol) and +#: the UI, so it gets its own constant rather than a bare string literal. +CHIRP_PROTOCOL: str = "Chirp" + #: Current clamp protocol type names. CURRENT_PROTOCOLS: list[str] = [ "Step", "Ramp", "Pulse Train", "Sinusoidal", - "Chirp", + CHIRP_PROTOCOL, "Noise", ] diff --git a/patch_sim_ui/components/metrics/__init__.py b/patch_sim_ui/components/metrics/__init__.py index 17977158..5b8972d3 100644 --- a/patch_sim_ui/components/metrics/__init__.py +++ b/patch_sim_ui/components/metrics/__init__.py @@ -21,7 +21,7 @@ import reflex as rx -from patch_sim.constants import CURRENT_CLAMP +from patch_sim.constants import CHIRP_PROTOCOL, CURRENT_CLAMP from patch_sim_ui.components.metrics.ap_panel import _ap_metrics_tab, _spike_row from patch_sim_ui.components.metrics.burst_panel import _burst_row from patch_sim_ui.components.metrics.calcium_panel import _ca_transient_row @@ -127,8 +127,7 @@ def _expanded_panel() -> rx.Component: rx.cond( ProtocolState.clamp_mode == CURRENT_CLAMP, rx.cond( - # "Chirp" is the protocol-type name from CURRENT_PROTOCOLS. - ProtocolState.protocol_type == "Chirp", + ProtocolState.protocol_type == CHIRP_PROTOCOL, rx.text("Impedance Analysis", size="4", weight="bold"), rx.text("AP Analysis", size="4", weight="bold"), ), @@ -153,8 +152,7 @@ def _expanded_panel() -> rx.Component: rx.cond( ProtocolState.clamp_mode == CURRENT_CLAMP, rx.cond( - # "Chirp" is the protocol-type name from CURRENT_PROTOCOLS. - ProtocolState.protocol_type == "Chirp", + ProtocolState.protocol_type == CHIRP_PROTOCOL, impedance_tab(), _ap_metrics_tab(), ), diff --git a/patch_sim_ui/plotting/impedance.py b/patch_sim_ui/plotting/impedance.py index 8cf6087b..ce620adb 100644 --- a/patch_sim_ui/plotting/impedance.py +++ b/patch_sim_ui/plotting/impedance.py @@ -6,7 +6,7 @@ def build_impedance_figure(imp_data: dict) -> go.Figure: - """Build a Plotly impedance-profile figure from serialised chirp results. + """Build a Plotly impedance-profile figure from serialized chirp results. Renders the impedance magnitude (left y-axis, teal) and phase (right y-axis, orange, dashed) plotted against frequency (Hz). When a resonance @@ -19,7 +19,7 @@ def build_impedance_figure(imp_data: dict) -> go.Figure: ``resonance_frequency`` (pre-formatted string, ``"—"`` when none). Returns: - A Plotly :class:`go.Figure` ready to be serialised and sent to the UI. + A Plotly :class:`go.Figure` ready to be serialized and sent to the UI. """ freqs = imp_data["frequencies"] magnitude = imp_data["magnitude"] diff --git a/patch_sim_ui/state/_analysis_format.py b/patch_sim_ui/state/_analysis_format.py index c8abf14e..192a5a35 100644 --- a/patch_sim_ui/state/_analysis_format.py +++ b/patch_sim_ui/state/_analysis_format.py @@ -493,13 +493,15 @@ def _compute_impedance_data( Returns: On success, a dict with keys ``frequencies``, ``magnitude``, ``phase`` - (lists), ``resonance_frequency``, ``quality_factor``, ``peak_impedance`` - (pre-formatted strings), and ``units`` — plus a ``caption`` string when - the profile was recovered from a spike-free sub-window rather than the - full chirp window. On failure, a dict with a single - ``unavailable_reason`` key carrying a human-readable sentence. Returns - an empty dict only when the sweep itself is too short to attempt - analysis. + (lists), ``resonance_frequency``, ``quality_factor`` (pre-formatted + strings, ``"—"`` when there is no resonance), ``peak_impedance`` (the + band-maximum ``|Z|`` as a pre-formatted string — always populated; for + a resonant cell this is ``|Z(f_R)|``, otherwise the ≈ DC / input + impedance), and ``units`` — plus a ``caption`` string when the profile + was recovered from a spike-free sub-window rather than the full chirp + window. On failure, a dict with a single ``unavailable_reason`` key + carrying a human-readable sentence. Returns an empty dict only when the + sweep itself is too short to attempt analysis. """ time_arr = np.asarray(sweep.time, dtype=float) if len(time_arr) < 2: @@ -541,11 +543,11 @@ def _compute_impedance_data( return {"unavailable_reason": reason} if profile.magnitude_mohm is not None: magnitude = profile.magnitude_mohm - peak = profile.peak_impedance_mohm + peak = profile.max_impedance_mohm units = "MΩ" else: magnitude = profile.magnitude - peak = profile.peak_impedance + peak = profile.max_impedance units = "kΩ·cm²" result: dict[str, Any] = { "frequencies": profile.frequencies, diff --git a/patch_sim_ui/state/_sweep_executor.py b/patch_sim_ui/state/_sweep_executor.py index 825f9c30..c95da34f 100644 --- a/patch_sim_ui/state/_sweep_executor.py +++ b/patch_sim_ui/state/_sweep_executor.py @@ -15,7 +15,7 @@ import patch_sim import patch_sim.channels -from patch_sim.constants import CURRENT_CLAMP, VOLTAGE_CLAMP +from patch_sim.constants import CHIRP_PROTOCOL, CURRENT_CLAMP, VOLTAGE_CLAMP from patch_sim_ui import constants from patch_sim_ui.api import traces from patch_sim_ui.plotting import TraceVisibility, build_figure @@ -100,12 +100,13 @@ def _compute_simulation( stimulus_step: Stimulus step size for analysis range. pre_stimulus_duration: Pre-stimulus duration (ms) for analysis windows. stimulus_duration: Stimulus duration (ms) for analysis windows. - protocol_type: Protocol-type name (e.g. ``"Step"``, ``"Chirp"``). Used - to decide which protocol-specific analyses to run. - start_frequency: Chirp sweep start frequency (Hz); only used when - ``protocol_type == "Chirp"``. - end_frequency: Chirp sweep end frequency (Hz); only used when - ``protocol_type == "Chirp"``. + protocol_type: Protocol-type name (e.g. ``"Step"``, + :data:`~patch_sim.constants.CHIRP_PROTOCOL`). Used to decide which + protocol-specific analyses to run. + start_frequency: Chirp sweep start frequency (Hz); only used for the + chirp protocol. + end_frequency: Chirp sweep end frequency (Hz); only used for the chirp + protocol. Returns: A :class:`_SimResult` containing sweeps, figure token, and all @@ -250,7 +251,8 @@ def _compute_simulation( burst_metrics, burst_summary = _compute_burst_data( ap_result, np.asarray(result["time"]) ) - # "Chirp" is the protocol-type name from constants.CURRENT_PROTOCOLS. + # The chirp protocol is always a single-sweep current clamp, so the + # impedance analysis only needs to be wired into this branch. impedance_data = ( _compute_impedance_data( sweep, @@ -260,7 +262,7 @@ def _compute_simulation( end_frequency, neuron.area_cm2, ) - if protocol_type == "Chirp" + if protocol_type == CHIRP_PROTOCOL else {} ) return _SimResult( diff --git a/tests/unit/test_impedance.py b/tests/unit/test_impedance.py index e4c02419..7f7717d4 100644 --- a/tests/unit/test_impedance.py +++ b/tests/unit/test_impedance.py @@ -100,6 +100,8 @@ def test_pure_resistor_flat_magnitude() -> None: assert result.resonance_frequency is None assert result.peak_impedance is None assert result.quality_factor is None + # max_impedance is always populated, even with no resonance. + assert result.max_impedance == pytest.approx(resistance, rel=1e-6) mag = np.asarray(result.magnitude) phase = np.asarray(result.phase) assert mag == pytest.approx(resistance, rel=1e-6) @@ -161,6 +163,8 @@ def test_synthetic_resonance_detects_f_r() -> None: assert result.peak_impedance is not None assert result.quality_factor is not None assert result.quality_factor > 0.3 + # When a resonance is reported, the resonance peak *is* the band maximum. + assert result.peak_impedance == pytest.approx(result.max_impedance) # The interior peak should clearly exceed the band edges. mag = np.asarray(result.magnitude) assert result.peak_impedance > mag[0] @@ -218,6 +222,37 @@ def test_band_restriction() -> None: assert freqs.max() <= 50.0 +def test_interior_stimulus_dropout_does_not_spike_the_estimate() -> None: + """A near-zero interior stimulus FFT bin is dropped, not divided through. + + The voltage carries the full chirp's content (response of a pure resistor + R), but the *stimulus* trace has one in-band FFT bin zeroed out — so a naive + ``V̂ / Î`` would blow up to infinity at that bin and steal the resonance + peak. Instead the analysis keeps only the longest contiguous run of usable + bins, so ``|Z|`` stays flat ≈ R, finite, with no spurious resonance. + """ + resistance = 25.0 # kΩ·cm² + time, current = _chirp_and_time() + voltage = _BASELINE_MV + resistance * current # response to the *full* chirp + # Knock out a single in-band FFT bin of the recorded stimulus only. + spectrum = np.fft.rfft(current) + band_freqs = np.fft.rfftfreq(current.size, d=1.0 / _FS_HZ) + in_band = np.flatnonzero((band_freqs >= _F_START) & (band_freqs <= _F_END)) + spectrum[int(np.median(in_band))] = 0.0 + notched_current = np.fft.irfft(spectrum, n=current.size) + + result = analyze_impedance( + time, voltage, notched_current, 0.0, _window_end(time), _F_START, _F_END + ) + + assert result is not None + mag = np.asarray(result.magnitude) + assert np.all(np.isfinite(mag)) + assert mag == pytest.approx(resistance, rel=1e-6) + assert result.resonance_frequency is None + assert result.max_impedance == pytest.approx(resistance, rel=1e-6) + + def test_returns_none_short_window() -> None: """A chirp window shorter than the minimum returns None.""" time, current = _chirp_and_time() @@ -276,8 +311,8 @@ def test_returns_none_when_spikes_throughout() -> None: voltage = np.full_like(time, _BASELINE_MV) # Spike every 50 ms (= _REFRACTORY_SAMPLES at this 1 kHz sample rate, so # each event is detected) across the entire trace. - for centre_ms in np.arange(50.0, _DURATION_MS, 50.0): - spike_window = (time >= centre_ms) & (time < centre_ms + 5.0) + for center_ms in np.arange(50.0, _DURATION_MS, 50.0): + spike_window = (time >= center_ms) & (time < center_ms + 5.0) voltage[spike_window] = 30.0 result = analyze_impedance( @@ -344,8 +379,8 @@ def test_impedance_unavailable_reason_messages() -> None: # Spikes throughout → "fired throughout" message. busy_voltage = np.full_like(time, _BASELINE_MV) - for centre_ms in np.arange(50.0, _DURATION_MS, 50.0): - sp = (time >= centre_ms) & (time < centre_ms + 5.0) + for center_ms in np.arange(50.0, _DURATION_MS, 50.0): + sp = (time >= center_ms) & (time < center_ms + 5.0) busy_voltage[sp] = 30.0 busy_reason = patch_sim.impedance_unavailable_reason( time, busy_voltage, current, 0.0, _window_end(time), _F_START, _F_END @@ -369,8 +404,10 @@ def test_absolute_conversion_with_area() -> None: assert with_area is not None and without_area is not None assert without_area.magnitude_mohm is None + assert without_area.max_impedance_mohm is None assert with_area.magnitude_mohm is not None assert len(with_area.magnitude_mohm) == len(with_area.magnitude) expected = resistance / area_cm2 / 1000.0 assert np.asarray(with_area.magnitude_mohm) == pytest.approx(expected, rel=1e-6) + assert with_area.max_impedance_mohm == pytest.approx(expected, rel=1e-6) assert with_area.area_cm2 == area_cm2 diff --git a/tests/unit/test_passive_properties.py b/tests/unit/test_passive_properties.py index ba62aa99..bb277555 100644 --- a/tests/unit/test_passive_properties.py +++ b/tests/unit/test_passive_properties.py @@ -112,8 +112,8 @@ def test_is_subthreshold_flat_trace() -> None: # --------------------------------------------------------------------------- -def _inject_spike(voltage: np.ndarray, time: np.ndarray, centre_ms: float) -> None: - """Inject a brief synthetic spike around ``centre_ms`` into a flat voltage trace. +def _inject_spike(voltage: np.ndarray, time: np.ndarray, center_ms: float) -> None: + """Inject a brief synthetic spike around ``center_ms`` into a flat voltage trace. The spike is a 1 ms-wide +30 mV pulse — sharp enough to read as a threshold-crossing event for ``analyze_aps`` at the default @@ -122,9 +122,9 @@ def _inject_spike(voltage: np.ndarray, time: np.ndarray, centre_ms: float) -> No Args: voltage: Voltage trace to modify in place. time: Time axis aligned with ``voltage``. - centre_ms: Centre of the spike in ms. + center_ms: Center of the spike in ms. """ - mask = (time >= centre_ms) & (time < centre_ms + 1.0) + mask = (time >= center_ms) & (time < center_ms + 1.0) voltage[mask] = 30.0 @@ -138,7 +138,7 @@ def test_longest_subthreshold_run_no_spikes_returns_whole_trace() -> None: def test_longest_subthreshold_run_one_spike_returns_trailing_gap() -> None: """One early spike leaves a long trailing spike-free segment past the AHP guard.""" time, voltage = _flat_trace(duration_ms=500.0) - _inject_spike(voltage, time, centre_ms=50.0) + _inject_spike(voltage, time, center_ms=50.0) run = longest_subthreshold_run(time, voltage) @@ -159,8 +159,8 @@ def test_longest_subthreshold_run_spikes_throughout_returns_none() -> None: # Spaced 5 ms apart, guard windows (~22 ms wide) fully overlap into one # merged excision interval — and we start at 1 ms and extend past the end # of the trace so the leading and trailing gaps both collapse to empty. - for centre in np.arange(1.0, 305.0, 5.0): - _inject_spike(voltage, time, centre_ms=float(centre)) + for center in np.arange(1.0, 305.0, 5.0): + _inject_spike(voltage, time, center_ms=float(center)) assert longest_subthreshold_run(time, voltage) is None