Directly monitoring target cells’ activities under stimulation would be immensely valuable for the interpretation of TMS effects, but only a few in vitro studies exist so far ^{19, 20} . However, computational studies can explore the effects of the electromagnetic fields on neural activation in the intact brain by simulating models of neural stimulation in silico. In early computational models, straight axonal fibers were considered numerically and the response of neurons induced by the external field was modeled by means of the cable equation ^{21, 22} . Later models investigated the role of cell morphology using multi-compartmental modeling ^{23– 25} . The responses of cortical neurons vary depending on not only the neuronal morphology but also orientation relative to the induced electric field and stimulus amplitude ^{6, 26, 27} . Therefore, recent models have used anatomical information on cortical folding that induces a wide range of field orientations and have applied the calculated electric field from the head model to the neuronal models ^{28, 29} .

Here, we use an advanced multi-scale modeling approach that combines a high-resolution head model with detailed multi-compartmental neuron models. We construct an anatomically realistic head model based on MRI and calculate the external currents that affect neurons via the TMS-induced electric field with high accuracy. We concentrate on the hand knob area of the motor cortex that is the predominant target of many TMS studies ^{16, 30} . A multitude of layer 5 and layer 3 pyramidal neurons (L5/L3 PNs) is incorporated on the basis that they might be primary activators of the corticospinal tract and provide the main input to the direct pathway ^{24, 28, 31, 32} . We estimate the target area of activation as a function of coil orientation as well as the stimulation intensities required to activate different types of pyramidal neurons. Finally, we predict the precise sites where the neurons initiate their action potentials. We find that action potentials are mostly generated at the axon initial segment of these cells.

In order to study the cellular effects of TMS in the brain we employed a multi-scale computational modeling approach combining a volume conductor head model with detailed neuronal models of cortical pyramidal neurons. The motor cortex, especially the hand area, was considered as a cortical target location. The volume conductor head model was used to simulate the stimulus-induced electric fields; it was based on the SimNIBS v1.1 software pipeline ^{14, 33} . The precise impact of these fields on different neural targets was evaluated using multi-compartmental models of pyramidal neurons embedded into the head model. This allowed us to predict differences in individual neurons’ susceptibilities to TMS depending on neuron placement and coil orientation.

Figure 3 depicts the magnitude of the electric fields (| E → |, top row) and the orthogonal component of the electric fields ( E ⊥ , bottom row) for different coil orientations. All calculations were performed for a rate of change of the coil current of 1 A/μs. Electric fields had higher magnitudes in the precentral and postcentral gyrus and focused on the top of the gyri, regardless of coil orientations. We observed only slight changes in the absolute field strengths depending on coil orientation. In contrast, the orthogonal component of electric fields ( E ⊥ ) showed different spatial patterns compared to the electric field magnitude. High strengths of E ⊥ were found on the walls of the gyri and strongly depended on coil orientation. Furthermore, while the spatial extent of | E → | was the same for the standard orientation and the +180 degree orientation, the sign of E ⊥ in the +180 degree orientation was reversed due to the reversed sign of the induced electric fields. Interestingly, the maximum value of | E → | depended on coil orientation; it was lowest in the base coil orientation and highest at +135 degrees. However, the maximum values of E ⊥ were highest for the standard coil orientation and lowest at +90 degrees.

To assess the neuronal activations as a function of coil orientation, we determined the excitation threshold required to cause action potentials of L5/L3 PNs. For each coil orientation, we kept increasing the stimulator output until a neuron generated an action potential or we reached a maximum rate of change of the current defined as 171 A/μs. Our focus is on the excitability for a stimulation intensity corresponding to 67 A/μs, as this value corresponds to the average motor threshold for the Magstim 200 stimulator connected to the coil modeled ^{18, 29, 49, 51} . The excitability of L5/L3 PNs was predicted either by the direct estimation of the electric field ( E ⊥ map in the Figure 4(c)) or by simulating the induced depolarization and firing of the detailed neuronal models (threshold maps in Figure 4(d,e)). The color of the threshold maps represents the stimulator output necessary to activate the corresponding cell and the estimated excitable area in the E ⊥ maps. The blue colored areas indicate an excitability in the opposite direction, because the head model was linear with respect to the electric field. As shown in Figure 4(a), we virtually divided the precentral and postcentral gyrus to better visualize the excitability in the walls of the gyrus.

In L5 PNs, we observed that the predicted excitability depended on coil orientation for both E ⊥ ( Figure 4(b)) and threshold maps ( Figure 4(c)). For the base orientation and +45 degrees, a high excitability was predominantly observed in the wall of the precentral gyrus. In contrast, for orientations +90 degrees to +225 degrees we observed high excitability in the wall of the postcentral gyrus. Comparing these threshold maps to the E ⊥ maps, we see that the excited regions were quite well matched to the results from E ⊥ . Furthermore, in the standard direction, the spatial extents where L5 PNs were activated for stimulation intensities corresponding to the motor threshold seemed to enlarge with increasing coil rotation, while for the opposite direction of the coil current highly excitable areas shrank with increasing coil rotations.

Overall, the excitability of L3 PNs showed behavior comparable to that of L5 PNs ( Figure 4 (c,d)), but notable differences in threshold maps between L5 and L3 PNs were as follows: while L5 PNs in the top of the gyri were never activated, L3 PNs were excited in the top and also the wall of gyri. The excitable areas of L3 PNs caused by the +90 and +135 degree stimulations were relatively focused on the upper parts of the wall of the postcentral gyrus, while L5 PNs placed in the deeper parts of the sulcus were activated. Furthermore in L3 PNs, the excitability in the precentral gyrus and the postcentral gyrus was comparable and a bigger area was affected than for L5 PNs. The discrepancies between L5 and L3 PNs confirmed that the morphology and placement of neuronal models has an important impact beyond the position relative to the coil.

The percentage of excited neurons for a stimulation intensity at the motor threshold is shown in Figure 5. We separately analyzed neurons falling in the precentral and postcentral gyrus. When we focused on PN activations in the precentral gyrus, the highest percentage of excited neurons was observed at +90 degrees rather than the base orientation. For the base orientation PNs in the sulcal wall along the central sulcus were activated. At +90 degrees PNs were activated mostly in the opposite sulcal wall ( Supplementary Figure S3). For the postcentral gyrus, the maximum percentage of activated PNs was observed when the coil was oriented at +135 degrees, which is in agreement with the threshold maps in Figure 4.

Our multiscale modeling approach also allowed us to identify where on each cell an action potential was initiated. The majority of action potentials were initiated at the axon initial segment and others at the axon near the boundary between GM and WM for L5 PNs and at the middle and terminal points for L3 PNs ( Figure 6). In the base orientation, threshold stimulation elicited action potentials first at the initial segment for 90% of both the L5 and the L3 PNs. This fraction increased with increasing coil rotations up to 97% at +135 degrees in L5 PNs and up to 95% at +90 degrees in L3 PNs. Example plots of membrane potential dynamics induced by the threshold stimulus evoking action potentials are shown in Figure 7. We observe the propagation of the action potentials from the axon initial segment to the more distal parts of the neurons. In both L5 and L3 PNs, following the action potential at the initial segment, the soma was activated next, as it is closest to the initial segment. The terminal points of the axons were activated last as they are most distal from the axon initial segment. Since the axons of the simulated L5 PNs are quite long compared to that of the L3 PNs, the arrival of the action potential at the terminal point was substantially delayed. Similarly, dendrites of L5 PNs showed delayed activation while in the L3 PNs dendrites were occasionally activated early.

The PNs that were morphologically reconstructed had asymmetric dendritic trees that might affect the neuronal responses. We studied the impact of dendritic trees on threshold maps and the percentage of excited neurons for a stimulation intensity at the motor threshold by rotating them in steps of 90 degrees around the axis defined by their apical dendrite for a fixed base coil orientation. In the threshold maps for L5 PNs as shown in Figure 8(b), the highest variations of the thresholds caused by these rotations were observed in the boundary between the top of the precentral gyrus and the sulcus. However, the dependence on coil orientation in the threshold maps did not change and thus the L5 PNs in the precentral gyrus were activated consistently. Compared to the percentage of excited L5 PNs with randomly rotated dendritic trees (10.53% as shown in Figure 5), the fixed orientation of dendritic trees induced changes in the fraction of activated neurons of up to 2%. The threshold variations in the L3 PNs were hardly noticeable compared to those of the L5 PNs. The percentage of excited L3 PNs was 20.85% with randomly rotated dendritic trees, and when the orientations of dendritic trees were fixed it resulted in changes of at most 0.3%. In addition, the impact of scaling of PNs was investigated by reducing the dimension of the dendritic trees by 10%, 20% or 30%, as shown in Figure 8(b). For this, we scaled the length of all dendritic compartments while keeping their diameters identical ⁴⁶ . We observed consistently activated sites in the threshold maps and only slight changes of the percentage of excited L5 PNs of up to 0.3%. Thus, we found that the rotations of the dendritic trees had a bigger impact on PN excitability than that of scaling the dendritic trees. Overall, morphological changes in dendritic trees did not alter the spatial extent of activated regions much.

The detailed mechanisms through which TMS activates cortical cells and cortical circuits are still not fully understood. In this study, we used multi-scale computational modeling to predict cortical activation as a function of coil orientation in two different ways. First, we simply considered the strength of the component of the TMS-induced electric field that is orthogonal to the gray matter surface as suggested by the C ³-model ^{10, 41} . Second we developed a detailed computational modeling approach that combined an anatomically realistic head model with complex multi-compartment neuronal models of L5/L3 PNs and quantified their stimulation thresholds. A first major finding was the characterization of the induced electric fields and the thresholds of L5/L3 PNs as a function of coil orientations as shown in Figure 4. In addition, threshold variations according to different morphologies of PNs were observed.

The magnitude of the electric field was considered first, because the strength of the electric field is commonly used as a simple way to predict areas of neuronal activation ^{6, 7, 52} . We found that the magnitude of the TMS-induced field is focused on the top of the gyrus, which is in agreement with previous modeling studies ^{10, 12– 14} . However, the electric field magnitude showed little dependency on coil orientation ¹⁰ . Then, we investigated the directional electric field, specifically the orthogonal component that is perpendicular to the cortical surface, as this has been suggested to contribute most to the TMS-induced activation according to the C ³-model ^{10, 41, 53} . We found a strong dependence of the orthogonal field component on coil orientation, as shown in Figure 3. In contrast to the electric field magnitude, the highest orthogonal field values were found in the sulcal walls and never on the apex (or crown) of the gyrus.

While the analysis of TMS-induced electric fields has been widely addressed in the past, the incorporation of multi-compartment neuronal models has hardly been investigated. To permit a more detailed understanding of the biophysical mechanisms of TMS, a few previous modeling studies employed detailed neuronal models and calculated the membrane potential dynamics generated by the electromagnetic field. However, these attempts had various limitations. First, in early studies no anatomical information on large-scale brain morphology was applied ^{23– 25} . Rather than constructing a finite element head model, these studies applied a uniform electric field to the neuronal model. Even though such investigations achieved reasonable results regarding the neuronal activation, they did not consider the effects of the complex folding patterns of the cortex and the effects of tissue borders such as the borders between GM on the one hand and CSF or WM on the other hand. However, the importance of anatomically realistic head models has been shown convincingly ^{12– 14} . Furthermore, the impact of detailed brain anatomy has been considered in various methods of brain stimulation and substantial differences have been demonstrated by improving anatomical information related to the head model ^{54– 57} . Salvador et al. (2011) investigated neuronal responses using a simplified head model of a cortical sulcus with several types of neurons and found changes of the stimulation threshold depending on the pulse waveform and the coil orientation ²⁹ . However, the used head model had an approximated geometry restricted to the motor cortex and a full complex geometry, such as the hook-shaped hand knob, was not considered. Furthermore, the modeled coil orientation was limited to anterior to posterior and its reversal due to the simplified geometry of the head model. Most recently, Goodwin and Butson (2015) proposed a more realistic approach that integrates an anatomically realistic head model derived from MR images with detailed neuronal models ²⁸ . They considered the excitability of neurons as a function of coil orientation. However, in contrast to our results, excitability maps hardly showed a systematic dependence on coil orientation and activation thresholds were lower in the gyral crown. We speculate that this might be caused by the different morphology of PNs or the different way in which they calculated the external currents to simulate neuronal responses. Also, we considered two types of L3/L5 PNs spread over a wider region of the cortex. Finally, as another major result, we also established the site of action potential initiation and found that most PNs are activated at the axon initial segment and action potential initiations at other parts of the neuron are comparatively rare.

The threshold maps we calculated demonstrate acute sensitivity to coil orientation, but different spatial extents of activation were observed according to the different morphologies of the PNs. In L5 PNs, activation thresholds were low in the sulcal walls, matching predictions based on the orthogonal component of the electric field. The excitation in the sulcal cortical surface was consistent with the well-established columnar neuronal orientation and functional organization of the cortex and functional imaging studies ^{41, 58} . The L3 PNs had different morphology with shorter axons than the L5 PNs such that they were located completely within the gray matter. Similar to the L5 PNs, the coil orientation had a significant impact on the responses of the L3 PNs, but the precise patterns of the threshold maps differed between the L3 and L5 PNs. As L3 PNs in the gyral crown and the upper parts of the sulcal wall were predominantly activated, as would be expected based on their proximity to the stimulation coil ⁵⁹ .

The neural response to TMS is composed of a direct (D) and several indirect (I) waves. The D-wave is thought to be produced by direct activation of L5 PNs as we have modeled it here and is followed by I-waves that are thought to be generated by synaptic excitation and/or re-excitation of L5 pyramidal cells with longer latencies ^{59, 60} , presumably via pyramidal cells in superficial cortical layers L2 and L3. According to Di Lazzaro et al. (2004), at the lowest stimulation intensity to evoke neuronal responses, an I-wave is elicited, and with increasing stimulation intensity, the earlier, small D-wave is produced ³⁰ . This indicates that thresholds for eliciting I-waves are lower than those for eliciting D-waves ⁶¹ . In this work, we explored the excitation thresholds of both L3 and L5 PNs and found as another major result that the percentage of excited neurons for a stimulation intensity at the motor threshold was about two times higher for L3 PNs than for L5 PNs. Furthermore, the activation of the L3 PNs was consistently higher than that of the L5 PNs for the full range of stimulation intensities. Mainen and Sejnowski (1996) compared the dendritic structure of L3 and L5 PNs with a common distribution of ion channels and found that a smaller intracellular current injection was necessary to activate L3 PNs compared to L5 PNs ⁴² . Thus, the morphology of L3 PNs might result in their higher excitability in response to TMS. In addition, the lower stimulation thresholds of the L3 PNs are consistent with lower stimulation intensities required to produce I-waves ^{50, 62} , and the higher stimulation intensities required to produce D-waves.

For motor cortex, there is evidence that the base coil orientation induces the highest motor evoked potential ^{63, 64} , while we found that the +90 degrees orientation activates the highest number of PNs across the whole part of the precentral gyrus falling within the ROI ( Figure 5(b)). At +90 degrees PNs were activated mostly in the curved part of the precentral gyrus ( Supplementary Figure S2), and according to the spatial distributions of thresholds ( Figure 4) the sulcal wall along the central sulcus was targeted better at the base coil orientation than for other coil orientations. Thus, the base coil orientation is recommended for M1 by our model, congruent with previous research ^{63, 64} , and +135 degrees should be ideal to stimulate the postcentral gyrus.

The question of the precise initiation site of action potentials is a central issue for understanding the physiological effects of TMS. According to our study, the dominant initiation site leading to action potentials is the axon initial segment in both L5 and L3 PNs. In previous modeling studies, the action potentials were initiated at the axons crossing the boundary between gray matter and white matter, where the conductivity changes abruptly ⁴⁹ and at the bending parts of the axon due to charge accumulation ²⁹ . However Goodwin and Butson (2015) combined a realistic head model with detailed PN models and observed that most action potentials were initiated at or close to the axon hillock just adjacent to the axon initial segment ²⁸ . Pashut et al. (2011) have also argued for action potential initiation at the axon initial segment ²⁴ . Furthermore, this is consistent with previous studies arguing that action potentials giving rise to the D-wave might be initiated close to the soma and/or axon initial segment ^{65, 66} . In our study, L5 PN action potentials were only rarely initiated at the axon where it crosses the boundary between gray matter and white matter, where tissue conductivity changes abruptly ⁴⁹ . It will be important to verify these results with more realistic neuron (in particular: axon) models.

There are several limitations in our modeling study. A first limitation is that we have assumed isotropic conductivity, as is common in computational studies. Opitz et al. (2011) revealed that anisotropy might create hot spots in the WM with increased field strength that might affect neural excitation ¹² , and Seo et al. (2015) reported that anisotropy affected L5 PNs significantly while it had only minor impact on L3 PNs ⁴⁵ . Thus anisotropic conductivity might have significant effects on L5 PN axons running through the WM.

In this work, the L3 and L5 PNs were taken from cat visual cortex due to the lack of models for most human cortical cell types. Thus, while we lengthened PNs to fit the cortex, the uncertainties regarding morphology of neurons was not fully studied. Wu et al. (2016) incorporated a multitude of PNs with various stimuli ²⁵ and Salvador et al. (2010) constructed various types of neural structures including pyramidal neurons, interneurons, and association fibers ²⁹ ; they found that the excitability can be shaped by field orientation, pulse wave form, and diameter of neurons. In addition, changes in the electrical properties, such as membrane properties and ion channels had the largest influence on neuronal excitability ²⁵ . However, despite the uncertainty with regard to properties of PNs, we produced results matching both experimental studies and other computational studies ^{23, 25} that incorporated the same models of PNs.

While we observed the stimulation of neural activity in the superficial cortex close to the coil, TMS might also affect deep brain areas that cannot be stimulated directly. This can be explained on the basis of the propagation of action potentials along white matter fiber tracts. Recent studies modeled tractography-based white matter fiber tracts using diffusion tensor imaging (DTI) and observed activation of axon bundles ^{9, 11, 12} . Compared to fiber tracts in previous modeling, we modeled straightly stretched axons of L5 PNs inside the WM. Due to such limitations, the axons inside the WM occasionally passed through protruding parts of GM. Notwithstanding the fact that this intersection could affect the neuronal responses such as the action potential initiation or thresholds, most PNs initiated action potentials at the axon initial segment and coil orientation dependency observed in threshold maps was consistent with observations in previous studies. Further developments in tractography may improve detailed neuronal models and may lead to a deeper understanding of the TMS-induced brain activity propagations from the superficial cortex to distant brain regions.

Another limitation of our model is that the reconstructed PNs were synaptically isolated. For the L5 cells this means that we basically studied the generation mechanism of D-waves. The activation of L3 cells could be seen as a proxy for the generation of I-waves. A logical next step is to synaptically couple L3 and L5 cells as done in a recent model for the generation of D and I-waves using L5 PNs that were contacted by a pool of excitatory and inhibitory layer 2 and 3 neurons ⁵⁰ . This model successfully reproduced various characteristics of I-waves and highlighted the importance of the complex morphology of the L5 PNs for the generation of I-waves. An improvement would be to use the anatomical information on the activation of PNs as modeled here, as we found a clear difference in the threshold maps for L5 and L3 PNs based on their morphology. Therefore, in future work, we plan to incorporate synaptic connections between L3 and L5 PNs. We hope that this will bring us one step closer to a detailed understanding of the mechanisms through which TMS activates cortical circuits, paving the way for more precise and effective application of TMS based on individual brain morphology in clinical and basic research settings.

The data referenced by this article are under copyright with the following copyright statement: Copyright: � 2017 Seo H et al.

Data associated with the article are available under the terms of the Creative Commons Zero "No rights reserved" data waiver (CC0 1.0 Public domain dedication).

F1000Research: Dataset 1. Figure 3 input data, 10.5256/f1000research.9277.d132069 ⁶⁷

F1000Research: Dataset 2. Figure 4 input data, 10.5256/f1000research.9277.d132070 ⁶⁸

F1000Research: Dataset 3. Figure 5 raw data, 10.5256/f1000research.9277.d132071 ⁶⁹

F1000Research: Dataset 4. Figure 6 raw data, 10.5256/f1000research.9277.d132072 ⁷⁰