We derived computationally efficient average response models of different types of cortical neurons, which are subject to external electric fields from Transcranial Magnetic Stimulation. We used 24 reconstructions of pyramidal cells (PC) from layer 2/3, 245 small, nested, and large basket cells from layer 4, and 30 PC from layer 5 with different morphologies for deriving average models. With these models, it is possible to efficiently estimate the stimulation thresholds depending on the underlying electric field distribution in the brain, without having to implement and compute complex neuron compartment models. The stimulation thresholds were determined by exposing the neurons to TMS-induced electric fields with different angles, intensities, pulse waveforms, and field decays along the somato-dendritic axis. The derived average response models were verified by reference simulations using a high-resolution realistic head model containing several million neurons. The relative errors of the estimated thresholds between the average model and the reference model ranged between -3% and 3.7% in 98% of the cases, while the computation time was only a fraction of a second compared to several weeks. Finally, we compared the model behavior to TMS experiments and observed high correspondence to the orientation sensitivity of motor evoked potentials. The derived models were compared to the classical cortical column cosine model and to simplified ball-and-stick neurons. It was shown that both models oversimplify the complex interplay between the electric field and the neurons and do not adequately represent the directional sensitivity of the different cell types. The derived models are simple to apply and only require the TMS-induced electric field in the brain as input variable. The models and code are available to the general public in open-source repositories for integration into TMS studies to estimate the expected stimulation thresholds for an improved dosing and treatment planning in the future.

The extension of current models in the area of transcranial brain stimulation beyond the estimation of the electric fields is elementary to improve our understanding of the underlying stimulation processes. The key question is how the electric field modulates the behavior of neuronal structures. Earlier experimental studies showed that the depolarization threshold of isolated straight axons is inversely proportional to the cosine of the angle between the external current and the nerve fiber (Rushton, 1927). This led to the well-known cortical column cosine hypothesis (Fox et al., 2004), assuming that excitable neuronal elements, in particular axons, have a preferential orientation perpendicular to the cortical surface. At first glance, this model seems to be supported by the findings of Rudin and Eisenman (1954) and Ranck (1975), who consistently found that orthodromic currents are more effective than antidromic currents and especially transverse currents. Note, however, that the complex morphology of the neurons does not allow the generalization of observations made in single isolated axons to neuronal populations, because the orientations of the axon segments relative to the external electric field vary and have to be considered statistically. Accounting for these effects requires a model description across multiple scales. This involves first determining the electric field in the brain by solving Maxwell’s equations and then coupling it with detailed mesoscopic neuron models. Aberra et al. (2020) introduced a novel approach to simulate the effects of TMS in head models with morphologically realistic cortical neurons. These authors developed a multi-scale computational model that is capable of quantifying effects of different TMS parameters on the direct response of individual cortical neurons. They created digital representations of neurons that match the geometry and biophysical properties of mature human neocortical cells based on neuronal models of rodent cells from the Blue Brain Project (Markram et al., 2015). These models included a spatial representation of the neuronal compartments as well as experimentally validated electrophysiological parameters (Aberra et al., 2018). They were placed inside the gray matter of a realistic head model and the stimulation thresholds for the generation of action potentials were determined by coupling them with the TMS-induced electric fields. The results provide important mechanistic insights into TMS. However, a major limitation of this modeling approach is its high computational cost, which prevents most routine applications of the method in TMS studies. Moreover, a further challenge is that for estimating the overall threshold of a cortical group of neurons, the results of a large number of single simulated responses need to be determined and averaged. This calls for the development of simpler models of neural populations that still accurately account for the modulation of neuronal states through TMS-induced electric fields.

We developed a parsimonious model, which reproduces the effect of the electric field on cortical neurons with high accuracy for different pulse waveforms and geometric electric field parameters. We adapted and extended the approach of Aberra et al. (2020) to derive an average threshold model of layer 2/3 pyramidal cells (L2/3 PC), small, nested, and large basket cells in layer 4 (L4 SBC, L4 NBC, L4 LBC), and layer 5 pyramidal cells (L5 PC). We adapted the pipeline of Aberra et al. (2020) in Python and implemented additional improvements and extensions, such as support for SimNIBS 4 and the CHARM head modeling pipeline (Puonti et al., 2020). The code and associated example scripts are published in the open-source Python package TMS-Neuro-Sim (https://github.com/TorgeW/TMS-Neuro-Sim). Additionally, we determined estimators for the neuronal recruitment rate, which quantifies the relative number of neurons stimulated by TMS at a given stimulation intensity and field orientation.

Concurrently with this study, Aberra et al. (2023) developed a machine-learning model for efficient threshold determination. Both studies are comparable because the same method was used to calculate the stimulation thresholds (Aberra et al., 2020). However, the model presented in Aberra et al. (2023) included only five representations of L2/3 PC, L4 LBC, and L5 PC, each, which we found insufficient to derive accurate average threshold statistics at the population level. More importantly, we went beyond the mere approximation of the threshold estimation, by systematically characterizing the impact of the major parameters, which are the field angle and its gradient across grey matter (GM). Preserving the parameterization allowed us to perform a sensitivity analysis and to identify the most influential parameters of the models by determining so-called Sobol indices using a generalized polynomial chaos expansion (Weise, Poßner, et al., 2020). Moreover, the model was verified by comparing it to results of computationally expensive reference simulations, using a high-resolution realistic head model with a large number of realistically shaped neurons located within the motor cortex. Finally, we validated the model by comparing it with TMS experiments by Souza et al. (2022), who intensively investigated the directional sensitivity of motor evoked potentials using a novel multi-coil TMS transducer.

We also compared the results with those of the cortical column cosine (Fox et al., 2004) as well as a simplified ball-and-stick model (Bédard & Destexhe, 2008) adapted for TMS. It turned out that the stimulation properties differ significantly from detailed neurons and that a simplified modeling strategy is not appropriate in this context.

All data and code underlying the results presented in this paper, together with additional details including the average threshold models, the recruitment rate operators, and the neuron compartment models, are publically available in a repository (Weise, Worbs, et al., 2023), where we provide look-up tables, interpolators, and polynomial approximations for further use.

2.1 Neuron models

To derive the average neuron response models, we extended the set of neural compartment models by Aberra et al. (2020) from originally five neurons of each cell type to 24 L2/3 PC, 70 L4 SBC, 70 L4 NBC, 105 L4 LBC, and 30 L5 thick-tufted pyramidal cells (TTPC’s), taken from the Blue Brain Project (Ramaswamy et al., 2015). The cells originate from the somatosensory cortices of P14 male Wistar (Han) rats (Markram et al., 2015). They were stained with biocytin, visually recorded with a bright-field light microscope, and processed by the software Neurolucida (Williston, VT, USA). Shrinkage due to staining in the z-axis was corrected during the reconstruction. In an unraveling step, shrinkage in the xy-axis was corrected for with a method based on the centered moving window algorithm by smoothing and extending the reach of the branches while maintaining their overall length (Anwar et al., 2009). For branch repair, the cutting planes were first determined and the cut branches were then statistically regrown based on the intact branches. Because some resulting cell morphologies contain impoverished axonal/dendrite branching, a mix-and-match procedure was used to create cells with valid dendrite and axonal reconstructions. As a last step to increase morphological diversity, a cloning procedure was applied. The procedure assigns distributions to branch length and rotation while preserving the overall branching structure.

Because the cells provided were from rats, further modifications were necessary to obtain human-like neurons. We followed the procedure and parameters given in Aberra et al. (2018) to extend the set of neurons. First, the basal dendritic diameter, basal dendritic length, apical dendritic diameter, somatic diameter, and axonal diameter were scaled to create adult human-like neuron morphologies. Second, the axons were myelinated by registering nodes of Ranvier with a width of 1 μm, creating myelinated sections with a length (L) to diameter (D) ratio of L/D = 100 and myelinated axon terminals with L/D=70 (Hess & Young, 1949; Hursh, 1939; Waxman & Kocsis, 1995). And third, the ion channel properties were adapted according to the myelination (see Table 1 in Aberra et al., 2018). Figure 1 provides an overview of the cells used in the study. The average numbers of nodes per cell are 3,541 for L2/3 PC, 14,779 for L4 SBC, 13,091 for L4 NBC, 9,147 for L4 LBC, and 12,514 for L5 PC. For the compartment models, the neurons were discretized with a maximum compartment length of 20 μm. This resulted in an average number of compartments of 766 for L2/3 PC, 1,447 for L4 SBC, 1,762 for L4 NBC, 1,876 for L4 LBC, and 2,008 for L5 PC.

Fig. 1.

Example morphologies of L2/3 PC, L4 S/N/LBC, and L5 PC: The numbers below the cells indicate the corresponding IDs in the repository Weise, Worbs, et al. (2023). L4 BC are categorized in small basket cells (SBC), nested basket cells (NBC), and large basket cells (LBC). In total, the study includes 24 L2/3 PC, 70 L4 SBC, 70 L4 NBC, 105 L4 LBC, and 30 L5 thick-tufted pyramidal cells (TTPC’s), taken from the Blue Brain Project (Ramaswamy et al., 2015).

Fig. 1.

Example morphologies of L2/3 PC, L4 S/N/LBC, and L5 PC: The numbers below the cells indicate the corresponding IDs in the repository Weise, Worbs, et al. (2023). L4 BC are categorized in small basket cells (SBC), nested basket cells (NBC), and large basket cells (LBC). In total, the study includes 24 L2/3 PC, 70 L4 SBC, 70 L4 NBC, 105 L4 LBC, and 30 L5 thick-tufted pyramidal cells (TTPC’s), taken from the Blue Brain Project (Ramaswamy et al., 2015).

Close modal
Table 1.

Sobol indices of the electric field threshold models for L2/3, L4 S/N/LBC, and L5 PC for monophasic and biphasic pulse waveforms.

Cell typeL2/3 PCL4 SBCL4 NBCL4 LBCL5 PC
Wave-formMono-phasicBi-phasicMono-phasicBi-phasicMono-phasicBi-phasicMono-phasicBi-phasicMono-phasicBi-phasic
θ 0.919 0.710 0.991 0.994 0.974 0.985 0.962 0.974 0.951 0.971 
Δ|E˜| 0.052 0.247 0.001 0.001 0.001 0.001 0.003 0.002 0.037 0.017 
θ & Δ|E˜| 0.029 0.043 0.008 0.005 0.025 0.014 0.035 0.024 0.012 0.012 
Cell typeL2/3 PCL4 SBCL4 NBCL4 LBCL5 PC
Wave-formMono-phasicBi-phasicMono-phasicBi-phasicMono-phasicBi-phasicMono-phasicBi-phasicMono-phasicBi-phasic
θ 0.919 0.710 0.991 0.994 0.974 0.985 0.962 0.974 0.951 0.971 
Δ|E˜| 0.052 0.247 0.001 0.001 0.001 0.001 0.003 0.002 0.037 0.017 
θ & Δ|E˜| 0.029 0.043 0.008 0.005 0.025 0.014 0.035 0.024 0.012 0.012 

2.2 Coupling of electric fields into neuron models

The electric field E(z,t) caused by TMS generates an additional extracellular pseudo-potential φe(z,t) (Wang et al., 2018). It is coupled into the neurons’ cable equations by integrating the electric field component along the local longitudinal direction dz of each neuronal compartment, ranging from, for example, the initial point of a compartment z0 to the end of that compartment z:

φe(z,t)=z0zE(z,t)dz+φe(z0,t)
(1)

For the realistic head model simulations, the electric field is interpolated to the neurons’ segments using the superconvergent patch recovery approach (Zienkiewicz & Zhu, 1992).

2.3 Neuronal simulations

The stimulation behavior of the neurons is analyzed by calculating the transmembrane potential in each compartment using NEURON (Carnevale & Hines, 2006) following a similar approach as Aberra et al. (2020). The spatio-temporal dynamics of the transmembrane potential were modeled according to the Hodgkin-Huxley formalism. Detailed information about the ion channel parameters and membrane time constants can be found in the repository by Weise, Worbs, et al. (2023) and ModelDB (https://senselab.med.yale.edu/modeldb/ShowModel.cshtml?model=241165) by Aberra et al. (2018). The NEURON simulation was set up with a temperature of 37°C and an initial voltage of -70 mV for each compartment. The simulations were carried out over the course of 1 ms with time steps of 5 µs. The extracellular quasipotentials were scaled by the waveform and the amplitude of the TMS pulse. The used monophasic and biphasic waveforms were taken from a MagPro X100 stimulator (MagVenture A/S, Denmark) with a MagVenture MCF-B70 figure-of-eight coil (P/N 9016E0564) and were recorded using a search coil with a sampling rate of 5 MHz. The recordings were down-sampled to the simulation time steps and normalized to be applicable for scaling the extracellular potential. The cell thresholds are determined as the minimum electric field intensity needed to elicit action potentials in at least three compartments, using a binary search approach with a precision of 0.05 V/m. The simulation environment is implemented and published in the open-source Python package TMS-Neuro-Sim (https://github.com/TorgeW/TMS-Neuro-Sim), making use of the Python API of NEURON. The example dataset (Weise, Worbs, et al., 2023) contains the neuron models together with example scripts detailing the use of the implemented functions.

2.4 Average response model of cortical neurons

We exposed the model neurons to electric fields from different directions and strengths to examine their stimulation behavior in detail. We parameterized the electric field direction using spherical coordinates (Fig. 2a). The polar angle θ quantifies the angle between the electric field and the somato-dendritic axis (z-axis) and ranges from 0° to 180°. The azimuthal angle φ quantifies the electric field direction in the horizontal plane perpendicular to the somato-dendritic axis and ranges between 0° and 360°. The coordinate system is defined such that the soma is lying close to the center and the axon extends into the negative z-direction. Because of the comparatively large extension of the PC from the uppermost dendrites to the lowermost part of the axon, the decay of the electric field along the z-axis is not negligible. In simulations of a realistic head model, we found that the electric field can differ up to ±20% per mm over the somato-dendritic axis. A more detailed analysis of the underlying parameter distributions is given later in Section “Sensitivity analysis.” For this reason, we have added an additional parameter to the model, namely the relative change of the electric field magnitude per unit length Δ|E˜| measured in %/mm.

Fig. 2.

Neurons in an external electric field: (a) Parametrization of the TMS induced electric field relative to cortical neurons; (b) Example of an L2/3 PC, which is exposed to an external electric field with direction ϕ=0°, θ=135° and a field decay of Δ|E˜|=30%/mm. Note that the electric field is stronger in the upper part of the cell and decreases with depth, as is generally observed in the cortex.

Fig. 2.

Neurons in an external electric field: (a) Parametrization of the TMS induced electric field relative to cortical neurons; (b) Example of an L2/3 PC, which is exposed to an external electric field with direction ϕ=0°, θ=135° and a field decay of Δ|E˜|=30%/mm. Note that the electric field is stronger in the upper part of the cell and decreases with depth, as is generally observed in the cortex.

Close modal

The electric field at each location (x,y,z) at firing threshold is then given by:

E(z,ϕ,θ,Δ|E˜|)=Ethres(sinθcosϕsinθsinϕcosθ)max(Δ|E˜|100(zzsoma)+1,0)
(2)

where zsoma is the position of the soma on the z-axis (in mm). The soma will have an electric field magnitude of Ethres, which is to be found using the aforementioned binary search approach, while the magnitude for every other section of the cell is linearly interpolated based on their z-coordinate. An example of an external electric field distribution with ϕ=0°, θ=135° and Δ|E˜|=30%/mm is shown in Fig. 2b.

For the derivation of an average response model, all neurons were exposed to an external electric field with a polar angle θ (range: [0, 180]°, steps: 3°), an azimuthal angle ϕ (range: [0, 360]°, steps: 6°), and a relative change of the electric field along the somato-dendritic axis Δ|E˜| (range: [-100, 100] %/mm, steps: 10 %/mm) for both monophasic and biphasic pulse waveforms. After determining the electric field thresholds for each cell for all possible electrical field configurations, an average threshold model was derived by averaging the thresholds over all compartment models and over all azimuthal orientations ϕ, based on the assumption that the spatial locations and tangential orientations of the neurons in the cortex are random.

From the activation thresholds of the individual neurons, we determined the recruitment rate of the neurons in dependence of θ and Δ|E˜|. The recruitment rate estimates the relative number of neurons which were stimulated by TMS at a given stimulation intensity, with zero corresponding to no stimulation and one corresponding to stimulation of all neurons. To this end, we integrated the electric field thresholds along the electric field magnitude and smoothed the discrete behavior by fitting continuous sigmoid functions of the following type in the least squares sense:

f(E)=(1+er(θ,Δ|E˜|)(EE0(θ,Δ|E˜|)))1
(3)

where E denotes the electric field, r(θ,Δ|E˜|) the slope, and E0(θ,Δ|E˜|) the shift of the sigmoidal functions, which depend on, both, the polar angle θ and the relative change in electric field Δ|E˜|.

A threshold model was also created for a simplified ball-and-stick neuron model. Typically, ball-and-stick models are used for stimulation by weak electric fields in the context of transcranial electric stimulation (e.g., Aspart et al., 2016) and consist of one segment for the dendrites and one segment for the soma. Because the stimulation thresholds for TMS-induced electric fields of the dendrites are more than 10 times higher than those of the axons, the classical ball-and-stick model had to be modified for TMS. For this purpose, we integrated a straight axon instead of dendrites into the model and determined the stimulation thresholds as a function of the polar angle θ using the approach described by 0. For a similar approach, see also the Supplementary Material of that article. The parameters (axon length and diameter) were adapted such that the thresholds matched those of our complex model.

2.5 Sensitivity analysis

A sensitivity analysis of the derived threshold maps was conducted in terms of variations of the electric field parameters θ and Δ|E˜|. We derived a generalized polynomial chaos expansion of the threshold maps using the Python package pygpc (Weise, Poßner, et al., 2020) and determined the first- and second-order Sobol indices that quantify the fraction of the total variance of the threshold that stems from the variance of θ, from the variance of Δ|E˜|, and from the combination of both. The input distributions of both parameters were estimated from the electric field simulation of the high-resolution realistic head model. We extracted the polar angles θ and the changes in electric field magnitude Δ|E˜| in every surface element of layer 5 in the ROI and fitted uniform and beta distributions to the histograms (Fig. 3). We repeated the analysis for layer 2/3 and layer 4, and did not find any major differences in the parameter distributions, due to the close proximity of the layers. For the uncertainty analysis, we assumed that both parameters are uncorrelated.

Fig. 3.

Distribution of electric field parameters on layer 5 in a realistic head model: Histograms and fitted uniform and beta distributions of (a) the polar angle θ (uniform parameters: θmin=0°, θmax=180°; beta parameters: θmin=0°, θmax=180°; p=1.51, q=1.56) and (b) the relative change of the electric field magnitude Δ|E˜| (uniform parameters: Δ|E˜|min=15%, Δ|E˜|max=15%; beta parameters: Δ|E˜|min=30%, Δ|E˜|max=30%; p=13.86, q=13.78).

Fig. 3.

Distribution of electric field parameters on layer 5 in a realistic head model: Histograms and fitted uniform and beta distributions of (a) the polar angle θ (uniform parameters: θmin=0°, θmax=180°; beta parameters: θmin=0°, θmax=180°; p=1.51, q=1.56) and (b) the relative change of the electric field magnitude Δ|E˜| (uniform parameters: Δ|E˜|min=15%, Δ|E˜|max=15%; beta parameters: Δ|E˜|min=30%, Δ|E˜|max=30%; p=13.86, q=13.78).

Close modal

It is noted that the gPC surrogate model can be used as a very efficient threshold estimator, which can be written in matrix form as:

[Ψ(Θ,Δ|E˜|)]c=y(Θ,Δ|E˜|)
(6)

where [Ψ(Θ,Δ|E˜|)] is the gPC matrix of size [NROI×Nc]. Each row is calculated by inserting the electric field parameters defined in the vectors Θ and Δ|E˜| of length NROI in the gPC polynomials; c is the gPC-coefficient vector of length Nc, and y(Θ,Δ|E˜|) is the vector which contains the estimated thresholds in each point of the ROI. The current 15th-order gPCs with two parameters consist of Nc=136 coefficients.

2.6 Model verification

In order to verify the average response model, we conducted reference simulations using a high-resolution realistic head model, where we explicitly placed the neurons in the ROI and coupled the TMS-induced electric field into them. The head model was created using T1-, T2-, and diffusion-weighted MRI. The images were acquired on a 3 T MRI scanner (Siemens Skyra) with a 32 channel head coil using the same acquisition parameters as described in Weise, Numssen, et al. (2023). T1 and T2 images were used for tissue-type segmentation. Conductivity tensors in gray and white matter were reconstructed from diffusion-weighted images using the volume normalized mapping approach (dwi2cond, https://simnibs.github.io/simnibs/build/html/documentation/command_line/dwi2cond.html, Güllmar et al., 2010). The head model was generated using the headreco pipeline (Nielsen et al., 2018) utilizing SPM12 (https://www.fil.ion.ucl.ac.uk/spm/software/spm12/, Penny et al., 2011) and CAT12 (http://www.neuro.uni-jena.de/cat/, Gaser et al., 2022). A region of interest (ROI) was defined around the handknob area (FreeSurfer, http://surfer.nmr.mgh.harvard.edu/, Dale et al., 1999; Fischl et al., 1998) based on the fsaverage template. This covered parts of somatosensory cortex (BA1, BA3), primary motor cortex M1 (BA4), and dorsal premotor cortex PMd (BA6). The head model was refined in the ROI to provide accurate electric field values for the neuron models (Fig. 4). The final head model is composed of ~5∙106 nodes and ~29∙106 tetrahedra. The tetrahedra in the ROI have an average edge length of 0.45 mm and an average volume of 0.01 mm³. The model consists of six tissue types with the following electrical conductivities: white matter (0.126 S/m), gray matter (0.275 S/m), cerebrospinal fluid (1.654 S/m), bone (0.01 S/m), skin (0.465 S/m), and eyeballs (0.5 S/m) (Thielscher et al., 2015; Wagner et al., 2004). The entire process from MRI acquisition to the final head model can be reproduced in detail using the protocol of Weise, Numssen, et al. (2023) (steps 1-20) and details of the FEM are given in Saturnino et al. (2019).

Fig. 4.

Realistic head model with cortical layers and neurons: The model was constructed with SimNIBS v3.2.6 (Thielscher et al., 2015) using headreco (Nielsen et al., 2018). In the M1 region, the cortical layers 1-6 are generated and the mesh is refined to ensure highly accurate electric field profiles, which are coupled into the compartment models of cortical neurons. The bottom right inset shows an example of the magnitude of the electric field as color code and its orientation as white streamlines. The black arrow indicates the coil orientation.

Fig. 4.

Realistic head model with cortical layers and neurons: The model was constructed with SimNIBS v3.2.6 (Thielscher et al., 2015) using headreco (Nielsen et al., 2018). In the M1 region, the cortical layers 1-6 are generated and the mesh is refined to ensure highly accurate electric field profiles, which are coupled into the compartment models of cortical neurons. The bottom right inset shows an example of the magnitude of the electric field as color code and its orientation as white streamlines. The black arrow indicates the coil orientation.

Close modal

In order to place the neurons at the right locations in the cortex, we added cortical layers to the head model. The normalized depths of the six cortical layers range between 0 (gray matter surface, i.e., pia mater) and 1 (white matter surface) and were estimated from primate motor cortex slices (García‐Cabezas & Barbas, 2014). The normalized depths of the layer centers are 0.06 for layer 1, 0.4 for layer 2/3, 0.55 for layer 4, 0.65 for layer 5, and 0.85 for layer 6 (Aberra et al., 2020). We linearly interpolated between the white matter surface (1) and the gray matter surface (0) using the vertex positions of the two surfaces. To extract the cortical layers as isosurfaces from the 3D interpolation, we used a marching cubes algorithm (Fig. 3) (Lorensen & Cline, 1987). In every ROI surface element (size ~1 mm²), we placed all cells and rotated them from 0° to 360° in steps of 6°. This resulted in a total number of 12,947,040 L2/3 PC, 130,080,300 L4 S/N/LBC, and 15,760,800 L5 PC in the ROI to simulate. The total pure simulation time using 48 cluster nodes with 72 cores each (Intel Xeon Platinum 8360Y, 256 GB RAM) was approximately 40 days for both monophasic and biphasic waveforms.

The electric field calculations were conducted using SimNIBS v3.2.6 (Saturnino et al., 2019; Thielscher et al, 2015) using a regular figure-of-eight coil (MCF-B65, Magventure, Farum, Denmark), which is placed over the M1 region with an orientation of 45° towards the fissura longitudinalis. The angle θ of the electric field was calculated with respect to the surface normal of the cortical layers in the ROI. Likewise, the percentage change of the electric field magnitude between the WM and GM surfaces Δ|E˜| was calculated by extracting the field at a normalized depth of 10% of the distance between the current layer and the WM and GM surface, respectively, in order to avoid numerical inaccuracies close to the tissue boundaries. The simulation time of the electric field was relatively short compared to the NEURON simulations and took a few seconds.

To quantify the differences between the models, we determined the normalized root mean square error (NRMSE):

NRMSE=i=1NROI(yi,refyi)2i=1NROIyi,ref2
(3)

where yi,ref denotes the thresholds of the reference simulations in the i-th ROI element and yi the thresholds from the average model. Additionally, we calculated the mean absolute percentage error (MAPE) quantifying the prediction accuracy of the average models:

MAPE=1NROIi=1NROI|yi.refyiyi,ref|
(4)

and the coefficient of determination (R2) quantifying the proportion of the total variance explained by the average model:

R2=1i=1NROI(yiyi,ref)2i=1NROI(yiy¯)2
(5)

where y¯ is the mean of the average threshold model.

2.7 Comparison to directional sensitivity of motor evoked potentials

We compared the derived recruitment models to experimental observations from Souza et al. (2022). TMS experiments were conducted to investigate the orientation selectivity of neuronal excitability using a novel two-coil multi-channel TMS transducer for manipulating the electric field orientation. The advantage of the used two-coil TMS transducer is the possibility to precisely manipulate the pulse orientation electronically with high accuracy (~1°), without physically moving the transducer. To measure the effect of the electric field orientation on the motor evoked potential (MEP) amplitude, five single TMS pulses were applied to the abductor pollicis brevis (APB) muscle hotspot at each of 120 different pulse orientations (0–360°; steps of 3°) with a stimulation intensity of 110% of the resting motor threshold (rMT). The MEPs from the APB muscle were recorded from 11 subjects (mean age: 30 years, range 24-41; four women) using surface electromyography electrodes with a belly-tendon montage. TMS pulses had a trapezoidal monophasic waveform (timings: 60 µs of rising, 30 µs of hold, 43.2 µs of falling) and were delivered using custom power electronics. The interstimulus intervals were pseudo-randomized following a uniform distribution between 4 and 6 s. In two other subjects (ages 31 and 36 years; two men; right-handed), the experiment was repeated with stimulation intensities of 110%, 120%, 140%, and 160% rMT. The order of the orientations and intensities of the pulses was pseudo-randomized. A detailed description of the experimental procedure and TMS hardware is given in Souza et al. (2022).

We hypothesize that the sigmoidal-shaped recruitment rate correlates with the I/O curve of the MEPs, which can be also described by a sigmoid. We correlate the orientation sensitivity of the MEPs with the estimates from the recruitment rate at different electric field intensities at varying angles θ. Therefore, we determined the equivalent intensities of the electric field at 110%, 120%, 140%, and 160% rMT. Subsequently, we read out the recruitment rate at the four intensities and compared the courses of the recruitment rate with those of the MEPs.

3.1 Recruitment order and relative threshold ranges

For each cell type investigated, different stimulation thresholds were observed depending on the polar angle θ and the relative change of the electric field over the somato-dendritic axis Δ|E˜|. In Fig. 5a, an overview of the thresholds of all investigated cell types is shown together with the 95th percentiles of the confidence intervals of the means assuming a constant electric field along the somato-dendritic axis (Δ|E˜|=0%/mm). Congruent to that, Fig. 5b shows the threshold ranges, relative to the mean of L5 PCs. It is observed that L5 PCs are recruited first. The L4 LBCs have the second lowest thresholds followed by the L2/3 PCs and the L4 NBCs. The small basket cells are directly stimulated only at higher stimulation intensities. An analogous observation was also made for biphasic TMS pulses and the results are reported in Fig. S10 in the Supplementary Material.

Fig. 5.

Recruitment order and relative threshold ranges of pyramidal and basket cells for monophasic TMS excitation: (a) Average thresholds and 95th percentiles of the confidence intervals of the means assuming a constant electric field along the somato-dendritic axis (Δ|E˜|=0%/mm); (b) Threshold ranges of all investigated cell types relative to the mean of L5 PCs. The dots indicate the mean thresholds and the ranges stem from the variability across the polar angle θ from 0° to 180°.

Fig. 5.

Recruitment order and relative threshold ranges of pyramidal and basket cells for monophasic TMS excitation: (a) Average thresholds and 95th percentiles of the confidence intervals of the means assuming a constant electric field along the somato-dendritic axis (Δ|E˜|=0%/mm); (b) Threshold ranges of all investigated cell types relative to the mean of L5 PCs. The dots indicate the mean thresholds and the ranges stem from the variability across the polar angle θ from 0° to 180°.

Close modal

3.2 Stimulation behavior of L5 PCs

The detailed stimulation behavior for L5 PCs in case of monophasic excitation is shown in Fig. 6. For parameter combinations of particular interest, we illustrated the stimulation location on a representative neuron. Lowest thresholds can be observed when the electric field is parallel to the somato-dendritic axis. This effect is enhanced for positive electric field changes, that is, when the electric field increases from the dendrites to the lower parts of the axons.

Fig. 6.

Stimulation behavior of L5 PCs for monophasic excitation: (a) Threshold map in dependence of the polar angle θ and the relative change of the electric field over the somato-dendritic axis Δ|E˜|. The insets show the locations of excitation, and the red circles indicate the activated terminals. Blue arrows indicate the electric field direction and magnitude; (b) Thresholds of individual neurons for Δ|E˜|=0%/mm along the dashed line in (a). The blue area shows the 95th percentile of the confidence interval of the mean. The equivalent cortical column cosine model is y(θ)=y^|cos(θ)|1 with y^=233.66V/m (dashed line); the axon parameters of the equivalent ball-and-stick model are l=760μm and d=15μm (dotted line); (c) Recruitment rate for Δ|E˜|=0%/mm derived from the individual neuron activation in (b) by integrating over the electric field thresholds. The dashed line indicates the electric field intensity where the recruitment rate is 0.5. The recruitment rate for negative and positive field decays, that is, Δ|E˜|=20%/mm and Δ|E˜|=20%/mm is shown in Fig. S15 in the Supplementary Material.

Fig. 6.

Stimulation behavior of L5 PCs for monophasic excitation: (a) Threshold map in dependence of the polar angle θ and the relative change of the electric field over the somato-dendritic axis Δ|E˜|. The insets show the locations of excitation, and the red circles indicate the activated terminals. Blue arrows indicate the electric field direction and magnitude; (b) Thresholds of individual neurons for Δ|E˜|=0%/mm along the dashed line in (a). The blue area shows the 95th percentile of the confidence interval of the mean. The equivalent cortical column cosine model is y(θ)=y^|cos(θ)|1 with y^=233.66V/m (dashed line); the axon parameters of the equivalent ball-and-stick model are l=760μm and d=15μm (dotted line); (c) Recruitment rate for Δ|E˜|=0%/mm derived from the individual neuron activation in (b) by integrating over the electric field thresholds. The dashed line indicates the electric field intensity where the recruitment rate is 0.5. The recruitment rate for negative and positive field decays, that is, Δ|E˜|=20%/mm and Δ|E˜|=20%/mm is shown in Fig. S15 in the Supplementary Material.

Close modal

Compared to the behavior of the other cell types investigated, the L5 PCs have the lowest average thresholds (Fig. 6a). It can be observed that the electric field thresholds are highest when the electric field is approximately normal to the somato-dendritic axis (θ90°). Since in this case the electric field can approach from all azimuthal directions ϕ over which the average was taken, there are several potential stimulation sites. The thresholds are decreasing again when the electric field rotates further until it is pointing antidromic, that is, from the axons to the dendrites (θ=0°). In this case, the activation takes place at cortico-cortical axon branches pointing upwards. This effect is stable in terms of electric field changes along the somato-dendritic axis. It is noted that due to the geometrical relations of the two electric field angles θ and ϕ, the more parallel the electric field is to the somato-dendritic axis (θ0°, θ180°) the less the azimuthal direction ϕ plays a role in the stimulation behavior of the neurons. The thresholds for tangential electric fields (θ=90°) are about 15% higher compared to normal electric fields (θ=0° and θ=180°). The results of the individual neurons in Fig. 6b show that the variance of the thresholds increases with increasing θ. At θ=180°, a cluster of neurons can be identified that have very low stimulation thresholds. These are paralleled by a few neurons that have very high stimulation thresholds compared to this group. This affects the recruitment rate in Fig. 6c, whose 50% level (dashed line) is lower at θ=180° than at θ=0°. The most efficient way to stimulate L5 PCs is the application of electric fields with a polar angle of θ=0° and a negative change in electric field across the somato-dendritic axis (Δ|E˜|<0) or by applying electric fields with an angle of θ=180° together with a positive field change (Δ|E˜|>0). The stimulation behavior is much more diverse for transverse electric fields around θ90° due to the variety of azimuthal angles ϕ in which cortico-cortical and cortico-spinal connections can be stimulated.

The recruitment rate for Δ|E˜|=0%/mm is shown in Fig. 6c. We observed that the shape of the recruitment rate depends on the field decay Δ|E˜|. To illustrate this, the recruitment rate for negative and positive field decays, that is, Δ|E˜|=20%/mm and Δ|E˜|=20%/mm is shown in Fig. S15 in the Supplementary Material. Using the interpolators and data sets provided, the recruitment rate can be calculated for any field decays in the range between Δ|E˜|=[20,20]%/mm.

The results for biphasic excitation are shown in Fig. S9 in the Supplementary Material.

3.3 Stimulation behavior of L2/3 PCs

The results of the average response model of L2/3 PCs in case of a monophasic TMS pulse are shown in Fig. S1 in the Supplementary Material. In Fig. S1a, the electric field thresholds are shown as function of the polar angle θ and the relative change in electric field magnitude Δ|E˜|. The behavior of the 24 individual L2/3 neurons is shown in Fig. S1b for homogeneous electric fields, that is, Δ|E˜|=0%/mm (dashed line in Fig. S1a). The thresholds for tangential electric fields (θ=90°) are about 17% higher compared to normal electric fields (θ=0° and θ=180°). The recruitment rates are shown in Fig. S1c for homogeneous electric fields, that is, Δ|E˜|=0%/mm.

We provide the data of the threshold map from Fig. S1a, the individual neuron behavior from Fig. S1b and beyond, as well as the data of the recruitment rate from Fig. S1c in the associated dataset (Weise, Worbs, et al., 2023). In addition, we provide Python-based SciPy interpolators (Virtanen et al., 2020), whose simple usage is also explained in attached scripts.

The results for biphasic excitation are shown in Fig. S2 in the Supplementary Material.

3.4 Stimulation behavior of L4 BCs

3.4.1 Small basket cells

The results of the average response model of L4 SBCs in case of a monophasic excitation are shown in Fig. S3. A pronounced directional sensitivity can be observed also for this cell type. Again, lowest thresholds can be observed when the electric field is parallel to the somato-dendritic axis (θ=0° and θ=180°). The thresholds are about 17% higher when the external electric field is tangential to the cells (θ=90°). The thresholds are slightly affected if the electric field changes along the somato-dendritic axis (Δ|E˜|0%/mm). Compared to other cells, the average threshold is about 20% and 46% higher for L4 SBCs than for L2/3 PCs and L5 PCs, respectively. The results for biphasic excitation are shown in Fig S4 in the Supplementary Material.

3.4.2 Nested basket cells

The results of the average response model of L4 NBCs in case of a monophasic excitation are shown in Fig. S5. Their axonal arborization is distinct from pyramidal cells because they form intricate networks of branches that wrap around the soma of nearby pyramidal cells, forming a characteristic “basket” structure. Their axonal structure is generally more isotropic compared to pyramidal cells or SBCs and LBCs. Nevertheless, the thresholds for tangential electric fields are about 14% higher compared to normal electric fields (θ=0° and θ=180°). On average, the thresholds of L4 NBCs are 2% and 23% higher compared to L2/3 PCs and L5 PCs, respectively. The results for biphasic excitation are shown in Fig S6 in the Supplementary Material.

3.4.3 Large basket cells

The threshold results of L4 LBCs for monophasic excitation are shown in Fig. S7. Compared to PCs, LBCs exhibit a high degree of collateralization in their axonal tree. They can have multiple branches and collaterals that extend in different directions within the same cortical layer or across layers. A distinct directional sensitivity of the thresholds can be again observed together with an asymmetric modulation when the electric field changes along the somato-dendritic axis. On average, the thresholds of L4 LBCs are 9% lower than L2/3 PCs and 11% higher compared to L5 PC, respectively. Of all the basket cells investigated, the LBCs have the lowest thresholds. The average thresholds of LBCs are 24% and 10% lower compared to SBCs and NBCs, respectively. The results for biphasic excitation are shown in Fig. S8 in the Supplementary Material.

3.5 Sensitivity analysis

We used a 15th-order approximation to construct the surrogate models of the threshold maps using pygpc (Weise, Poßner, et al., 2020). The normalized root mean square deviation between the original model and the gPC approximation is 0.32% for L5 PC under monophasic excitation derived from 105 random samples. The accuracies of the gPC approximations of the L2/3 and L4 cells are very similar and given in the repository by Weise, Worbs, et al. (2023). The results of the sensitivity analysis of the threshold map of L5 PC with monophasic excitation are shown in Fig. 7. It can be seen that the surrogate model (Fig. S11b) almost perfectly resembles the behavior of the original model (Fig. 7a). The absolute differences between both is shown in Fig. 7c. The probability density distribution of the electric field threshold is shown in Fig. 7d under the assumption that the parameters θ and Δ|E˜| are beta distributed as in case of the realistic head model simulations (see Fig. 3 for parameter values). It can be observed that the distribution is u-shaped and bimodal because of the cyclic behavior of the electric field threshold over the polar angle θ. The results for L5 in case of monophasic excitation as well as for L2/3 PC and L4 S/N/LBC for both monophasic and biphasic excitation are given in the repository (Weise, Worbs, et al., 2023).

Fig. 7.

Results of the sensitivity analysis of the electric field threshold map of L5 PCs with monophasic excitation: (a) Original model of the threshold map of L5 PC with monophasic excitation; (b) gPC approximation (surrogate) of the original model; (c) Absolute difference between the original model and the gPC approximation; (d) Probability density of the electric field threshold for the original model and the gPC approximation using N=105 samples under the assumption that θ and Δ|E˜| are beta distributed (see Fig. 3 for parameters).

Fig. 7.

Results of the sensitivity analysis of the electric field threshold map of L5 PCs with monophasic excitation: (a) Original model of the threshold map of L5 PC with monophasic excitation; (b) gPC approximation (surrogate) of the original model; (c) Absolute difference between the original model and the gPC approximation; (d) Probability density of the electric field threshold for the original model and the gPC approximation using N=105 samples under the assumption that θ and Δ|E˜| are beta distributed (see Fig. 3 for parameters).

Close modal

The Sobol indices, that is, the fractions of the total variance, which originate from θ, Δ|E˜|, and the combination of both are given in Table 1. The polar angle θ has the strongest influence on the stimulation behavior for all investigated cell types. In contrast, the Sobol indices of Δ|E˜| are much lower, ranging between 2-5%, but the parameter significantly contributes to the increase of the accuracy of the overall model. There is even an exception in the L2/3 cells under biphasic excitation, where the influence reaches almost 25%. The relatively high influence of Δ|E˜| in this case can be also observed in the threshold map in Fig. S2. However, the analysis of the absolute Sobol coefficients has shown that the influence is similar to the monophasic case. But since the directional dependence in biphasic pulses with respect to θ is lower compared to monophasic pulses, the relative contribution of the field gradient is increased. Naturally, larger cells like L2/3 PC and L5 PC are affected more by Δ|E˜| measured in %/mm than smaller cells. The fact that this effect is observed only in the L2/3 PCs and not in the L5 PCs may be related to the deactivation of the main axons of the L5 PCs following the model of Aberra et al. (2020). This was motivated to account for the truncation of the main axons in the slicing process because the truncated terminals in these neurons were unrealistically close to their cell bodies (Aberra et al., 2020).

3.6 Verification

The average threshold model is compared against reference simulations using a high-resolution realistic head model. For the application of the average model, we extracted the electric field parameters θ and Δ|E˜| in every cortical element in the ROI on layer 2/3, 4, and 5 and determined the electric field thresholds (Fig. 5, Figs.S1-S9) by linearly interpolating the data between the sampling points. The approach is computationally very efficient, as the computation time is only a fraction of the one needed for the electric field computation. In the reference simulations, we calculated the stimulation thresholds for every neuron at every cortical location in the ROI separately by coupling the actual electric fields from the realistic head model into every neural compartment. Finally, we averaged the thresholds over all neurons and assigned the resulting average threshold to the ROI element. The resulting electric field threshold maps between the average threshold model and the reference simulations are shown in Fig. 8 for all cell types under investigation. The results for biphasic excitation are shown in Figure S12 in the Supplementary Material.

Fig. 8.

Comparison of electric field threshold maps (in V/m) for monophasic excitation determined using the average model and the reference model: The rows show the electric field threshold maps (in V/m) of the L2/3 PC, L4 S/N/LBC and the L5 PC between the average model (first column) and the reference model (second column). The last two columns show the absolute and relative difference between the models. The underlying electric field distribution and field direction is shown in Fig. 4. The results for biphasic excitation are shown in Fig. S12.

Fig. 8.

Comparison of electric field threshold maps (in V/m) for monophasic excitation determined using the average model and the reference model: The rows show the electric field threshold maps (in V/m) of the L2/3 PC, L4 S/N/LBC and the L5 PC between the average model (first column) and the reference model (second column). The last two columns show the absolute and relative difference between the models. The underlying electric field distribution and field direction is shown in Fig. 4. The results for biphasic excitation are shown in Fig. S12.

Close modal

For all stimulation conditions, the two models agree very well. The highest relative differences can be observed mainly at the gyral rims and the sulcal walls. Comparing the distributions and signs of the relative differences between monophasic and biphasic waveforms, it can be observed that they slightly depend on the stimulation waveform and the resulting current direction.

Additionally, we calculated the stimulation threshold maps when the TMS coil is located over the M1 region with a 45° orientation towards the fissura longitudinalis. For this, we determined the ratio between the electric field threshold map from Fig. 9 and the corresponding electric field distribution of this particular coil position, which was calculated assuming a normalized stimulation strength of 1 A/µs. This results in a map of the stimulation strength of the TMS stimulator in A/µs needed to stimulate the neurons with this particular coil position. Again, a high agreement between the average threshold model and the reference model can be observed. Note that the relative difference distributions in the last column of Fig. 9 are the same as for the electric field threshold maps from Fig. 8 since the electric field distribution is cut out in the error calculation. The analogous results for biphasic excitation are shown in Fig. S13.

Fig. 9.

Comparison of stimulation intensity threshold maps (in A/µs) for monophasic excitation determined using the average model and the reference model: The first two rows show the stimulation threshold maps (in A/µs) of the L2/3 PC and the last two rows of the L5 PC between the average model (first column) and the reference model (second column). The last two columns show the absolute and relative difference between the models. It is assumed that the TMS coil is located over the M1 area with an orientation of 45° towards the fissura longitudinalis. The maps indicate the stimulation strength of the TMS device in A/µs, which is required to stimulate this cortical area for this particular coil position and orientation. The underlying electric field distribution and field direction is shown in Fig. 4. The results for biphasic excitation are shown in Fig. S13.

Fig. 9.

Comparison of stimulation intensity threshold maps (in A/µs) for monophasic excitation determined using the average model and the reference model: The first two rows show the stimulation threshold maps (in A/µs) of the L2/3 PC and the last two rows of the L5 PC between the average model (first column) and the reference model (second column). The last two columns show the absolute and relative difference between the models. It is assumed that the TMS coil is located over the M1 area with an orientation of 45° towards the fissura longitudinalis. The maps indicate the stimulation strength of the TMS device in A/µs, which is required to stimulate this cortical area for this particular coil position and orientation. The underlying electric field distribution and field direction is shown in Fig. 4. The results for biphasic excitation are shown in Fig. S13.

Close modal

The histograms of the relative differences are shown in Fig. 10 together with the different error measures. The 1st and 99th percentile of the errors are given in Table 2 for monophasic and biphasic TMS pulses together with the average height and the radial extent of the cells.

Fig. 10.

Differences of the threshold maps between the average model and the reference model for monophasic excitation. Histograms of the relative difference between the reference model and the average threshold model over the ROI elements. Normalized root mean square error (NRMSE), mean absolute percentage error (MAPE), and coefficient of determination (R2) for L2/3 PC, L4 S/N/LBC, and L5 PC with monophasic excitation. The results for biphasic excitation are shown in Fig. S14.

Fig. 10.

Differences of the threshold maps between the average model and the reference model for monophasic excitation. Histograms of the relative difference between the reference model and the average threshold model over the ROI elements. Normalized root mean square error (NRMSE), mean absolute percentage error (MAPE), and coefficient of determination (R2) for L2/3 PC, L4 S/N/LBC, and L5 PC with monophasic excitation. The results for biphasic excitation are shown in Fig. S14.

Close modal
Table 2.

1st and 99th percentiles of the relative differences between the average and the reference model. Additionally, the average height of the cells is reported together with the average radial extent.

Cell typeL2/3 PCL4 SBCL4 NBCL4 LBCL5 PC
Wave-formMono-phasicBi-phasicMono-phasicBi-phasicMono-phasicBi-phasicMono-phasicBi-phasicMono-phasicBi-phasic
1st perc. (in %) -9.543 -3.217 -1.820 -1.206 -1.561 -0.872 -2.180 -1.168 -6.096 -2.384 
99th perc. (in %) 4.121 4.378 1.909 1.545 4.101 3.695 3.320 2.856 6.109 4.988 
Avg. height (in μm) 1518 (SD = 455) 741 (SD = 211) 783 (SD = 174) 1091 (SD = 246) 1844 (SD = 257) 
Avg. rad. extent (in μm) 865 (SD = 408) 392 (SD = 100) 511 (SD = 132) 607 (SD = 161) 1052 (SD = 269) 
Cell typeL2/3 PCL4 SBCL4 NBCL4 LBCL5 PC
Wave-formMono-phasicBi-phasicMono-phasicBi-phasicMono-phasicBi-phasicMono-phasicBi-phasicMono-phasicBi-phasic
1st perc. (in %) -9.543 -3.217 -1.820 -1.206 -1.561 -0.872 -2.180 -1.168 -6.096 -2.384 
99th perc. (in %) 4.121 4.378 1.909 1.545 4.101 3.695 3.320 2.856 6.109 4.988 
Avg. height (in μm) 1518 (SD = 455) 741 (SD = 211) 783 (SD = 174) 1091 (SD = 246) 1844 (SD = 257) 
Avg. rad. extent (in μm) 865 (SD = 408) 392 (SD = 100) 511 (SD = 132) 607 (SD = 161) 1052 (SD = 269) 

The distribution of relative differences is relatively symmetric and the means are close to zero. The remaining variance results from the inhomogeneity of the electric field across the neurons. This effect is the stronger the larger the cell is and is therefore particularly evident in L2/3 PCs and L5 PCs. Systematic field distortions in a particular direction across neurons, such as those occurring at the gyral crowns, are neglected and result in deviations from the exact reference model because in the simplified model, only the decay of the electric field across the somato-dendritic axis can be accounted for due to averaging over the azimuthal angle.

3.7 Comparison to directional sensitivity of motor evoked potentials

The predicted orientation sensitivity of the neurons is compared to the orientation sensitivity of MEPs. The polar plot in Fig. 11 shows the MEP amplitudes for different electric field angles and stimulation intensities together with the predictions of the recruitment rate from the average threshold model from L5 PC under monophasic excitation. To make both representations comparable, the data were normalized to their respective maximum values. The average threshold model closely matches the orientation sensitivity of MEPs, and the NRMSE between the experimental data and model predictions is 8.5%. It can be clearly observed that the directional sensitivity is more pronounced for low stimulation intensities close to rMT than for higher ones. The cortical column cosine model resembles the general behavior of the directional sensitivity of the MEPs at low stimulation intensities, but cannot represent the stimulation behavior in the transition to higher stimulation intensities. The ball-and-stick model was also not able to reflect the direction sensitivity over the investigated parameter range for both the incident angle θ and the stimulation intensity.

Fig. 11.

Comparison between directional sensitivity of motor evoked potentials and the derived recruitment rate from the theoretical neuronal response model: The plots show the directional sensitivity of measured MEPs as black crosses (exp) at different stimulation intensities with respect to resting motor threshold of subject 16 from Souza et al. (2022). The solid lines (sim) show the corresponding trajectory of the recruitment rate along θ assuming a constant electric field along the somato-dendritic axis (Δ|E˜|=0%/mm), and the dashed and the dotted lines show the predictions of the directional sensitivity of the MEPs according to the cortical column cosine model (cos) and the ball-and-stick model (b&s). The MEPs were normalized to their maximum values for comparability. The NRMSE between the experimental data and the recruitment rate model (sim) is 8.5%.

Fig. 11.

Comparison between directional sensitivity of motor evoked potentials and the derived recruitment rate from the theoretical neuronal response model: The plots show the directional sensitivity of measured MEPs as black crosses (exp) at different stimulation intensities with respect to resting motor threshold of subject 16 from Souza et al. (2022). The solid lines (sim) show the corresponding trajectory of the recruitment rate along θ assuming a constant electric field along the somato-dendritic axis (Δ|E˜|=0%/mm), and the dashed and the dotted lines show the predictions of the directional sensitivity of the MEPs according to the cortical column cosine model (cos) and the ball-and-stick model (b&s). The MEPs were normalized to their maximum values for comparability. The NRMSE between the experimental data and the recruitment rate model (sim) is 8.5%.

Close modal

In order to link the predicted electric field to actual neural activation, a range of different proposals with varying degrees of complexity have been put forward. The simplest approach just considers the magnitude of the electric field as a proxy for the activation strength (e.g., Weise, Numssen, et al. (2023), without any dependency on direction or local gradient of the field. This method disregards experimental observations and theoretical considerations showing that the activation threshold does indeed depend on the incidence angle between the field direction and the orientation of the axons (Ranck, 1975; Rudin & Eisenman, 1954; Rushton, 1927). This consideration led to the proposal of the cortical column cosine model (Fox et al., 2004), which is based on the assumption that axons aligned with the somato-dendritic axis (i.e., perpendicular to the cortical surface) dominate the stimulation process, and therefore predicts that only the projection of the electric field onto that axis has an effect. As a consequence, purely tangential fields would lead to no stimulation.

However, it is known that the axonal arbors of cortical cells are much branched and cover all directions (Aberra et al., 2020). In line with this, in earlier TMS motor mapping experiments (Numssen et al., 2021; Weise, Numssen, et al. (2020); Weise, Numssen, et al., 2023), we could show that the tangential field component does indeed have a substantial predictive power towards the resulting motor evoked potentials. In fact, it was even considerably more powerful than the radial component (i.e., the one aligned with the somato-dendritic axis), which can be understood in the light of the cortical geometry: At the gyral crowns, the field is largely tangential to the cortical surface (thus, having less impact), but its magnitude is much larger due to the greater proximity to the coil, thus overcompensating the former effect. In contrast, at the sulcal walls, located more distant to the coil, the field is normal, but much weaker, and therefore often does not effectively stimulate.

In order to obtain a more accurate account of the coupling between the electric field and the activation states of cortical neurons, detailed biological models based on realistic neuronal geometry and realistic Hodgkin-Huxley-like neural dynamics have been proposed (Aberra et al., 2018, 2020). These models have the potential to deliver a detailed and accurate picture of neuronal activation by TMS. However, they are computationally extremely demanding and therefore hardly suitable for routine applications, such as mapping or dosing. Moreover, the utilized neural geometries must be considered as samples of a distribution and do not account for any precise individual cortical architecture. This suggests that the predictions of these models should be representable in low-parametric models without much loss.

In our study, we attempt to bridge the gap between, on the one hand, the imprecise oversimplification of the magnitude, cortical column cosine, and ball-and-stick models and, on the other hand, the unwieldy and time-consuming biologically realistic modeling. The model we propose is as easily applicable as the former, while it very closely approximates the predictions of the latter. It determines the stimulation thresholds as functions of field angle with respect to the somato-dendritic axis, intensity, pulse waveform, and field decay along the somato-dendritic axis, and only requires the induced electric field as an input variable. Comparison with reference simulations with a detailed neuronal model yielded normalized root mean square errors of only 1.5-2.5%. It should be emphasized, however, that our model is not independent, but depends, for initial calibration, on a biologically realistic model based on the principal approach of Aberra et al. (2018, 2020).

The performance of the current model in terms of the observed differences to the reference simulations is very comparable to the convolutional neural network (CNN) model of Aberra et al. (2023). The MAPE scores of the average threshold model presented here and the CNN are 1.35% and 1.4% for L2/3 PCs, 0.961% and 1.1% for L4 LBCs, and 1.8% and 1.4% for L5 PCs, respectively. Both models are very efficient and reasonably fast. Nonetheless, the gPC approximation of the parameterized average threshold models is currently the most efficient description form, as it consists only of simple polynomial evaluations and a single matrix-vector multiplication (eq. (6)) to determine the thresholds for multiple stimulation locations. For example, the calculation of 105 thresholds for L5 PC takes 70 ms using a 2.80 GHz i7-7600U laptop CPU compared to about 1 min for the CNNs (Aberra et al., 2023).

Since our priority was to develop accurate descriptors of the thresholds for whole populations of certain cell types, we had to substantially increase the number of cells to reduce the confidence intervals of the means. For example, the mean 95th percentile confidence interval of the average threshold models of L5 PCs for monophasic excitation could be decreased from 85.5 to 29.2 V/m by using 30 instead of 5 cells.

The main advantage of the current model is its preservation of parameterization concerning the most important geometrical electric field parameters. This preservation allows for an easy transformation between neuron thresholds and electric field topographies. This way, one can assess the effective electric field, as opposed to relying on simpler proxies such as magnitude or normal components, both before and after conducting experiments within a specific target region.

Our model predicts a certain dependence of the stimulation threshold on the angle of incidence of the electric field, which is more pronounced for monophasic pulses than for biphasic pulses. However, compared to the predictions from simplified models such as the cortical column cosine model or ball-and-stick neurons, this dependence is much more moderate. This allows tangential electric fields of sufficient strength to contribute to the stimulation, as has been observed in previous experimental studies (Numssen et al., 2021; Weise, Numssen, et al., 2020). In particular, the L2/3 PC require 108%, L4 SBC require 113%, L4 NBC require 110%, L4 LBC require 114%, and L5 PC require 120% of the longitudinal stimulation strength (θ=0°) at θ=90° for monophasic excitation, respectively. For a comparison of the directional sensitivity profiles of our model and the cortical column cosine and ball-and-stick models, see Fig. 5, Figs. S1-S9.

These findings are confirmed by a comparison to the experimentally observed orientation sensitivity of MEPs (Fig. 11), where a difference of only 8.5% was observed between model predictions and experimental data. We observed a high directional sensitivity at low stimulation intensities close to the motor threshold, while at higher stimulation intensities the directional sensitivity rapidly decreases. While our model appears to be already a quite good predictor of the directional sensitivity observed by Souza et al. (2022), there are also deviations and limitations, which have to be addressed for a stringent validation. This is mainly explainable in the light of some important discrepancies between the assumptions underlying our model and those made by these authors. First, the results of Souza et al. are based on electric fields predicted using a spherical head model, while our model works with a realistic head model. Second, Souza’s report is based on the electric field direction with respect to the global coordinate system, while our angle definition is local and changes across the strongly curved cortical surface. Third, the location of the neuronal populations that mediate the relationship between stimulation and MEP is only roughly known in Souza et al. (2022). It may therefore be that the field angles at that location are different from those predicted for the assumed target spot. Accordingly, for an even better comparison, the currents in the multi-coil array would have to be optimized subject- and target-specifically to realize an ideal rotation of the electric field at a constant field strength in the target. This in turn requires precise knowledge of the target and thus a prior mapping of the motor cortex such as in Weise, Numssen, et al. (2023).

The major advantage of the presented model is its simplicity without sacrifice of realism. The availability of look-up tables of threshold maps and recruitment rates allows for the simple construction of interpolators and functions for computation. Alternatively, polynomial-based surrogate models based on generalized polynomial chaos (gPC; Weise, Poßner, et al., 2020) can be used for this purpose and provide high accuracy. Examples are given in the repository of Weise, Worbs, et al. (2023).

Importantly, the model is easy to adapt and refine, if more or better information about the neuronal geometry of particular tissues becomes available, using the provided scripts and simulation code (https://github.com/TorgeW/TMS-Neuro-Sim). Already in this study, we were able to distinguish between the stimulation thresholds and distributions among various distinct cell types. We observed that L5 PCs had the lowest thresholds compared to all other cell types studied, followed by L4 LBC and L2/3 PCs, which had 10% and 22% higher thresholds, respectively. This “library” of cellular stimulation profiles may be extended in the future. By comparison with experimentally observed stimulation profiles, such cell-specific sensitivity profiles may potentially allow for testing hypotheses about which cells are actually stimulated in particular experimental situations.

These traits allow for efficient implementation and extension of TMS models in the context of optimization, mapping, and dosing without the need to implement time-consuming and complicated neuron models. Especially in the field of cognitive TMS experiments, where an adequate dosing strategy is still subject to research, the gained knowledge could significantly not only contribute to the identification of the effectively stimulated regions but also exclude regions that are not eligible for stimulation due to the underlying electric field distribution and orientation relative to the cortex.

The threshold maps have revealed interesting parameter combinations of θ and Δ|E˜| that enable particularly effective stimulation. Here, it is an interesting observation that an increase of the electric field along the somato-dendritic axis of the neurons (Δ|E˜|>0) from the GM surface to the WM surface is as likely as a decrease and is usually in the range of ±20%/mm (Fig. 3b). Future optimization studies could be directed towards identifying coil positions and orientations that realize these parameter combinations in the targeted region. As a result, such new optimization strategies would have great potential to significantly enhance the overall efficacy of TMS and reduce the required dose. At the same time, the optimization criterion can be extended such that the electric field is oriented to prevent stimulation of other brain regions by targeting particularly high thresholds. This principled approach of multi-objective optimization was also taken up by Lueckel et al. (2022) in the framework of electric field and connectivity optimized TMS targeting.

Another area of application for the presented models is in the extension of existing mapping procedures (Numssen et al., 2021; Weise, Numssen, et al., 2020; Weise, Numssen, et al., 2023), as mentioned previously. Instead of the electric field magnitude, some kind of effective electric field could be used as a regressor for localization. It is noted that the integration of the stimulation thresholds into the analysis procedures occurs solely at the modeling level, thus improving the efficiency of the mapping procedures without increasing the experimental effort. Stronger even, based on the fact that we have an estimate of the stimulation threshold at every cortical location, we can successively exclude locations which are stimulated without an observable effect (e.g., MEP), and thereby even decrease the experimental effort.

4.1 Limitations of the study

The number of neurons available was limited. We were able to significantly expand the original dataset of Aberra et al. (2018, 2020), but especially for the calculation of the recruitment rate, a higher number of neurons would increase the model accuracy. This is especially true for incident angles where threshold variances are high.

Moreover, we limited the analysis to pyramidal cells in L2/3, L5, and basket cells in L4, which take a major part in generalized cortical circuits (Di Lazzaro et al., 2012). However, it is known that other cell types, like spiny stellates in L4, also play a major role in the stimulation of cortical microcircuits. The development of average threshold models for other cell types is straightforward using the tools provided in the repository Weise, Worbs, et al. (2023) and the Python package TMS-Neuro-Sim (https://github.com/TorgeW/TMS-Neuro-Sim) if the appropriate morphologies and parameterizations are available.

A limitation of the current model is the absence of neuronal activity of the cells. Shirinpour et al. (2021) showed that temporal characteristics of neuronal activity and synaptic inputs play a major role to estimate neural responses especially in the case of rTMS. The development of a statistical average threshold model taking into account neuronal activity requires a statistical description of this activity itself and should be investigated in the future.

In the modeling, we also neglected the effect of the presence of the neurons and other cells on the external electric field. While for the macroscopic field estimation, these structures are already accounted for through the (macroscopically acquired) tissue conductivity, at the microscale, the presence of low conducting membranes might cause local deviations from that macroscopically predicted field, which may have an effect on the actual stimulation of neurons. Following Aberra et al. (2018, 2020), another potential source for the general overestimation of thresholds may be related to the sodium channel models of the axons.

4.2 Future work

Insights into the stimulation behavior of neurons are essential to develop realistic coupling models for downstream neuronal mass models along the lines of Montbrió et al. (2015) or Jansen and Rit (1995), which in turn could be used to model the dynamic processes of entire populations of neurons, such as the D- and I-wave dynamics in the motor cortex (Di Lazzaro et al., 1998, 2012; Ziemann, 2020).

In further follow-up studies, the degree to which the spatial fine-structure of the electric field is affected by the high membrane resistance of the neurons should be investigated. The resulting change in the electric field distribution may have a non-negligible influence on the local electric field angles and magnitudes, which in turn change the stimulation thresholds. However, this will require very detailed volume conductor models of whole cortical columns or at least geometric information about the neuron surfaces and their position with respect to each other. It is expected that this type of model will require high computational power to solve and is far from being routinely used in daily TMS experiments and that it will lead to an anisotropic macroscopic conductivity profile as well as a potentially modified sensitivity profile due to local electric field fluctuations. The former can be estimated with diffusion-weighted MRI (Güllmar et al., 2010), but its influence on the stimulation behavior on a micro- and mesoscopic scale is yet unknown. The goal of such a study could be the derivation of a new generation of low-parametric models, in a similar sense as in this study, in order to be able to apply the gained knowledge in practice.

All data underlying the results presented in this paper, together with additional details including the average threshold models, the recruitment rate operators, and the neuron compartment models, are publically available in a repository (Weise, Worbs, et al., 2023), where we provide look-up tables, interpolators, and polynomial approximations for further use (https://doi.org/10.17605/OSF.IO/C8J35).

The code and associated example scripts are published in the open-source Python package TMS-Neuro-Sim (https://github.com/TorgeW/TMS-Neuro-Sim).

K.W., T.W. developed the theory, conceived and planned the numerical experiments. T.W. wrote the code and carried out the threshold simulations. W.V.G., A.T.J. curated and contributed the neuron data from the Blue Brain Project. T.W. post-processed the data and conducted the verification study. V.H.S. contributed the experimental data. K.W. determined the recruitment rate operators, performed the sensitivity analysis, and compared the results to the experimental findings of V.H.S., K.W., T.W., A.T.J., T.R.K., V.H.S., W.V.G., and A.T.J. contributed to the interpretation of the results. K.W., A.T., and T.R.K. supervised the project. K.W. took the lead in writing the manuscript. K.W., T.R.K. obtained funding for this project. All authors provided critical feedback and helped to shape the research, analysis, and manuscript.

This work was partially supported by the German Science Foundation (DFG) (grant number WE 59851/2 to K.W. KN 588/10-1 to T.R.K.), Lundbeckfonden (grant no. R244-2017-196 and R313-2019-622) and the NVIDIA Corporation (donation of one Titan Xp graphics card to K.W.). V.H.S. is funded by the Academy of Finland (decision No 349985). A.J. and W.V.G are supported by funding to the Blue Brain Project, a research center of the École polytechnique fédérale de Lausanne (EPFL), from the Swiss government’s ETH Board of the Swiss Federal Institutes of Technology.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00036.

Aberra
,
A. S.
,
Lopez
,
A.
,
Grill
,
W. M.
, &
Peterchev
,
A. V.
(
2023
).
Rapid estimation of cortical neuron activation thresholds by transcranial magnetic stimulation using convolutional neural networks
.
NeuroImage
,
275
,
120184
. https://doi.org/10.1016/j.neuroimage.2023.120184
Aberra
,
A. S.
,
Peterchev
,
A. V.
, &
Grill
,
W. M.
(
2018
).
Biophysically realistic neuron models for simulation of cortical stimulation
.
Journal of Neural Engineering
,
15
(
6
),
066023
. https://doi.org/10.1088/1741-2552/aadbb1
Aberra
,
A. S.
,
Wang
,
B.
,
Grill
,
W. M.
, &
Peterchev
,
A. V.
(
2020
).
Simulation of transcranial magnetic stimulation in head model with morphologically-realistic cortical neurons
.
Brain Stimulation
,
13
(
1
),
175
189
. https://doi.org/10.1016/j.brs.2019.10.002
Anwar
,
H.
,
Riachi
,
I.
,
Hill
,
S.
,
Schurmann
,
F.
, &
Markram
,
H.
(
2009
).
An approach to capturing neuron morphological diversity
. In
De Schutter
E.
(Ed.),
Computational neuroscience: Realistic modeling for experimentalists
(pp.
211
232
).
MIT Press Scholarship
Online. https://doi.org/10.7551/mitpress/9780262013277.003.0010
Aspart
,
F.
,
Ladenbauer
,
J.
, &
Obermayer
,
K.
(
2016
).
Extending integrate-and-fire model neurons to account for the effects of weak electric fields and input filtering mediated by the dendrite
.
PLoS Computational Biology
,
12
(
11
),
e1005206
. https://doi.org/10.1371/journal.pcbi.1005206
Bédard
,
C.
, &
Destexhe
,
A.
(
2008
).
A modified cable formalism for modeling neuronal membranes at high frequencies
.
Biophysical Journal
,
94
(
4
),
1133
1143
. https://doi.org/10.1529/biophysj.107.113571
Carnevale
,
N. T.
, &
Hines
,
M. L.
(
2006
).
The NEURON book
.
Cambridge University Press
. https://doi.org/10.1017/cbo9780511541612
Dale
,
A. M.
,
Fischl
,
B.
, &
Sereno
,
M. I.
(
1999
).
Cortical surface-based analysis: I. Segmentation and surface reconstruction
.
NeuroImage
,
9
(
2
),
179
194
. https://doi.org/10.1006/nimg.1998.0395
Di Lazzaro
,
V.
,
Profice
,
P.
,
Ranieri
,
F.
,
Capone
,
F.
,
Dileone
,
M.
,
Oliviero
,
A.
, &
Pilato
,
F.
(
2012
).
I-wave origin and modulation
.
Brain Stimulation
,
5
(
4
),
512
525
. https://doi.org/10.1016/j.brs.2011.07.008
Di Lazzaro
,
V. D.
,
Restuccia
,
D.
,
Oliviero
,
A.
,
Profice
,
P.
,
Ferrara
,
L.
,
Insola
,
A.
,
Mazzone
,
P.
,
Tonali
,
P.
, &
Rothwell
,
J.
(
1998
).
Effects of voluntary contraction on descending volleys evoked by transcranial stimulation in conscious humans
.
The Journal of Physiology
,
508
(
2
),
625
633
. https://doi.org/10.1111/j.1469-7793.1998.625bq.x
Fischl
,
B.
,
Dale
,
A. M.
,
Sereno
,
M. I.
,
Tootell
,
R. B. H.
, &
Rosen
,
B. R.
(
1998
).
A coordinate system for the cortical surface
.
NeuroImage
,
7
(
4
),
S740
. https://doi.org/10.1016/s1053-8119(18)31573-8
Fox
,
P. T.
,
Narayana
,
S.
,
Tandon
,
N.
,
Sandoval
,
H.
,
Fox
,
S. P.
,
Kochunov
,
P.
, &
Lancaster
,
J. L.
(
2004
).
Column‐based model of electric field excitation of cerebral cortex
.
Human Brain Mapping
,
22
(
1
),
1
14
. https://doi.org/10.1002/hbm.20006
García‐Cabezas
,
M. Á.
, &
Barbas
,
H.
(
2014
).
Area 4 has layer IV in adult primates
.
European Journal of Neuroscience
,
39
(
11
),
1824
1834
. https://doi.org/10.1111/ejn.12585
Gaser
,
C.
,
Dahnke
,
R.
,
Thompson
,
P. M.
,
Kurth
,
F.
, &
Luders
,
E.
(
2022
).
CAT-a computational anatomy toolbox for the analysis of structural MRI data
.
bioRxiv
. https://doi.org/10.1101/2022.06.11.495736
Güllmar
,
D.
,
Haueisen
,
J.
, &
Reichenbach
,
J. R.
(
2010
).
Influence of anisotropic electrical conductivity in white matter tissue on the EEG/MEG forward and inverse solution. A high-resolution whole head simulation study
.
NeuroImage
,
51
(
1
),
145
163
. https://doi.org/10.1016/j.neuroimage.2010.02.014
Hess
,
A.
, &
Young
,
J. Z.
(
1949
).
Correlation of internodal length and fibre diameter in the central nervous system
.
Nature
,
164
,
490
491
. https://doi.org/10.1038/164490a0
Hursh
,
J. B.
(
1939
).
Conduction velocity and diameter of nerve fibers
.
American Journal of Physiology
,
127
,
131
139
. https://doi.org/10.1152/ajplegacy.1939.127.1.131
Jansen
,
B. H.
, &
Rit
,
V. G.
(
1995
).
Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns
.
Biological Cybernetics
,
73
(
4
),
357
366
. https://doi.org/10.1007/bf00199471
Lorensen
,
W. E.
, &
Cline
,
H. E.
(
1987
).
Marching cubes: A high resolution 3D surface construction algorithm
.
ACM SIGGRAPH Computer Graphics
,
21
(
4
),
163
169
. https://doi.org/10.1145/37402.37422
Lueckel
,
M.
,
Radetz
,
A.
,
Yuen
,
K.
,
Mueller-Dahlhaus
,
F.
,
Kalisch
,
R.
, &
Bergmann
,
T. O.
(
2022
).
TU-176. E-field-and connectivity-optimized TMS targeting: A pilot TMS-fMRI validation at the single-subject level
.
Clinical Neurophysiology
,
141
,
S30
S31
. https://doi.org/10.1016/j.clinph.2022.07.080
Markram
,
H.
,
Muller
,
E.
,
Ramaswamy
,
S.
,
Reimann
,
M. W.
,
Abdellah
,
M.
,
Sanchez
,
C. A.
,
Ailamaki
,
A.
,
Alonso-Nanclares
,
L.
,
Antille
,
N.
,
Arsever
,
S.
,
Kahou
,
G. A. A.
,
Berger
,
T. K.
,
Bilgili
,
A.
,
Buncic
,
N.
,
Chalimourda
,
A.
,
Chindemi
,
G.
,
Courcol
,
J.-D.
,
Delalondre
,
F.
,
Delattre
,
V.
,…
Schürmann
,
F.
(
2015
).
Reconstruction and simulation of neocortical microcircuitry
.
Cell
,
163
(
2
),
456
492
. https://doi.org/10.1016/j.cell.2015.09.029
Montbrió
,
E.
,
Pazó
,
D.
, &
Roxin
,
A.
(
2015
).
Macroscopic description for networks of spiking neurons
.
Physical Review X
,
5
(
2
),
021028
. https://doi.org/10.1103/physrevx.5.021028
Nielsen
,
J. D.
,
Madsen
,
K. H.
,
Puonti
,
O.
,
Siebner
,
H. R.
,
Bauer
,
C.
,
Madsen
,
C. G.
,
Saturnino
,
G. B.
, &
Thielscher
,
A.
(
2018
).
Automatic skull segmentation from MR images for realistic volume conductor models of the head: Assessment of the state-of-the-art
.
NeuroImage
,
174
,
587
598
. https://doi.org/10.1016/j.neuroimage.2018.03.001
Numssen
,
O.
,
Zier
,
A. L.
,
Thielscher
,
A.
,
Hartwigsen
,
G.
,
Knösche
,
T. R.
, &
Weise
,
K.
(
2021
).
Efficient high-resolution TMS mapping of the human motor cortex by nonlinear regression
.
NeuroImage
,
245
,
118654
. https://doi.org/10.1016/j.neuroimage.2021.118654
Penny
,
W. D.
,
Friston
,
K. J.
,
Ashburner
,
J. T.
,
Kiebel
,
S. J.
, &
Nichols
,
T. E.
(Eds.). (
2011
).
Statistical parametric mapping: The analysis of functional brain images
.
Elsevier
. https://doi.org/10.1016/b978-012372560-8/50024-3
Puonti
,
O.
,
Van Leemput
,
K.
,
Saturnino
,
G. B.
,
Siebner
,
H. R.
,
Madsen
,
K. H.
, &
Thielscher
,
A.
(
2020
).
Accurate and robust whole-head segmentation from magnetic resonance images for individualized head modeling
.
Neuroimage
,
219
, 117044. https://doi.org/10.1016/j.neuroimage.2020.117044
Ramaswamy
,
S.
,
Courcol
,
J. D.
,
Abdellah
,
M.
,
Adaszewski
,
S. R.
,
Antille
,
N.
,
Arsever
,
S.
,
Atenekeng
,
G.
,
Bilgili
,
A.
,
Brukau
,
Y.
,
Chalimourda
,
A.
,
Chindemi
,
G.
,
Delalondre
,
F.
,
Dumusc
,
R.
,
Eilemann
,
S.
,
Gevaert
,
M. E.
,
Gleeson
,
P.
,
Graham
,
J. W.
,
Hernando
,
J. B.
,
Kanari
,
L.
,… &
Markram
,
H.
(
2015
).
The neocortical microcircuit collaboration portal: A resource for rat somatosensory cortex
.
Frontiers in Neural Circuits
,
9
,
44
. https://doi.org/10.3389/fncir.2015.00044
Ranck
Jr,
J. B.
(
1975
).
Which elements are excited in electrical stimulation of mammalian central nervous system: A review
.
Brain Research
,
98
(
3
),
417
440
. https://doi.org/10.1016/0006-8993(75)90364-9
Rudin
,
D. O.
, &
Eisenman
,
G.
(
1954
).
The action potential of spinal axons in vitro
.
The Journal of General Physiology
,
37
(
4
),
505
538
. https://doi.org/10.1085/jgp.37.4.505
Rushton
,
W. A. H.
(
1927
).
The effect upon the threshold for nervous excitation of the length of nerve exposed, and the angle between current and nerve
.
The Journal of Physiology
,
63
(
4
),
357
. https://doi.org/10.1113/jphysiol.1927.sp002409
Saturnino
,
G. B.
,
Madsen
,
K. H.
, &
Thielscher
,
A.
(
2019
).
Electric field simulations for transcranial brain stimulation using FEM: An efficient implementation and error analysis
.
Journal of Neural Engineering
,
16
(
6
),
066032
. https://doi.org/10.1088/1741-2552/ab41ba
Shirinpour
,
S.
,
Hananeia
,
N.
,
Rosado
,
J.
,
Tran
,
H.
,
Galanis
,
C.
,
Vlachos
,
A.
,
Jedlicka
,
P.
,
Queisser
,
G.
, &
Opitz
,
A.
(
2021
).
Multi-scale modeling toolbox for single neuron and subcellular activity under transcranial magnetic stimulation
.
Brain Stimulation
,
14
(
6
),
1470
1482
. https://doi.org/10.1016/j.brs.2021.09.004
Souza
,
V. H.
,
Nieminen
,
J. O.
,
Tugin
,
S.
,
Koponen
,
L. M.
,
Baffa
,
O.
, &
Ilmoniemi
,
R. J.
(
2022
).
TMS with fast and accurate electronic control: Measuring the orientation sensitivity of corticomotor pathways
.
Brain Stimulation
,
15
(
2
),
306
315
. https://doi.org/10.1016/j.brs.2022.01.009
Thielscher
,
A.
,
Antunes
,
A.
, &
Saturnino
,
G. B.
(
2015
),
Field modeling for transcranial magnetic stimulation: A useful tool to understand the physiological effects of TMS?
In
IEEE EMBS 2015
,
Milano, Italy
.
IEEE
. https://doi.org/10.1109/embc.2015.7318340
Virtanen
,
P.
,
Gommers
,
R.
,
Oliphant
,
T. E.
,
Haberland
,
M.
,
Reddy
,
T.
,
Cournapeau
,
D.
,
Burovski
,
E.
,
Peterson
,
P.
,
Weckesser
,
W.
,
Bright
,
J.
,
van der Walt
,
S. J.
,
Brett
,
M.
,
Wilson
,
J.
,
Millman
,
K. J.
,
Mayorov
,
N.
,
Nelson
,
A. R. J.
,
Jones
,
E.
,
Kern
,
R.
,
Larson
,
E.
,… &
Van Mulbregt
,
P.
(
2020
).
SciPy 1.0: Fundamental algorithms for scientific computing in Python
.
Nature Methods
,
17
(
3
),
261
272
. https://doi.org/10.1038/s41592-020-0772-5
Wagner
,
T.
,
Gangitano
,
M.
,
Romero
,
R.
,
Théoret
,
H.
,
Kobayashi
,
M.
,
Anschel
,
D.
,
Ives
,
J.
,
Cuffin
,
N.
,
Schomer
,
D.
,… &
Pascual-Leone
,
A.
(
2004
).
Intracranial measurement of current densities induced by transcranial magnetic stimulation in the human brain
.
Neuroscience Letters
,
354
(
2
),
91
94
. https://doi.org/10.1016/s0304-3940(03)00861-9
Wang
,
B.
,
Grill
,
W. M.
, &
Peterchev
,
A. V.
(
2018
).
Coupling magnetically induced electric fields to neurons: Longitudinal and transverse activation
.
Biophysical Journal
,
115
,
95
107
. https://doi.org/10.1016/j.bpj.2018.06.004
Waxman
,
S. G.
, &
Kocsis
,
J. D.
(
1995
).
The axon: Structure, function, and pathophysiology
.
Oxford University Press
. https://doi.org/10.1056/nejm199512143332422
Weise
,
K.
,
Numssen
,
O.
,
Kalloch
,
B.
,
Zier
,
A. L.
,
Thielscher
,
A.
,
Haueisen
,
J.
,
Hartwigsen
,
G.
, &
Knösche
,
T. R.
(
2023
).
Precise motor mapping with transcranial magnetic stimulation
.
Nature Protocols
,
18
(
2
),
293
318
. https://doi.org/10.1038/s41596-022-00776-6
Weise
,
K.
,
Numssen
,
O.
,
Thielscher
,
A.
,
Hartwigsen
,
G.
, &
Knösche
,
T. R.
(
2020
).
A novel approach to localize cortical TMS effects
.
NeuroImage
,
209
,
116486
. https://doi.org/10.1016/j.neuroimage.2019.116486
Weise
,
K.
,
Poßner
,
L.
,
Müller
,
E.
,
Gast
,
R.
, &
Knösche
,
T. R.
(
2020
).
Pygpc: A sensitivity and uncertainty analysis toolbox for Python
.
SoftwareX
,
11
,
100450
. https://doi.org/10.1016/j.softx.2020.100450
Weise
,
K.
,
Worbs
,
T.
,
Kalloch
,
B.
,
Souza
,
V. H.
,
Jaquier
,
A. T.
,
Van Geit
,
W.
,
Thielscher
,
A.
, &
Knösche
,
T. R.
(
2023
).
Direction sensitivity of cortical neurons towards TMS induced electric fields—Data and code. osf.io
. https://doi.org/10.17605/OSF.IO/C8J35
Ziemann
,
U.
(
2020
).
I-waves in motor cortex revisited
.
Experimental Brain Research
,
238
(
7
),
1601
1610
. https://doi.org/10.1007/s00221-020-05764-4
Zienkiewicz
,
O. C.
, &
Zhu
,
J.
(
1992
).
The superconvergent patch recovery (SPR) and adaptive finite element refinement
.
Computer Methods in Applied Mechanics and Engineering
,
101
(
1–3
),
207
224
. https://doi.org/10.1016/0045-7825(92)90023-d

Author notes

*

Contributed equally

This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International (CC BY 4.0) license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. For a full description of the license, please visit https://creativecommons.org/licenses/by/4.0/legalcode.

Supplementary data