Skip to content

Clinical examples: SCI-like iEMG via motoneuron PIC up-regulation - #20

Merged
RaulSimpetru merged 64 commits into
mainfrom
example/pic-mechanistic-sci
Jun 28, 2026
Merged

Clinical examples: SCI-like iEMG via motoneuron PIC up-regulation#20
RaulSimpetru merged 64 commits into
mainfrom
example/pic-mechanistic-sci

Conversation

@RaulSimpetru

@RaulSimpetru RaulSimpetru commented Jun 28, 2026

Copy link
Copy Markdown
Member

Summary

Adds a new Clinical & Pathology example gallery (examples/04_clinical/) that reuses MyoGen's standard simulation pipeline and changes only the descending drive and motoneuron excitability (the persistent inward current, PIC) to reproduce clinically recognizable EMG signatures — demonstrated on spinal-cord-injury (SCI)-like intramuscular EMG.

The peripheral model (muscle geometry, electrode position, MU pool) is held fixed; the SCI discharge phenotypes emerge from up-regulating the NERLab motoneuron PIC (dendritic L-type Ca²⁺ via gamma, somatic persistent Na⁺ via nap_factor), not from hand-sculpting the drive.

What's included

  • pic_protocols.py — reusable NERLab PIC protocol helpers (single-cell pool builder, NaP scaling, ramp-hysteresis, after-discharge, pool spasm, iEMG synthesis, plateau-restricted ISI CV2).
  • 03_sci_iemg_mechanistic.py — the main mechanistic figure: three discharge phenotypes (voluntary modulation → loss of derecruitment → modulation-to-spasm) from one pool, varying only the PIC; single-cell PIC bistability panel; plateau-restricted ISI CV2.
  • 04_simulate_pic_spasticity.py — compact PIC-spasticity simulation.
  • README.rst + docs wiring — registers the gallery (docs/source/examples.rst, conf.py).
  • tests/test_pic_spasticity.py — 6 tests covering NaP scaling/caps, ramp hysteresis (control vs SCI), after-discharge, pool spasm, iEMG synthesis, and cyclic drive returning to zero at troughs.

Notes

  • Self-contained; reuses the existing simulation pipeline (no core API changes).
  • Underpins the SCI/PIC supplementary analysis in the MyoGen manuscript.
  • Planning/design docs used during development were intentionally kept out of this branch.

RaulSimpetru and others added 30 commits June 1, 2026 13:46
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
…MG example

Reproduces the SCI discharge phenotypes (self-sustained firing, loss of
derecruitment, modulation->spasm) from motoneuron PIC up-regulation with a fixed
physiological drive, instead of hand-sculpting the descending command. NERLab
model; reuses the example-10 clamp helpers and realistic iEMG noise.
…caled noise)

Rework the mechanistic SCI iEMG figure so its discharge statistics hold up to
scrutiny, addressing a peer-style critique of the earlier version.

run_pool (pic_protocols):
- Support the Powers2017 motoneuron model alongside NERLab. Powers2017's mAHP
  and dendritic Ca PIC put the self-sustained discharge at a physiological
  ~6-8 Hz (Gorassini 2004 ~5.2 Hz); NERLab's plateau floor (~12-16 Hz) is too
  high. Model-dependent NaP knob (lambda_factor vs scale_nap) and spike
  threshold (rest ~-71 vs ~0 mV).
- Inject an independent OU membrane-noise current per motoneuron whose amplitude
  scales with the descending drive and decays to an intrinsic floor when the
  drive withdraws (mn_noise, noise_floor). This makes voluntary firing irregular
  (ISI CV ~10%) while a drive-off self-sustained spasm, paced by the intrinsic
  PIC, is regular (~4%), reproducing the Gorassini spasm signature. Without it
  the PIC discharge is artifactually clock-like (CV <1%). Replaces the built-in
  Gfluctdv mechanism, whose Scop random source does not generate.
- Configurable descending-drive pool size (dd_n) and optional Gamma drive
  (dd_shape).

Analysis:
- population_rate: sliding-window estimator (smoother envelopes).
- population_cv: new sliding-window per-unit median ISI CV with inter-burst-gap
  exclusion, so the CV reflects within-burst regularity.

Example 03:
- Use Powers2017 with drive-scaled noise; recompute rate/CV per render.
- Honest docstrings/captions: the PIC is bistable (a gamma sweep shows no
  small-PIC self-sustained regime; ~22 nA is the dendritic plateau, larger than
  somatic voltage-clamp estimates that underread the distal dendrite); the
  second-model choice, the open-loop (no afferent) design, and the
  inhibition-dependent (non-self-terminating) spasm are stated as caveats.
The CV log axis was pinned to a fixed 0.2-20% window with empty space below
the data. Derive the shared range from the actual CV values (min/max across all
conditions, padded) so the axis spans only where the data live, making the
drive-off spasm collapse to ~3% read more clearly.
- Draw a dashed 'command off -> spasm' line at the moment the voluntary command
  stops (SPASM_ONSET_S, tied to the drive cutoff) on the modulation->spasm panel
  of all three figures, so the modulation->spasm transition is explicit.
- Scale the ISI CV y-axis per row (each condition's own min->max) instead of one
  shared range, so every panel's CV fills its own axis.
The rate envelope was scaled from an absolute 0, so conditions that never
silence (loss of derecruitment, spasm) sat compressed near the top. Scale each
panel's envelope from its OWN min (aligned with the EMG zero line) to its max
(the EMG peak), matching the per-row treatment of the CV axis, so the rate
modulation fills each panel. Right axis still reports true pps.
…e ISI CV

Replace the three separate figures with one composite that tells the whole
story: a single-cell PIC-mechanism band (Vm + dendritic Ca PIC current; brief
pulse latches the bistable plateau, inhibition switches it off) over a
drive / raster / iEMG / ISI-CV grid with columns = the three conditions under
the SAME descending drive, only the PIC state varied. The drive row makes the
"same input, different PIC" argument explicit; the CV row is anchored to the
Gorassini 2004 self-sustained-firing band (5.4 +/- 1.6 %).

pic_protocols: add single_cell_pic_mechanism() (Powers2017 near the bistable
threshold) and expose mahp_factor / dd_weight__uS knobs in run_pool.

Revert the mAHP discharge-rate tuning. Halving mAHP to lift the (low ~6 pps)
rate also destroyed the ISI statistics -- mAHP is what regularises firing, so
cutting it pushed the CV to ~26% and the spasm collapse to 14.7%. The diagnostic
spasticity signature is the low ISI CV, not the rate (Gorassini 2004: the
self-sustained rate even rises with added drive), so keep mAHP intact and accept
the low rate. Restores voluntary CV ~14%, loss-of-derecruitment ~10%, spasm
collapsing to ~3.5% (Gorassini band) when the drive withdraws.
- Combine the ISI CV onto the iEMG panel (purple, twin log axis, with the
  Gorassini 2004 band) and drop the separate CV row -- more compact, and the
  CV collapse reads directly against the continuing EMG. Frees vertical space
  for taller rasters.
- Add a y-axis buffer to the rasters so the first/last motor-unit marker rows
  are no longer clipped.
Combine the drive into the raster panel (grey, right axis), so the raster
directly shows firing tracking the command in the healthy column and persisting
through the drive troughs under loss-of-derecruitment/spasm. The composite is
now four rows: single-cell mechanism (Vm + Ca PIC), raster(+drive),
iEMG(+ISI CV).
…ron)

Switch from Powers2017 to NERLab so the figure uses the model the manuscript
describes, removing the "why two models?" problem. SCI spasticity is induced by
up-regulating the PIC with gamma (dendritic Ca) + NaP (somatic), which NERLab
needs together to self-sustain; healthy = gamma 0.5/NaPx1, SCI = gamma 1.3/NaPx5.

- single_cell_pic_mechanism(): NERLab-aware (gamma + nap_factor). The NERLab
  cell gives a clean latch + working inhibition off-switch and a physiological
  ~-9 nA dendritic Ca PIC (vs Powers2017's ~-22 nA).
- Retune the membrane noise for NERLab (mn_noise 3.0 -> 0.4): NERLab's voluntary
  ISI variability comes from the descending drive itself (~24-32%), so only a
  small noise floor is needed; the drive-off spasm CV then collapses to ~5.8%,
  inside the Gorassini 2004 band (5.4 +/- 1.6%).
- Honest caveats rewritten for NERLab: realistic voluntary rate (~23 Hz) but the
  self-sustained spasm (~11 Hz) is ~2x real spasms (NERLab's brief AHP), which
  cannot be lowered without killing self-sustain or contradicting the SCI
  literature (which reduces, not strengthens, the AHP).
…density

The opaque black-edged markers read as a regular grid; lower alpha and drop the
edges so overlapping spikes blend into density, matching the actual irregular
firing (ISI CV 24-32%).
facecolor none + colored edge (donut), s=5, alpha 0.9 -- cleaner than filled
dots and matches the original SCI example's marker style.
RaulSimpetru and others added 28 commits June 14, 2026 20:48
…ack edge)

Restore the original example's marker spec (s=9, facecolor=rank color,
edgecolors=black, linewidth=0.2, alpha=0.85) instead of the hollow rings, so the
composite matches the established SCI figure look.
Match the MyoGen ISI-CV figure style (examples/02_finetune/05): scienceplots
science+nature, white background, sans-serif, top/right spines off (re-enabled
per twin axis), thick axes, editable SVG/PDF fonts (svg.fonttype none,
pdf.fonttype 42), 300 dpi. Use font_scale 1.0 (not the comparison figure's 2)
and a bit more grid spacing since this composite has many more panels.
Importing myogen applies the seaborn default theme (lavender #EAEAF2 axes, no
spines) as an import side effect, overriding the science+nature style set at the
top -- hence the grey panels and missing left/bottom spines. Refactor the style
into apply_pub_style() and re-assert it right before the figure is built (also
force white facecolor + left/bottom spines explicitly).
…ight

Match the other MyoGen figures: sns.despine(trim=True, offset=2) on the main
axes (twin axes keep their right spine), and savefig with bbox_inches='tight'.
The comparison figure's axes.linewidth=2.0 was too thick for this dense
composite; 1.0 is proportional. (Spine thickness is linewidth, not dpi; the
SVG/PDF are vector and resolution-independent.)
Remove the right spine everywhere, including the drive/CV twin axes; keep their
right-hand tick labels so the secondary scales still read without a spine line.
Per-row split, each its own figure with its own figsize (like the original SCI
example): (1) single-cell mechanism (Vm + currents), (2) rasters + drive, (3)
iEMG + ISI CV. Shared despine_fig/save_fig helpers; saves
sci_mechanistic_{mechanism,raster,iemg}.{svg,pdf}.
…files

rasterized=True on the raster scatters and on the dense line traces (Vm, the Ca
PIC and NaP currents, the iEMG) -- they embed as 300-dpi bitmaps while axes,
text and the CV/drive lines stay vector. Renders fast in Affinity; mechanism PDF
9.1 -> 0.11 MB, iEMG 2.3 -> 0.18 MB. save_fig pins dpi=300 for the bitmaps.
Per request, leave the iEMG as crisp vector; only the raster scatters and the
single-cell Vm/current traces stay rasterized.
The project convention is a fixed Nature column width: 7.18 in (180 mm) double-
column, with deliberate per-figure heights. Set mechanism 7.18x3.0, raster
7.18x2.8, iEMG 7.18x2.5 (was arbitrary 11-13 in wide).
The SCI iEMG/drive figures used sns.despine(..., offset=5); match it (was 2).
The big filled donut scatter looked clumpy, esp. at 7.18 in. Switch to eventplot
(linelengths=0.8, linewidths=0.7), rainbow by first-spike order, rasterized --
the project's standard raster style, much cleaner for dense firing.
Apply the MyoGen eventplot raster style (linelengths=0.8, linewidths=0.7,
rainbow by first-spike order, rasterized) to 01_simulate_sci_iemg.py (was donut
scatter) and 02_pic_loss_of_derecruitment.py (was marker='|' scatter), for a
consistent raster look across the clinical examples.
- Remove all hardcoded fontsize= overrides (they were tuned for the old 13-in
  width and oversized at 7.18 in); scale fonts via sns paper context font_scale
  0.7 so they're proportional to the compact panels (final typography in post).
- Halve the raster figure height (7.18 x 2.8 -> 1.4).
- Drop the 'command off -> spasm' text annotation; keep the dashed vertical line.
The embedded raster images are lossless FlateDecode -- there is no encoding
'quality' to change, only resolution. The eventplot ticks are line art, so
rasterizing them just blurs thin lines; leave them vector (one LineCollection,
~same size, perfectly crisp). Only the ~220k-sample Vm/current traces stay
rasterized at 600 dpi.
…oire)

Too many spikes for vector ticks (overlapping thin lines clutter/artifact), so
rasterize them after all -- but at 1200 dpi (RASTER_DPI) so the dense thin ticks
don't alias/moire. Still lossless FlateDecode.
- Add a mean discharge-rate figure (sci_mechanistic_rate, shared y): voluntary
  ramps to ~0 at troughs, loss of derecruitment stays elevated, spasm drops at
  t=4 s but sustains ~10 pps (not zero) -- the key Gorassini comparison, was
  computed but never plotted.
- Share the iEMG y-axis across conditions so amplitudes are comparable.
- Clamp each CV panel's lower y-bound so the Gorassini band is visible in EVERY
  condition (healthy/loss stay above it, spasm drops into it), not just the spasm.
- Shorten the over-long single-cell title to a stub (caption goes in post).
Kept separate figures, rainbow rasters, and CV-on-iEMG per-row (the explicit
prior choices Codex didn't know about).
CV line 1.4 -> 0.9; replace the dense CV tick candidate set with a sparse
1-2-5-per-decade scheme so the small iEMG panels show ~2-3 ticks, not a dozen.
…cted ISI CV2

- 03_sci_iemg_mechanistic.py: PIC phenotypes (voluntary modulation, loss of
  derecruitment, modulation->spasm) from up-regulating gamma + NaP; per-cycle
  plateau-restricted CV2 markers (median across MUs), noise calibration
  (mn_noise=2.0, noise_floor=0.06).
- pic_protocols.py: NERLab-only shared helpers, consolidated imports, dead-code
  removal (Powers2017/mahp/lambda machinery).
- rename 01 -> 04_simulate_pic_spasticity.py; tidy 02; gallery-standard formatting.
- tests pass (6/6).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…comments

- Remove 01_simulate_sci_iemg.py (old drive-sculpting approach, superseded by the
  mechanistic 03) and 02_pic_loss_of_derecruitment.py (= condition 2 of 03).
- Keep the non-overlapping pair: 03 (pool-level, all phenotypes + iEMG + ISI CV)
  and 04 (single-cell PIC mechanism + spasm, covered by the tests).
- Trim 03 module docstring (48->23 lines) + verbose code-comment blocks, and 04
  docstring. ruff clean, py_compile OK, 6/6 tests pass.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
- Rename pic_protocols.py -> _pic_protocols.py so sphinx-gallery skips it as a
  helper (filename_pattern executes every .py and requires a title); update the
  two example imports and the test import.
- README: lead with motoneuron PIC excitability (the primary manipulation),
  descending drive second, with the peripheral model held fixed -- matching the
  manuscript's SCI/PIC framing. 6/6 tests pass.
…strate & cache provenance

- Truly silence the descending command after offset (zero its synaptic noise too),
  so the self-sustained discharge is paced by the intrinsic PIC alone (no residual
  drive). Verified: pool still self-sustains 40/40 MUs post-offset (CV2 median ~6%).
- Reseed identically before each pool run so descending-drive connectivity and OU
  membrane noise form the SAME substrate across conditions; only the PIC varies.
- Embed a hash of all sim params in the results cache filename so stale caches
  can't silently leak into a figure (no manual .pkl deletion needed).
- Single-cell mechanism uses the manuscript SCI regime (gamma=1.3, nap=5) -> peak
  L-type Ca ~10.1 nA (within the 5-15 nA range).
- Axis label CV -> CV2; mark example 02 (gamma 0.2->1.5) as illustrative, not the
  manuscript parameterization; exclude _pic_protocols.py from sphinx-gallery.
Code review (Codex) follow-ups:
- brief_command_drive: truly silent after offset (zero its noise too), matching
  the example's modulation_then_silence — post-offset firing is PIC-paced.
- per_mu_cv2: form a CV2 pair only when BOTH adjacent ISIs are within max_isi, so
  a dropped long gap never fabricates a spurious small-small pair. Medians are
  unchanged (voluntary ~20%, self-sustained ~6%).
- Cache provenance: results/MUAP/single-cell cache keys now include a hash of the
  _pic_protocols.py source AND the generated drive waveforms, so helper edits or
  drive changes can't silently reuse stale figures.
- Reword 'intrinsic PIC alone' -> 'PIC under a low membrane-noise floor'.

Tests (manuscript regression at the exact gamma 0.5->1.3 / nap 1->5 regime):
- Sustained single-cell L-type Ca ~10 nA, measured over the self-sustained
  plateau window (not the transient peak); 5-15 nA bound + 8-13 nA drift guard.
- Pool self-sustains vs baseline post-offset; seed before drive (deterministic),
  assert baseline stops (base_after <= 1).
…tion_mode

The IntramuscularEMG simulator returns monopolar per-contact channels; the
electrode's differential_matrix is a separate, optional tool that is not
auto-applied. Remove the misleading differentiation_mode='consecutive' arg
(it equals the default and had no effect on the output) and add a comment so
the illustrative single-trace iEMG isn't mistaken for a differential montage.
@RaulSimpetru
RaulSimpetru merged commit 5fcdd8a into main Jun 28, 2026
12 checks passed
@RaulSimpetru
RaulSimpetru deleted the example/pic-mechanistic-sci branch June 28, 2026 19:08
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant