Simulation of transcranial magnetic stimulation in head model with morphologically-realistic cortical neurons

Abstract

Background: Transcranial magnetic stimulation (TMS) enables non-invasive modulation of brain activity with both clinical and research applications, but fundamental questions remain about the neural types and elements TMS activates and how stimulation parameters affect the neural response.

Objective: To develop a multi-scale computational model to quantify the effect of TMS parameters on the direct response of individual neurons.

Methods: We integrated morphologically-realistic neuronal models with TMS-induced electric fields computed in a finite element model of a human head to quantify the cortical response to TMS with several combinations of pulse waveforms and current directions.

Results: TMS activated with lowest intensity intracortical axonal terminations in the superficial gyral crown and lip regions. Layer 5 pyramidal cells had the lowest thresholds, but layer 2/3 pyramidal cells and inhibitory basket cells were also activated at most intensities. Direct activation of layers 1 and 6 was unlikely. Neural activation was largely driven by the field magnitude, rather than the field component normal to the cortical surface. Varying the induced current direction caused a waveform-dependent shift in the activation site and provided a potential mechanism for experimentally observed differences in thresholds and latencies of muscle responses.

Conclusions: This biophysically-based simulation provides a novel method to elucidate mechanisms and inform parameter selection of TMS and other cortical stimulation modalities. It also serves as a foundation for more detailed network models of the response to TMS, which may include endogenous activity, synaptic connectivity, inputs from intrinsic and extrinsic axonal projections, and corticofugal axons in white matter.

Keywords: Finite element method; Motor cortex; Neuron models; Simulation; Transcranial magnetic stimulation.

Figure 1.. Threshold dependence on cell type and direction of uniform E-field.

A) Coordinate system shown for example cell (L2/3 PC). Somato-dendritic axis aligned to z-axis, with polar angle θ and azimuthal angle ϕ. B) Uniform E-field directions represented as normal vectors on sphere centered at origin. Thresholds were calculated for 398 directions spanning the sphere, and each threshold value was represented as a point on the sphere (white dot) corresponding to the E-field vector $\stackrel{\rightharpoonup }{E}$ C) 3D threshold–direction map projected into 2D using Mollweide projection. White dot indicates threshold value for example vector $\stackrel{\rightharpoonup }{E}$ in B, crossed circle represents E-field pointing into the page, and circle with dot represents E-field pointing out of the page. Bottom: Recorded MagPro X100 Monophasic TMS waveform. D) Threshold–direction maps for all cell types and their virtual clones, i.e. models of the same cell type with stochastically varied morphologies, normalized to the minimum for each clone. Within cell type, threshold–direction maps for each clone are ordered by minimum threshold, with lowest minimum threshold at the bottom. White star denotes minimum threshold orientation. Corresponding cell morphology plotted to the right of each map with same color scheme as in A. Black arrow points in direction of minimum threshold orientation, matching white star. All scale bars are 250 μm. E) Minimum threshold for 5 clones of each cell type, grouped by layer.

Figure 2.. Embedding populations of cortical neuron models in FEM models of TMS induced E-fields.

A) Scalp and gray matter meshes are shown with the overlying TMS coil outline. The coil center and orientation are given by the green sphere and arrow, the hand knob region populated with neurons is indicated in red, and the putative hand muscle representation used in Figure 7 is shown in blue. B) Model neurons located in the crown of the pre-central gyrus between the gray matter and white matter surfaces. One clone from each layer is shown with color corresponding to layer (shown in C); five co-located model neuron populations (virtual clones of each cell type) are simulated in each layer. The 2D analysis plane (red) is used to visualize threshold data in Figure 5 and also shown in C–D. The plane is parallel to coil orientation (45° relative to midline). C) Cortical layers used to place and orient model neurons shown in 2D analysis plane, extracted by intersecting the analysis plane with the layer surfaces. D) Neural populations from B visualized within the 2D analysis plane. Neurons in all five layers are plotted with their respective layer colors (left). L2/3 (middle) and L5 (right) PC populations are plotted with axon, apical dendrites, and basal dendrites colored separately.

Figure 3.. Layer-specific spatial distribution of activation thresholds correlate better with the E-field magnitude than normal component.

A) Magnitude of simulated E-field (normalized across layers) on layer surfaces for L1–L6 arranged adjacent to each other. B) Component of E-field normal to layer surfaces (normalized across layers). Positive values indicate E-field pointing out of surface and negative values indicate inward E-field. C) Median thresholds (across 5 clones and 6 rotations) for monophasic P–A simulation. D) Inverse threshold of each cell plotted against E-field magnitude (top) and absolute value of normal component (bottom) at soma, both normalized to maximum within layer (in A). Each plot includes R² value for linear regression of inverse threshold with corresponding E-field metric.

Figure 4.. TMS activates axonal terminations aligned to local E-field direction.

Axonal arbors for single population of L2/3 PCs, L4 LBCs, and L5 PCs with directly activated branch colored from AP initiation point (terminal) to proximal branch point for monophasic P–A (green) and A–P (magenta) stimulation. Dendrites not shown. Somas indicated by black dots.

Figure 5.. Layer-specific spatial distribution of activation varies with TMS pulse waveform and direction, shown in a cross-section of the hand knob (as in Figure 2C).

Median thresholds for L1–L6 on analysis plane through pre-central gyrus, parallel to coil handle and near coil center for A) monophasic, B) half sine, and C) biphasic stimulation in the P–A and A–P directions. Arrows indicate direction of initial phase of E-field waveform. Note that biphasic stimulation conditions are plotted in opposite order to group stimuli by the direction of their dominant waveform phase. D) Percent difference in median thresholds between P–A and A–P current directions, indicated in title. Regions where thresholds for both P–A and A–P were above 230 A/μs are colored gray.

Figure 6.. Simpler neuron model with straight axon morphology results in higher thresholds and different spatial distribution of activation for L2/3 PCs.

A) Median thresholds for L2/3 PCs on layer surface for monophasic P–A stimulation (top; compare to Figure 3C) and A–P stimulation (middle). Placement of L2/3 PCs with straight axons visualized in gyral crown between the gray matter and white matter surfaces (similar to Figure 2B) and zoomed in view of single model neuron (bottom). Straight axon versions of L2/3 PCs were co-located with original, realistic axon population and only differed in their axon morphology. B) Median thresholds for L2/3 PC on analysis plane (same as Figure 5) for monophasic P–A (top) and A–P (bottom) stimulation. C) Percent difference in median thresholds between P–A and A–P current directions. Regions where thresholds for both P–A and A–P were above 230 A/μs are colored gray.

Figure 7.. Activation thresholds within putative hand muscle representation in hand knob.

Model thresholds for each pulse waveform and direction combination are shown in log scale, with each Tukey boxplot (outliers excluded) describing statistics of thresholds from 5 clones and 6 rotations at each position within the hand muscle representation. Boxplots of experimental motor threshold (MT) data are included in gray (12 subjects) [41]. Hand muscle representation in L5 surface is marked in blue (right), with green arrow indicating center of TMS coil and direction of induced current (for P–A stimulation).

Figure 8.. Strength–duration time constants of model neurons match experimental measures.

A) Median time constant (across 5 clones and 6 rotations) for layers with lowest activation thresholds, estimated for each model neuron using their activation thresholds for 30, 60, and 120 μs cTMS pulses. Putative hand muscle representations are outlined in black. B) Strength–duration curves using median neuronal population threshold within the model hand muscle representation and experimentally measured mean (± SD) motor thresholds [14]. cTMS waveforms are shown in the bottom left corner. C) Experimentally estimated time constant [14] (mean ± SD) and model time constants for a range of cutoff threshold percentiles (2.5–50%) within the hand muscle representation.