Cortical stimulation · methods

Experiments worth trying

The reference run uses a 40 µm disk, a 10 µA cathodic leading phase lasting 200 µs, a 50 µs gap and an equal anodic recovery phase. It solves a seeded sample of 512 cells in a 600 µm cube.

1. Size versus selectivity. Save a single-contact run as a comparison. Change only contact diameter and rerun. Inspect delivered current, charge density, membrane response and the sampled recruitment curve together. Smaller contacts do not automatically produce more selective neural activation.
2. Move the return. Compare the remote-return single contact with a local bipolar pair and the surround-return preset at the same total leading cathodic current. Return contacts can also excite neural processes. Field confinement and neural selectivity are different measurements.
3. Steer a fixed current. Choose Current steering. Keep amplitude and contact spacing fixed; move the share sent to E2. Compare the resulting recruitment pattern using the same anatomy seed.
4. Find the first excited segment. Select a spiking cell, use Inspect cell, and compare soma distance with neurite distance. The wire marker identifies the compartment with the first upward 0 mV crossing. Axonal recruitment can precede the somatic trace, consistent with the mechanism explored by Histed et al.
5. Separate material from geometry. Switch between PtIr and iridium-oxide circuit examples. When compliance is not reached, the delivered current and tissue field should agree. When compliance clips a pulse, its actual amplitude, charge and neural response change.
6. Check a conclusion. Double the solved sample, reduce time step from 5 to 2.5 µs, and reduce segment length from 20 to 10 µm. A result that changes substantially under refinement is not ready to interpret quantitatively.

1. The extracellular electric field

The medium is an unbounded, homogeneous, resistive volume conductor. Conductivity is a positive diagonal tensor Σ = diag(σ_x, σ_y, σ_z); σ_x = σ_y, and the depth/tangential control sets σ_z/σ_x. The default isotropic value, 0.276 S/m, is a modeling choice used in Aberra et al., not a measurement of the individual simulated tissue.

The anisotropic Green function follows the volume-conductor formulation discussed by Hindriks et al. Equal-area quadrature nodes distribute each contact’s total current across its finite source surface. Disks and annuli use area-uniform radial sampling; rectangles use a regular area grid; spheres use approximately equal-area surface points. Sources superpose linearly, so signed electrode weights implement steering and local returns.

Flat contacts are ideal thin sheets immersed in tissue with both faces exposed and an insulated rim. Their electrochemical area is twice their projected area. The field uses the collapsed source sheet; the 3D pad thickness is a rendering aid. A sphere uses its full 4πr² surface. These are prescribed uniform current distributions, not equipotential metal surfaces. Edge crowding, insulating substrates, shanks and redistribution of surface current are absent. The more complete finite-element boundary treatment in Joucla & Yvert is the appropriate comparison.

A small quadrature core regularizes each patch: its radius is half the radius of an equal-area circular patch. This recovers that patch’s central mean inverse distance and decreases as quadrature is refined. It is a numerical integration device, not a biological distance cutoff. Contact-intersecting cells are excluded, and field pixels within 2 µm of a contact are masked.

Geometry is stored in µm, input current in µA, potential in mV, electric field in V/m and current density in A/m². All conversions occur explicitly. The displayed section is y = 0 at the maximum delivered leading-phase drive; its field is static during recruitment playback. The playback animates neural recruitment, not electromagnetic propagation. A remote ideal return at infinity supplies any nonzero net contact current.

2. Material, polarization and delivered current

Each electrode has an effective parallel capacitance and charge-transfer resistance. The tissue access matrix is computed by averaging the same field kernel across receiving contact surfaces. This keeps its full-space geometry consistent with the field calculation.

The PtIr example uses 30 µF/cm² and 100 Ω·cm². Reported PtIr DBS interface ranges include approximately 12–47 µF/cm² and 40–290 Ω·cm² under particular pulse conditions (Wei & Grill; Howell et al.). Applying one equivalent pair to a microcontact remains an assumption. Iridium oxide (2,000 µF/cm²; 1,000 Ω·cm²) and TiN (1,000 µF/cm²; 500 Ω·cm²) are explicitly illustrative circuit examples. They are not fitted measurements or interchangeable coating specifications.

The command is biphasic: leading current I for duration t_p, a gap, then −I/k for duration kt_p. Integrating partial time intervals avoids charge errors when an edge lies between solver steps. Frequency is limited so phases cannot overlap. If a contact reaches voltage compliance, every channel is scaled by a common factor for that interval. This preserves the chosen current ratios and exact local return balance. It models a coordinated controller, not a specific stimulator.

Different clipping during the two phases can leave nonzero delivered charge. The contact table reports this imbalance rather than silently forcing it to zero. Between pulses, each interface decays through its internal charge-transfer branch; no active discharge circuit is added. The ideal-source option bypasses interface and compliance constraints.

These contextual measurements are not thresholds in the simulation. Charge-storage capacity, reversible charge-injection capacity and tissue injury are different quantities. The model does not predict electrolysis, electrode degradation or tissue damage.

3. Active, branched membrane dynamics

Every solved neuron includes an ellipsoidal soma, an apical trunk and branches, three basal processes, a 40 µm axon initial segment, a descending axon and branching collaterals. Axons are unmyelinated and their ends are sealed. Processes may extend outside the displayed soma volume; they remain part of the calculation.

This is the discrete extracellular cable formulation underlying the activating-function approach (Rattay). A spatially constant φ cancels exactly. A uniform electric field can still polarize the ends of a sealed cable. Consequently, neither |E| nor extracellular voltage alone is used as an activation threshold.

Gate kinetics are independently implemented from the equations used by Pospischil et al. and checked against the authors’ ModelDB implementation. For u = V − V_T, rates in ms⁻¹ are:

Temperature is fixed at 36 °C. Sodium reversal is +50 mV; potassium reversal is −100 mV; specific capacitance is 1 µF/cm². The original single-compartment models use an equivalent membrane area, not a realistically sized anatomical soma. This lab uses physical soma dimensions and a newly specified compartment distribution, so it is not a reproduction of those fitted cell models.

The spatial allocation, axonal threshold shift and morphology are model assumptions. Leak conductance is 0.1 mS/cm². Each compartment’s leak reversal is set so the steady-state gates balance at −70 mV, eliminating artificial initialization spikes. Intracellular resistivity defaults to 150 Ω·cm. Adjacent half-compartment axial resistances add in series.

Gate updates use exponential Rush–Larsen steps with linear interpolation in a 0.1 mV voltage table. Voltage and axial coupling use backward Euler with a tree-structured Hines elimination. The default step is 5 µs and maximum neurite segment is 20 µm; 50 µs pulses automatically use at most 1 µs. These resolutions are informed by Aberra et al., whose much more detailed reconstructed models are not replicated here.

A spike event is the first upward crossing of 0 mV in a soma, AIS or axonal compartment. Dendritic crossings do not independently mark recruitment. Soma and axon times are retained separately. The timing raster shows first events; pulse-train responses remain visible in the voltage traces. Excursions beyond |V_m| = 200 mV are flagged as outside the reliable range of this reduced membrane model.

4. Anatomical scale and sampling

Garcia-Marin et al. report 99,587 neurons/mm³ in human V1 L2/3, with 722 µm layer thickness beneath a 281 µm L1. The default cube spans 350–950 µm depth and contains round(99,587 × 0.6³) = 21,511 somata. Their estimated 13% inhibitory fraction is not a direct human GABA count. The experiment maps that estimate to a single FS-like class; real inhibitory populations are more diverse. The paper’s 2024 correction concerns its processing-module summary table.

Benavides-Piccione et al. report a 154 µm² mean maximum soma profile for human BA17 L3a pyramidal cells. The model’s 14 µm profile-equivalent diameter is derived as 2√(154/π), not a measured spherical diameter. Somata are ellipsoids; the 18% lognormal size dispersion, shape ratios and 0.85 inhibitory size multiplier are explicit assumptions. The same study places proximal somatic axons near 1.25 µm, tapering to roughly 0.6 µm; apical and basal processes are also kept on micrometer scales.

The 40 µm AIS is an idealized choice informed by human cortical measurements, including approximately 39.5 µm mean length in Ostos et al. Their human data are temporal cortex, not V1. Apical direction points toward pia with modest jitter. Randomized orientation is a sensitivity experiment, not an alternative anatomical claim.

A seeded rejection process places variable-sized somata without overlap, using a conservative bounding-sphere clearance. This is not a measured cortical point pattern, column reconstruction or connectome. Cell sizes and classes are drawn before placement, and larger somata are placed first to avoid redrawing difficult-to-place cells as smaller ones. The run reports accepted mean diameter, class fraction and achieved density. If the placement limit is reached, it reports the smaller achieved population instead of claiming the requested count. The characteristic spacing ρ^−1/3 is about 21.6 µm; it is not nearest-neighbor distance. The run reports a separate nearest-neighbor estimate from a 256-cell placement audit. The full density population is shown with physical soma size; dimming does not enlarge cells.

A deterministic random subset receives full membrane calculations. Increasing sample count preserves the earlier subset. Cells whose soma or sampled neurite paths intersect a contact or its 2 µm clearance are excluded and counted. The remaining unsolved background has no implied spike state. Recruitment bars are sample fractions, not exact population counts; Wilson intervals describe finite sampling only. Soma-distance extent is measured within the sampled soma volume and is not an all-brain activation radius.

5. Numerical checks and reproducibility

The committed scientific tests check the solver against analytical cases and invariants. Passing them establishes numerical behavior under those cases. It does not validate this reduced model’s absolute thresholds against human experiments.

Those tolerances concern the committed benchmark, not every contact arrangement. Near threshold, a small numerical or biological change can alter whether a neuron fires. Use the refinement controls on any conclusion that matters. Geometry, source integration and membrane equations are separated in the source, with computations performed in a cancellable browser worker.

Run JSON exports the completed configuration, model version, anatomical summary and interface accounting. Neuron CSV includes positions, type, soma size, nearest soma/neurite distances and first response times. Trace CSV exports the selected cell’s soma/AIS voltages and every contact’s current and interface voltages. Traces are sampled approximately every 25 µs; spike events are detected at the integration time step. Loading JSON restores configuration without silently starting a new run.

6. What can and cannot be concluded

This is a mechanistic, reduced simulation with literature-based anatomical scale. It can show how signed sources superpose, how contact size changes charge density, how compliance changes a waveform, and why neurite geometry can make recruitment spatially unintuitive. Its absolute neuron counts and thresholds are not calibrated predictions for a human implant.

• The tissue has no CSF, white-matter boundary, cortical folding, vasculature or electrode encapsulation. Conductivity is spatially constant and frequency independent.
• Uniform surface sources omit edge crowding and the insulating structures of commercial probes. The sheet geometry is not a flush contact on a one-sided insulated substrate.
• Local idealized arbors omit many branches, dendritic spines, long-range axons and myelin. Artificial sealed terminals can affect threshold and initiation site. Full reconstructed morphologies are needed for stronger anatomical claims.
• The anatomical references are human; the channel kinetics derive from simplified cortical model classes and are not fitted to these human V1 cells. Spatial conductance distributions require independent calibration.
• There are no synapses, recurrent excitation/inhibition, spontaneous activity, ephaptic feedback or plasticity. Density changes the represented population and spacing; it does not change excitability through network coupling.
• Electrochemical circuit constants are linear equivalents. Nonlinear Faradaic reactions, diffusion, water-window limits, heating and injury are not calculated.
• The visual field is a section through the volume, while neuron responses depend on their full modeled arbor. Interpreting a heat-map contour as an activation boundary would be incorrect.

Near-contact predictions deserve particular caution: the tool flags selected cells with neurites within 30 µm. The explicit finite contacts remove a point-source singularity, but do not resolve membrane-scale extracellular geometry. Detailed stimulation studies also treat electrode–neurite proximity as a modeling limitation.

7. Extending this experiment to optogenetics

The anatomy generator and membrane solver are separate from the extracellular drive. A later optical experiment can reuse the same cells while replacing that drive with a light-transport calculation, wavelength-dependent scattering and absorption, cell-specific opsin expression, and photocycle conductances. This release implements electrical stimulation only. The useful comparison will be recruited cell identity and timing under matched anatomical assumptions, with optical heating and expression uncertainty included explicitly.

Primary sources

References anchor specific equations or measurements; citing a detailed model does not mean that this reduced implementation reproduces it.

1. Garcia-Marin, Kelly & Hawken (2024). Neuronal composition of processing modules in human V1. 10.1093/cercor/bhad512.
2. Benavides-Piccione et al. (2024). Key morphological features of human pyramidal neurons. 10.1093/cercor/bhae180.
3. Ostos et al. (2023). GABAergic innervation of the soma and axon initial segment in human and mouse neocortex. 10.1093/cercor/bhac314.
4. Rattay (1986). Analysis of models for external stimulation of axons. 10.1109/TBME.1986.325670.
5. Pospischil et al. (2008). Minimal Hodgkin–Huxley type models for different classes of cortical and thalamic neurons. 10.1007/s00422-008-0263-8.
6. Aberra, Peterchev & Grill (2018). Biophysically realistic neuron models for simulation of cortical stimulation. 10.1088/1741-2552/aadbb1.
7. Histed, Bonin & Reid (2009). Direct activation of sparse, distributed populations of cortical neurons by electrical microstimulation. 10.1016/j.neuron.2009.07.016.
8. Joucla & Yvert (2009). Improved focalization of electrical microstimulation using microelectrode arrays: a modeling study. 10.1371/journal.pone.0004828.
9. Hindriks et al. (2017). Linear distributed source modeling of local field potentials recorded with intra-cortical electrode arrays. 10.1371/journal.pone.0187490.
10. Wei & Grill (2009). Impedance characteristics of deep brain stimulation electrodes in vitro and in vivo. 10.1088/1741-2560/6/4/046008.
11. Howell, Naik & Grill (2014). Influences of interpolation error, electrode geometry, and the electrode–tissue interface…. 10.1109/TBME.2013.2292025.
12. Cogan et al. (2009). Sputtered iridium oxide films for neural stimulation electrodes. 10.1002/jbm.b.31223.
13. Leung et al. (2015). In vivo and in vitro comparison of the charge injection capacity of platinum macroelectrodes. 10.1109/TBME.2014.2366514.
14. Weiland, Anderson & Humayun (2002). In vitro electrical properties for iridium oxide versus titanium nitride stimulating electrodes. 10.1109/TBME.2002.805487.