Compute impedance profile from chirp current-clamp protocol#365
Merged
Conversation
…mp runs 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).
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.
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.
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).
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.
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.
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).
…l defined
`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.
- 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).
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
patch_sim.analyze_impedance/ImpedanceProfile(newpatch_sim/analysis/impedance.py): an FFT-based membrane impedanceZ(f) = V̂(f) / Î(f)computed over the swept band of a chirp current-clamp stimulus. Reports|Z(f)|(kΩ·cm², plus absolute MΩ when the neuron declares a surface area),∠Z(f)in degrees, the resonance frequencyf_R(interior peak of|Z|that rises above the low-frequency edge), and the quality factorQ = f_R / FWHM(−3 dB half-power width via linear interpolation). ReturnsNonefor degenerate inputs (mismatched lengths, non-finite values, invalid band, too-short window, too few in-band FFT bins).numpy.fft.rfft; only DC removal is applied before transforming — a linear chirp already starts at zero phase and has an approximately flat band spectrum, so a Hann window is counter-productive here (it crushes the low/high-frequency content and inflates the band-edge estimate). This is a deliberate deviation from the "Hann window" suggestion in the issue.f_R/Q/ peak-|Z|summary plus a|Z(f)|(left axis) and∠Z(f)in degrees (right axis) plot. The profile is computed in the single-sweep current-clamp branch of the sweep executor and threaded through_SimResult→AnalysisState. Other protocols / clamp modes keep their existing AP-metrics / I-V views.Nonecases, area conversion); 2 integration tests against real simulations (structural sanity on the squid axon; thalamic relay has far higher low-frequency|Z|than the squid); 2 e2e tests (Frequency Response preset populatesimpedance_data; a non-chirp protocol leaves it empty). The e2erun_flowhelper now also threadsprotocol_type/ chirp frequencies into_compute_simulation.Closes #185.
Test plan
f_R/Q/ peak-|Z|summary plus the|Z(f)|+∠Z(f)plot.