[Track] multiscale-coupling: SimNIBS E-field → NEURON neuron response#125
Open
idossha wants to merge 16 commits into
Open
[Track] multiscale-coupling: SimNIBS E-field → NEURON neuron response#125idossha wants to merge 16 commits into
idossha wants to merge 16 commits into
Conversation
Add optional tit.microscale module that couples SimNIBS TI/mTI E-fields to NEURON multicompartment neurons under the quasi-static approximation. - config.py: MicroscaleConfig, NeuronModelSpec with carrier/dt validation - field_sampler.py: mm<->um conversion, uniform + path quasipotentials, SimNIBS-native field sampling (NEURON-free, pure NumPy core) - NEURON added to pyproject + both SimNIBS Dockerfiles (lazy-imported) - 16 unit tests for the pure-math core
- models.py: authored ball-and-stick cell (built-in hh + extracellular, no vendored assets); Cell wrapper with 3D segment geometry; model registry. Realistic Blue Brain/Aberra cells are CC-BY-NC-SA -> not vendored; users register their own via register_model(). - coupling.py: two-carrier extracellular drive (envelope emerges from active channels), uniform-field quasipotentials, simulate_response/polarization_map/ find_threshold (NEURON lazy-imported); pure waveform + spike math. - field_sampler.py: rotation_align + place_morphology orientation geometry. - metrics.py: targets CSV + response/polarization npz writers. - __main__.py: response|threshold mode dispatch. - paths: pm.microscale()/microscale_sim(); config_io: MicroscaleConfig _type. - 15 more unit tests (31 total); full suite green (3 pre-existing flex fails).
- gui/extensions/microscale.py: Microscale extension (subject/sim/model select, response|threshold mode, carrier freqs, target list) -> simnibs_python -m tit.microscale in a background thread, mirroring the Source extension. - README.md: physics (Ve = -E.l_c, two-carrier drive), mm<->um/units, licensing rationale (realistic cells are CC-BY-NC-SA, not vendored), usage, and Wang 2022 validation-target thresholds.
- find_threshold: fix early-return reporting 0.0 (now reports a positive floor when the cell already fires at the smallest probe); switch to geometric bisection (uniform relative precision) matching the docstring. - simulate_response: remove Vector.play handles and zero e_extracellular in a finally block so a run can't perturb later runs (matters for find_threshold's repeated calls); exclude an initial settling window from the spike count to drop the init transient. - polarization_map: measure field-OFF baseline on the SAME cell (rest first, then field), not a second coexisting cell in NEURON's global state; zero e_extracellular in finally. - neuron made an OPTIONAL dependency (project.optional-dependencies.microscale) so pip install . doesn't require a neuron wheel; pinned >=8.2,<8.3 (pyproject + both Dockerfiles) to keep numpy at 1.26.x. - 3 more tests for find_threshold logic (37 total).
…build Verified end-to-end inside a local idossha/simnibs:v2.5.0 image (NEURON 8.2.7, numpy 1.26.4 pinned). The coupling machinery checks out: Vector.play drives e_extracellular (+/-200 mV swing), polarization map is orientation-sensitive (perpendicular field -> 0), and the cell fires. Fix: under a (quasi-)uniform field the lowest-threshold elements are axon terminals, not the soma (Aberra 2020). simulate_response now counts activation at the spike-initiation site (a spike in any compartment) instead of only at soma(0.5); the somatic trace is still returned. After the fix the TI sweep is monotonic (0/0/119/143 spikes) and find_threshold converges (amp~51). Adds container/blueprint/Dockerfile.simnibs.v250local: layers neuron + the working-tree tit onto idossha/simnibs:v2.4.0 for local testing (NOT pushed).
…naming Found by running the real CLI against a real simulation (sub-101/L_Insula_scalar): - _load_pair_meshes looked in high_Frequency/ but the simulator organizes the pair meshes under high_Frequency/mesh/ -> fixed the path. - _build_ball_stick read sec.L for the second pt3dadd AFTER pt3dclear() had zeroed it, collapsing the apical dendrite and axon to ~0 length (the whole morphology shrank to the soma, under-driving the coupling ~60x). Use explicit length constants; verified max|ve| now 89.4 mV vs 90 mV hand calc for a 150 V/m field (was 1.5 mV). - response vs threshold modes wrote the same _response.npz; threshold now writes _threshold.npz so a threshold run can't clobber response spikes.
New tit/microscale/viz.py renders four artifacts per target, all headless matplotlib (no pyvista/VTK; GIF via pillow since the image has no ffmpeg): - morphology (3D, colored by part), - neuron embedded in a cropped cortical-surface patch (colored by TI_max), - TI E-field 3D vectors around the target with the neuron for scale, - animated GIF clip: per-compartment Vm + the oscillating E-field drive over time (envelope-tracking is clearly visible). - simulate_response gains return_traces=True (v_all, ve_t, world_um) for the clip. - Cell.section_spans() groups segments per section for rendering. - new 'viz' CLI mode; render_target() orchestrates everything from a real sim. - render_target globs for *_TI_central.msh (montage name can differ from the sim folder, e.g. L_Insula vs L_Insula_scalar) and skips the cortex panel if absent. - pure helpers (section_polylines, crop_surface_patch, grid_around) unit-tested (40 tests total). Generated + verified on sub-101/L_Insula_scalar.
Raises the visualization to the Shirinpour 2021 / Wang 2023 / Mirzakhalili 2020 standard (assessed from the papers in the local academic library). - new tit/microscale/morphology.py: MorphologySpec/SectionSpec; pyramidal_l5() generates a branched L5 cell (soma, basal dendrites, apical trunk+tuft, AIS, myelinated axon with nodes of Ranvier; ~50 sections, region-tagged, deterministic); load_swc() parses standard SWC reconstructions. - models.py: build_from_spec() builds NEURON cells from a spec (nseg by length, hh+extracellular, region tags carried to section_spans); register l5_pyramidal as the DEFAULT model; keep ball_stick; load_swc_cell() for user reconstructions. - viz.py: Shirinpour region palette (red soma / orange basal / green apical / blue tuft / purple axon / highlighted AIS+nodes); diameter-scaled neurites; 100um scale bar; plot_morphology gains a color-by-value mode (Vm / quasipotential with a diverging colormap); new *_quasipotential.png panel (Fig 2F style); cortex panel tightened (4mm patch, neuron drawn on top + soma marker). - 4 new tests (morphology structure + SWC parse); 39 total. Generated + verified on sub-101/L_Insula_scalar.
Two features (delegated to and built by parallel sub-agents, integrated here). POPULATION (tit/microscale/population.py, the standard approach for subthreshold polarization per Aberra 2018/2020, Seo & Jun 2017, Shirinpour 2021): - analytic_polarization_map(): first-order ΔVm = coupling x E_normal for every cluster vertex (coupling 0.27 mV/(V/m), Radman 2009 / Bikson 2004). - run_population(): unconnected population = clones x azimuthal rotations on a representative subsample of cluster vertices; characterizes the distribution around the analytic central estimate. No synaptic connectivity (quasi-uniform justification in the module docstring). - select_cluster/azimuths/load_cluster_surface helpers; PopulationConfig; place_morphology gains azimuth_deg (rotate about the cortical normal); metrics population npz/csv writers; 'population' CLI mode. - Verified on sub-101: analytic ΔVm over 246,384 vertices + NEURON subsample. ROTATING TI FIELD (viz.py): the instantaneous resultant E1*sin(w1 t)+E2*sin(w2 t) rotates over the beat cycle; the old clip drew a fixed-direction arrow (wrong). - instantaneous_field() pure helper; animate_response() now draws the true rotating arrow; new plot_field_hodograph() proves the rotation (Lissajous in the E1-E2 plane); simulate_response(return_traces) returns e1_vec/e2_vec. - new *_hodograph.png artifact in render_target. Config/exports: PopulationConfig registered in config_io; __init__ exports the new public API (run_population, analytic_polarization_map, instantaneous_field, plot_field_hodograph, ...). 56 microscale tests, full suite green.
Adds a populated-cortex render matching the Aberra et al. 2018/2020 figure: many neurons placed along a localized cortical cross-section, oriented to each site's cortical normal, colored by neurite TYPE (axon red / apical blue / basal green). Confirmed by research (no-license on the Aberra repo/morphologies -> procedural L5 default + load_swc for user cells; Blender/BlenderSpike not worth it -> matplotlib; LineCollection-per-type + 2D projection per Mirzakhalili MIT techniques). - viz.py: NEURITE_COLORS_ABERRA + COLOR_SCHEMES; plot_population_in_cortex() (2D projection, LineCollection per neurite type, gyral outline, scale bar, legend); render_population_cortex() localizes to a patch around the field focus, slabs the thinnest in-plane axis, places cells along the gyral section. - population.py: sample_cortical_strip() (robust median-centered slab) and place_spec_world() (pure-geometry placement; no NEURON). run_population now also emits *_population_cortex.png. - 2 new pure tests (strip selection, spec placement); 58 microscale tests. Verified on sub-101: 60 L5 cells along the gyral cross-section, neurite-typed.
The populated-gyrus cells previously floated over faint dots with no anatomical context. Now they are embedded in the actual cortical-surface cross-section: - load_cluster_triangles() returns the central-surface connectivity. - render_population_cortex crops the surface (vertices+normals+TI_normal+tris in one consistent indexing) around the field focus, filters ribbon triangles to the slab, and draws the cortical ribbon as a translucent PolyCollection colored by TI_normal, with the neurons placed on it. - plot_population_in_cortex gains surface_tris/surface_scalar (projected surface mesh) and a region_label box (space + patch + coloring context). - Verified on sub-101: 60 L5 cells embedded in the gyral ribbon (subject space), ribbon colored by TI_normal (field focus visible). 58 tests green.
…ortical) Generalizes the populated-region figure to arbitrary regions in any space, so cells are embedded in a named anatomical context (no longer floating): - population.select_region(): atlas label (DK40/a2009s/HCP_MMP1, nearest-neighbour transfer to the remeshed TI surface), MNI/subject sphere on the cortical surface, or a subcortical MNI/subject sphere on the TI VOLUME mesh (radial orientation). MNI->subject via simnibs.mni2subject_coords; volume mesh globbed (montage-vs-folder naming). scipy cKDTree for the atlas snap. - viz.render_population_region(): renders a region dict; surface regions as the embedded gyral ribbon (colored by TI_normal), subcortical as the nucleus point cloud with radial cells; region-context label. - Verified on sub-101: left insula (DK40), MNI sphere r=10mm @ L-insula [-38,6,2], and left thalamus (MNI sphere r=8mm @ [-10,-19,8], subcortical volume). - Exports select_region/render_population_region; 58 tests green.
Scope per request: L5 pyramidal only; regions = atlas label / sphere / binary mask, all GM-restricted, subject space (fsaverage scaffolded, next iteration). - config.RegionSpec: clean typed region (kind atlas|sphere|mask, space, validated). - population.select_region rewritten: GM cortical surface only (dropped the subcortical-volume path); atlas via nearest-neighbour annot transfer, sphere (MNI/subject center), binary-NIfTI mask sampled at vertices (GM inclusion). - viz.plot_population_3d: clean 3D render — smooth lit cortical patch colored by TI_normal (magma) with L5 cells embedded as Line3DCollections by neurite type, oriented to the cortical normal; colorbar, legend, scale bar, high DPI (220). Avoids the 2D-projection self-overlap of a folded gyrus. - render_population_region: RegionSpec-based, windows a patch around the field focus so the ~1 mm cells are clearly visible, 3D render. - RegionSpec validation tests (59 total). Verified on sub-101: insula (DK40), MNI sphere r=10mm, and binary mask — all clean/professional.
- _atlas_vertex_mask: exact label match first; reject ambiguous substring matches (e.g. 'temporal' -> 4 labels) instead of silently unioning them. - _mask_vertex_mask: document the subject-space assumption. - RegionSpec: validate sphere center is a 3-tuple and snap_mm > 0. - plot_population_3d: drop dead code; render_population_region: drop the unused thickness_mm param. - render_population_cortex now delegates to render_population_region (3D), via a sphere at the field focus — retires the duplicate 2D entry point from the auto figure (plot_population_in_cortex stays as the 2D option). - 59 tests green; insula/sphere/mask all regenerate clean.
Update microscale coupling module (config, models, coupling, metrics, population, viz, GUI extension), add population/meso run scripts, and add microscale-coupling wiki page.
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
Adds an optional, research-grade
tit.microscalemodule that couples themacroscale TI/mTI electric field the simulator already computes to NEURON
multicompartment neuron models — letting researchers see how neurons in an ROI
respond to the exposure, accounting for field intensity and orientation,
not just magnitude. This is the single-neuron (NEURON) scale; whole-brain
neural-mass / TVB coupling is deliberately out of scope (it discards orientation).
Motivated by a feasibility question (TVB/NEURON-style multiscale coupling). The
approach follows the validated literature recipe (Aberra 2020, Wang 2022,
Shirinpour 2021) under the quasi-static approximation; the field→compartment
formula
Ve(c) = −E·l_cwith two-carrier superposition matches Wang et al. 2022verbatim.
What's included
config.py—MicroscaleConfig,NeuronModelSpecwith carrier/dt(Nyquist) validation.
field_sampler.py— NEURON-free core: mm↔µm conversion, uniform + pathquasipotentials, orientation geometry (
rotation_align,place_morphology),SimNIBS-native field sampling.
models.py— authored ball-and-stick cortical cell (built-inhh+extracellular, no vendored assets);Cellwrapper + model registry.coupling.py— two-carrier extracellular drive (envelope emerges fromactive channels),
simulate_response/polarization_map/find_threshold(NEURON lazy-imported); pure waveform + spike-count math.
metrics.py— targets CSV + response/polarization.npzwriters.__main__.py—response | thresholdmode dispatch(
simnibs_python -m tit.microscale config.json).Microscaleextension (tit/gui/extensions/microscale.py).pm.microscale()/microscale_sim();MicroscaleConfig_typein
config_io;neuronadded topyproject+ both SimNIBS Dockerfiles(lazy-imported).
tit/microscale/README.md(physics, units, licensing, validationtargets).
Licensing note
The realistic Blue Brain / Aberra cortical morphologies are CC-BY-NC-SA
(non-commercial, share-alike) and are not vendored. The shipped default is an
authored, license-free ball-and-stick cell; users register their own realistic
cells via
register_model().Testing
orientation geometry, carrier waveform, spike counting, config validation,
config_io round-trip, PathManager wiring, metrics IO).
test_opt_flex.py::TestRunFlexSearchfailures are pre-existing on
main, unrelated to this change).blackclean.Not verified locally (requires the container)
The actual NEURON cable solve runs inside the SimNIBS Linux container with real
subject meshes. This PR wires it up and documents it, but the in-container
end-to-end run (and the
neuronimage rebuild) must be executed on the image —see Phase 1 task notes in the track. NEURON is lazy-imported, so the rest of the
package (and the test suite) works without it.
Follow-ups
already supports it; only the integration step changes).