Directional sensitivity of cortical neurons towards TMS-induced electric fields

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.

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.

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$⁠:

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

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.

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$⁠.

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

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:

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 ⁰. 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.

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.

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:

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.

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∙10⁶ nodes and ~29∙10⁶ 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).

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):

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:

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

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

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.

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.

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.

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.

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.

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.

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.

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.

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 10⁵ 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).

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).

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.

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.

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.

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.

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.

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 10⁵ 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.

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.

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.