A multi-scale computational model of the effects of TMS on motor cortex

The detailed biophysical mechanisms through which transcranial magnetic stimulation (TMS) activates cortical circuits are still not fully understood. Here we present a multi-scale computational model to describe and explain the activation of different pyramidal cell types in motor cortex due to TMS. Our model determines precise electric fields based on an individual head model derived from magnetic resonance imaging and calculates how these electric fields activate morphologically detailed models of different neuron types. We predict neural activation patterns for different coil orientations consistent with experimental findings. Beyond this, our model allows us to calculate activation thresholds for individual neurons and precise initiation sites of individual action potentials on the neurons’ complex morphologies. Specifically, our model predicts that cortical layer 3 pyramidal neurons are generally easier to stimulate than layer 5 pyramidal neurons, thereby explaining the lower stimulation thresholds observed for I-waves compared to D-waves. It also shows differences in the regions of activated cortical layer 5 and layer 3 pyramidal cells depending on coil orientation. Finally, it predicts that under standard stimulation conditions, action potentials are mostly generated at the axon initial segment of cortical pyramidal cells, with a much less important activation site being the part of a layer 5 pyramidal cell axon where it crosses the boundary between grey matter and white matter. In conclusion, our computational model offers a detailed account of the mechanisms through which TMS activates different cortical pyramidal cell types, paving the way for more targeted application of TMS based on individual brain morphology in clinical and basic research settings.

This article is included in the gateway. Transcranial magnetic stimulation (TMS) is a neurostimulation and neuromodulation technique that noninvasively activates neurons in the brain 1,2 . It generates a time varying magnetic field using a coil above the head, which induces an electric field in the brain that can be of sufficient magnitude to depolarize neurons beyond their firing threshold. In recent years, TMS has been widely tested as a tool for diagnosis and treatment for a broad range of neurological and psychiatric disorders [3][4][5] . Although the efficacy of TMS has been demonstrated, there remains great uncertainty regarding the factors influencing the affected brain areas and relevant circuits.
To provide a better understanding of the biophysical mechanisms behind TMS, several computational studies have been performed to try to reveal the effects of a number of parameters contributing to variable outcomes. The majority of models try to predict the brain regions influenced by TMS exclusively based on stimulus-induced electric fields 6,7 . While early studies utilized spherical models of the human head, in recent years high-resolution volume conduction models of the head have been developed from human magnetic resonance imaging (MRI) to improve the accuracy of calculated electric fields [8][9][10][11][12][13][14][15] . These models have revealed that the geometry of the volume conduction model, such as complex gyral folding patterns, is one of the key parameters determining the induced electric fields. In addition, computational studies have been extended by connecting numerical results with experimental observations to show the agreement of computed electric fields with physiological measurements [16][17][18] .
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][24][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.

Methods
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 stimulusinduced 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.

Volume conductor model
The simulated effects of TMS depend not only on stimulation parameters but also on the anatomical information specified in the volume conductor model. To calculate the precise electric field, a volume conductor head model for TMS that reflected T1-weighted and T2-weighted magnetic resonance (MR) images was constructed using the example dataset provided by SimNIBS v1. 1

Field calculations
The TMS-induced electric field was calculated based on the SimNIBS v1.1 pipeline. Briefly, the electric field, primary electric field was directly determined by the coil geometry and the secondary electric field was caused by charge accumulations at tissue interfaces. Using magnetic dipoles to model the TMS coil, the primary electric field was calculated directly without the volume conductor model and then used as input for the secondary electric field calculation via a finite element method using the GetFEM++ library and MATLAB 33,40 . The Magstim 70 mm figure-8 coil was represented by magnetic dipoles positioned above the hand knob ( Figure 1) and the stimulator output was set to 1 A/μs. The base coil orientation was defined relative to the direction of the central sulcus such that the electric field induced was in the posterior to anterior direction (the yellow arrow in Figure 1(c)). Then, three additional coil orientations were tested by rotating in steps of 45 degrees and reversed orientations were simulated by changing the sign of the current through the coil.
To investigate the TMS-induced cellular effects, we quantified the magnitude of the electric field | | E  and the orthogonal component of the electric field E E n ⊥ = ⋅   to the gray matter surface, where n  is the normal vector for the boundary surface element. The orthogonal component was expected to contribute to TMS-induced brain activation by the theoretical cortical column cosine model of TMS efficacy (C 3 -model) 10,41 .

Multi-compartmental neuronal models
We adapted existing multi-compartmental models of layer 5 and 3 pyramidal neurons (L5/L3 PNs) from cat visual cortex 42 using the NEURON simulation software 43 . The electrical properties were unchanged from the original models. Briefly, a high density of fast, inactivating voltage-dependent Na + channels were present in the axon hillock and axon initial segment, and a low density of these channels was present in the soma and dendrites. Slow Ca 2+ -dependent K + channels and high threshold Ca 2+ channels were located in the soma and dendrites. Except for the dendrites, fast K + channels were present.
L5/L3 PNs were combined virtually with the head model and modified to accommodate the irregular geometry of the cortex 28,29,44-47 , as shown in Figure 2. To reduce superfluous computations, we preselected a region of interest (ROI) of 50 × 50 × 50 mm 3 around the hand knob and then placed L5/L3 PNs in each triangular element comprising the gray matter surface. The multi-compartmental models of PNs consisted of a series of compartments connected by resistors. Each compartment was further discretized into segments of equal length to allow for accurate numerical simulation. The center points of each segment were extracted and used to calculate the necessary changes to neuron geometry, as described below. The dendritic trees were lengthened or shortened by re-scaling the lengths of the compartments according to the local dimensions of the cortex such that dendrites reached layer 1 and the orientation was perpendicular to the cortical surface 46,48 . Since the morphology of the dendritic trees was not symmetric and it might influence the neuronal activation, L5/L3 PNs had randomly rotated dendritic trees at different locations. The axons of L5 PNs were defined to curve beyond the boundary between GM and WM in the direction of the internal capsule (Supplementary Figure S1) as described previously 45 . Further adjustments of L5 PNs' geometry were as follows (illustrated in Supplementary Figure S2): each dendritic tree was oriented such that its principal axis would align with the normal vector of its associated triangular surface element. The bending part of the axon was calculated according to the normal vector of the surface element. The arc length of the axon bend was set to 0.6 mm when the z-component of the normal vector was positive and otherwise the arc length was 0.3 mm (compare Supplementary Figure S2). Note that when we varied the angle and arc length of the axon bend, it usually did not alter the activation threshold. The axons of L3 PNs were defined to terminate in layer 5/6 within the GM. Altogether, a total of 10,888 L5 PNs and 10,888 L3 PNs was constructed. This process was implemented in MATLAB (MathWorks, Natick, MA, USA).

Computation of neuronal activation induced by stimulation
The membrane potentials induced by stimulation were approximated by adding an external current source I ext to the cable model 2,21,22,24,25 : 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    Sites include the axon initial segment (iseg) and the boundary between gray matter and white matter (boundary). Additionally, the terminal part (terminal) and middle point of the axon (middle) for L3 PNs were considered. Most action potentials are first evoked at the axon initial segment of L5 PNs (96.31±1.72%) and L3 PNs (92.76±2.42%). The remaining number of L5 PNs show action potential initiation at the axon near the boundary between gray matter and white matter. Only few L5 PNs (0.49±0.14%) initiate action potentials simultaneously at the axon initial segment and the GM-WM boundary. For L3 PNs, middle (1.05±0.78%) and terminal points (2.08±1.45%) of axons are also activated occasionally.
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 46 . 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.

Discussion
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 3 -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 10 . 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 3 -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 29 . 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 28 . 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 59 .
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 Iwave is elicited, and with increasing stimulation intensity, the earlier, small D-wave is produced 30 . This indicates that thresholds for eliciting I-waves are lower than those for eliciting D-waves 61 . 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 42 . 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 49 and at the bending parts of the axon due to charge accumulation 29 . 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 28 . Pashut et al. (2011) have also argued for action potential initiation at the axon initial segment 24 . 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 . The presence of a highly increased density of fast voltage-gated Na + and K + channels leads to high intrinsic excitability of the axon initial segment 67 . Additionally, the displaced charges from across the entire dendritic tree can accumulate at the axon initial segment due to the neuronal morphology. This means that the initiation of action potentials in the initial segments is not necessarily a direct effect of TMS at that particular site, but rather the consequence of membrane properties and the morphology of PNs. 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 49 . 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. 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. 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 25 . 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.
An important limitation of our study is that we have not included realistic axon morphologies. A thorough modeling of the complex morphologies of local axons could reveal additional sites of activation. In particular, axon terminals whose orientations are aligned with the electric field may be particularly easy to stimulate.
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 50 . 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.

Competing interests
No competing interests were disclosed.   The revised manuscript was substantially improved. I have two remaining points, which I suggest to clarify in form of a minor revision of the article:

Grant information
At the TBS conference in Göttingen 2016, Peterchev and colleagues presented modelling work (P90 in the abstract handbook) related to this study, and found the axon terminals to have the lowest thresholds. While their work is not yet published as full paper, it is worth adding a statement to the discussion of the study limitations here, e.g. that axon terminals of incoming axonal projections were not modelled, but might be further neural elements that are well excitable by TMS.
I am still puzzled by the way in which the effect of the TMS field on the neural compartments was modelled (page 5 of the revised manuscript). This seems to depend on the assumption of rotationally symmetric compartments of negligible diameter (as the change of E is only evaluated along the compartment, and the transverse component of E is neglected). In contrast, the soma has a more complex geometry with several branching points at which the axons and the dendritic tree start. As the validity of the main conclusions (i.e. that the axon initial segment close to the 2. tree start. As the validity of the main conclusions (i.e. that the axon initial segment close to the soma is most easily stimulated by TMS) depends on an accurate modelling of the TMS effects on this region of the neuron, I would ask the authors to elaborate further on this question. A visualization of the compartmentalization might help to clarify this issue.
No competing interests were disclosed.

Competing Interests:
I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
Author Response 28 Apr 2017 , Gwangju Institute of Science and Technology, Korea, South Hyeon Seo 1. Thank you for your careful suggestion. We added the statement to the discussion of the study limitations, as follows:

An important limitation of our study is that we have not included realistic axon morphologies. A thorough modeling of the complex morphologies of local axons could reveal additional sites of activation.
In particular, axon terminals whose orientations are aligned with the electric field may be particularly easy to stimulate.
2. The compartments are indeed all modeled as cylinders, but with finite diameters. Importantly, that the site of action potential initiation is the axon initial segment does not imply that the electric field necessarily has a strong direct effect on this compartment. There are two effects contributing to action potential initiation at the axon initial segment. First, in response to TMS displaced charges from across the entire dendritic tree can accumulate at the axon initial segment due to the neuronal morphology. Second, due its membrane properties, the axon initial segment is easily excitable. We have clarified these aspects in the discussion, as follows: 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 (2011) have also argued for et al. 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 . The presence of a highly increased density of fast voltage-gated Na and K channels leads to high intrinsic excitability of the axon initial segment . Additionally, the displaced charges from across the entire dendritic tree can accumulate at the axon initial segment due to the neuronal morphology. This means that the initiation of action potentials in the initial segments is not necessarily a direct effect of TMS at that particular site, but In our study, L5 PN rather the consequence of membrane properties and the morphology of PNs. 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. 1.

4.
verify these results with more realistic neuron (in particular: axon) models.
No competing interests were disclosed. Competing Interests: 24  Please provide a more detailed description in the Methods of how you obtained the primary and secondary electric fields. For instance, you state that the former is directly determined from the coil and head geometries, whilst the latter is calculated using FEM. Why this distinction? Cannot FEM with appropriate boundary conditions be used for a single/combined electric field calculation?
Please clarify more clearly the nature of the magnetic dipoles you used to simulate the TMS coil i.e. how many and how were they distributed on the scalp? These dipoles are shown in the panel of Figure 1b as yellow arrows, but I find the figure very unclear, with minimal description in the text.
Can you speculate on how the assumption of isotropic conductivity in WM would affect your results? Current is more likely to flow along WM fiber tracts than across these, potentially altering your electric fields and activation results.
In the abstract, replace "corctial" with "cortical", and your field calculation methods, replace "in anterior to posterior" with "in anterior to posterior". the No competing interests were disclosed. Competing Interests: 1.

3.
No competing interests were disclosed. Competing Interests:

I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.
Author Response 31 Jan 2017 , Gwangju Institute of Science and Technology, Korea, South Hyeon Seo The contribution of the proposed approach is combining TMS-induced electric field with the pyramidal neurons. All the process to calculate stimulus-induced electric field using the head model was based on the SimNibs software pipeline, and thus we just briefly introduced the procedure related to field calculations. According to the SimNibs pipeline, they used magnetic dipoles to model the TMS coils as it allows to easily determine the magnetic vector potential of a coil (A-field; the primary electric fields). Furthermore, magnetic vector potential is unaffected by the conductivity in the head model and thus it is straight-forward to calculate. The secondary electric field arises from the charge accumulations at conductivity discontinuities and thus was numerically determined using FEM. For clarity, we now wrote that the field calculation followed the SimNibs pipeline and rephrased corresponding sentences as follows: "The TMS-induced electric field was calculated based on the SimNibs v1.1 pipeline. Briefly, the electric field, E→=−∂A→∂t−∇ϕ=−Ep→−Es,→ consisted of primary (Ep→) and secondary (Es→) electric fields. The primary electric field was directly determined by the coil geometry and the secondary electric field caused by charge accumulations at tissue interfaces. Using magnetic dipoles to model the TMS coil, the primary electric field was calculated directly without the volume conductor model and then used as input for the secondary electric field calculation via a finite element method using the GetFEM++ library and MATLAB ." As already responded to the previous comment, we simulated the field following SimNibs software and thus the magnetic dipoles representing the TMS coil was provided in SimNibs. In this paper, the magstim 70 mm figure-8 coil was modelled as two circular disks (radius r = 5 cm) which are divided into 10 rings each (Thielscher and Kammer, 2004.) In previous study, Opitz Because of this, one could speculate that the anisotropic conductivity may increase activation of L5 PNs. We rephrased these points in the discussion section as follows: "There are several limitations in our modeling study. A first limitation is that we have assumed isotropic conductivity, as is common in computational studies. Seo and coworkers combine realistic calculations of the electric field that is induced by TMS with multi-compartmental neural models in order to reveal how and at which sites TMS generates neural activity. Only a few prior studies have targeted this topic in detail, and the presented study contributes relevant new insights into the putative stimulation mechanisms of TMS. Clearly, on the long run, multi-scale modeling approaches as presented by the authors should be superior to approaches based only on field calculations in explaining the biophysics of TMS. As such, I am very supportive of the study and the topic in general. In the following, I suggest changes to improve the clarity of the presentation of the results and to better align them with the known electrophysiological findings for magnetic stimulation of the motor cortex. In addition, I suggest extending the section in the discussion that deals with the study limitations. Given the novelty of the overall approach, it is relevant to educate the readers on the uncertainties involved in the modeling process, and to point toward putative improvements that could be taken up in future studies.
Comments on the methods and results section (sorted by order of occurrence, not relevance): Please specify the origin of the MR data and -if relevant -the underlying ethics approval (page 2, volume conductor model) The A field is not dependent on the head model, only coil geometry (page 2, field calculations) How much did the re-scaling of the neural models affect the results (page 3, multi-compartmental neuronal models)? This should be tested for a few selected cases. The question how well neuronal models for other species and brain areas can be transferred to human sensorimotor cortex is highly relevant. In that respect, it would be good to know how sensitive the results are towards changes of the local dimensions of the neurons.
Why were the axons of the L5 PNs modelled to curve towards the corpus callosum rather than the internal capsule (page 3)? How much could this have influenced the results?
Please give more details on the method used for placing the neurons and axons (page 3). Was it fully automated?
Adding an external current source to model the TMS field (page 3): How well is this justified for Adding an external current source to model the TMS field (page 3): How well is this justified for non-cable like structures such as the soma?
Angular frequency of 30 kHz (page 3): This is too high, 3 kHz would be more realistic. Maybe just a typo? Figure 1C: It would be helpful to indicate the region in which the neurons were placed in this plot. Figure 3: It is surprising to see that only slight changes in the field strength in dependence on coil orientation were observed. This is in contrast to prior results. The peak field strength might stay relatively constant, but the spatial pattern of the gyral crowns seeing high field strengths should clearly change. Maybe this is only a scaling issue in the plot, or related to the question on where the field was read out (on the CSF-GM surface or within GM?). Please clarify, and compare to prior studies (in particular those using the same FEM method -i.e. simnibs) Figure 3 and definition of coil orientations in the text: It would be easier to define the orientations according to the current direction induced in the brain. It is somewhat counterintuitive that that an arrow pointing anteriorly (i.e., in the base orientation) represents a current flowing in the opposite direction, as stated on page 2. Physiological experiments show that M1 has the lowest threshold for monophasic pulses when currents are flowing from posterior to anterior. The results in Figure 4 indicate that this should be the case for the base orientation. This, however, is in contrast to the statement in the methods that current direction was anterior-posterior in that case. Please clarify. A schematic plot of all coil (or better current) directions would be helpful to prevent confusion. In a related manner, please also add information on the meaning of blue and red colors for the lower half of Figure 3 (normal component). Do red colors indicate inflowing or outflowing currents? Figure 4: It would be good to be able to see the fields and excitabilities on the medial part of the hand knob. For example, for 270° (L5 PNs), the hand knob seems mostly non-excited. However, this impression might simply result from not being able to see the medial part.
Figures 4 and 5, and related results: I strongly suggest to differentiate between M1 and S1, and present separate results for both, with a focus on M1. From Figure 5, it seems that coil orientation 135° activates the largest number of neurons. However, Figure 4 indicates that those neurons are mostly in S1. For this reason, it would make more sense to present separate subfigures for M1 and S1 in Figure 5. Figure 6, and related results: It seems that axon bends of L5 PNs were never the site where stimulation occurred most easily. This is in contrast to prior hypotheses, and also (simplified) modeling studies. Can you indicate how much higher the excitation thresholds of the bends are? What is the reason that the bends were not excited? By which angle did the axons bend, and which bend radius was modelled? Figure 8: It would be good to replace the results using the coil orientation which stimulates the hand knob most strongly, rather than S1.
Comments on the Discussion section: This is one of only few studies that combine field calculations with morphologically realistic models of neurons. It shows very encouraging results. Given its pioneering character, however, I would appreciate if the uncertainties involved in the modeling process could be more systematically discussed.I would like to the uncertainties involved in the modeling process could be more systematically discussed.I would like to emphasize that this suggestion is not meant as a critique of quality of the study. The match between modeling and experimental results is promising, but should not be misunderstood as a proof that the model is "correct".
In addition to limitations in modeling the neurons, uncertainties regarding the tissue conductivities (resulting in uncertainties regarding the field distribution) should be mentioned. The neuron models were taken from cat visual cortex. Even in humans, the histology of visual and motor cortex differs substantially. While the visual cortex has a thick layer 4, the motor cortex is dominated by layer 3 and 5 (including the large "Betz" cells). The study results demonstrate a strong impact of neural features on the results, raising the question on how well a transfer between species and areas is possible. Here, the size of the neurons was rescaled to fit to the target cortical area. I am curious on how much is known about other differences, such as differences in axon and soma diameters, or channel densities, that might systematically impact the results, but where not taken into account? In a related manner, the study of Salvador (2011) et al. hints towards putative additional structures that might be affected, such as short-and intermediate connections in WM (their a1 and a2 association fibres), and terminals of incoming fibre projections. It would be worth mentioning them as well.
I am curious to understand which features of the L3 PNs resulted in their better excitability to TMS? One could intuitively assume that the larger and far-projecting L5 PNs might be the better targets. Discussing this might also help to better reveal the important neural features that dominate the cortical excitability to TMS. Thank you very much for your valuable comments, which have helped us greatly to improve the manuscript. We reply to your comments point by point below.

Please specify the origin of the MR data and -if relevant -the underlying ethics approval (page 2, volume conductor model)
We have used the SimNIBS pipeline to construct the head model and to calculate electric field distributions using the example dataset provided by SimNIBS. We rephrased the corresponding sentences as follows: "(At the beginning of Methods section) The volume conductor head model was used to simulate the stimulus-induced electric fields; it was based on the SimNIBS v1.1 software pipeline ." "(Volume conductor model in Method section) To calculate the precise electric field, a volume 14,33 "(Volume conductor model in Method section) To calculate the precise electric field, a volume conductor head model for TMS that reflected T1-weighted and T2-weighted magnetic resonance (MR) images was constructed using example dataset provided by SimNIBS v1.1 under the ethical approval " The A field is not dependent on the head model, only coil geometry (page 2, field calculations) We rephrased the sentence as follows: "The primary electric field was directly determined by the coil geometry and the secondary electric field caused by charge accumulations at tissue interfaces." How much did the re-scaling of the neural models affect the results (page 3, multi-compartmental neuronal models)? This should be tested for a few selected cases. The question how well neuronal models for other species and brain areas can be transferred to human sensorimotor cortex is highly relevant. In that respect, it would be good to know how sensitive the results are towards changes of the local dimensions of the neurons. We re-scale the dendritic trees while keeping axon and soma identical and observe the spatial distributions of thresholds (Figure 8(b)). We found only slight changes in threshold maps and the percentage of excited PNs by reducing the scale of dendritic trees, while the impact of rotation of dendritic trees was relatively bigger. However, in the presence of morphological changes of dendritic trees, the walls of the precentral gyrus were consistently targeted and the spatial extent of activated regions did not change much. We rephrased the corresponding sentences as follows: "(Results) In addition, the impact of scaling of PNs was investigated by reducing the dimension of the dendrite trees by 10%, 20% or 30%, as shown in Figure 8(b). For this, we simply 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." Why were the axons of the L5 PNs modelled to curve towards the corpus callosum rather than the internal capsule (page 3)? How much could this have influenced the results? The corpus callosum was a typo. We tilted the triangular elements comprising the gray matter defined in the ROI toward the internal capsule (Supplementary Figure S1). We rephrased these points as follows: "(Multi-compartmental neuronal models in Method section) The axons of L5 PNs were defined to curve beyond the boundary between GM and WM in the direction of the internal capsule (Supplementary Figure S1)."

Please give more details on the method used for placing the neurons and axons (page 3). Was it fully automated?
It is not fully automated. When the different neuronal models or different head models are applied, we need to change the parameters to adapt it. We rephrased corresponding paragraphs to include more details as follows: "L5/L3 PNs were combined virtually with the head model and modified to accommodate the irregular geometry of the cortex , as shown in Figure 2. To reduce superfluous   14   46 28, 29, 44-47 irregular geometry of the cortex , as shown in Figure 2. To reduce superfluous computations, we preselected a region of interest (ROI) of 50 × 50 × 50 mm 3 around the hand knob and then placed L5/L3 PNs in each triangular element comprising the gray matter surface. The multi-compartmental models of PNs consisted of a series of compartments connected by resistors. Each compartment was further discretized into segments of equal length to allow for accurate numerical simulation. The center points of each segment were extracted and used to calculate the necessary changes to neuron geometry, as described below. The dendritic trees were lengthened or shortened by re-scaling the lengths of the compartments according to the local dimensions of the cortex such that dendrites reached layer 1 and the orientation was perpendicular to the cortical surface . Since the morphology of the dendritic trees was not symmetric and it might influence the neuronal activation, L5/L3 PNs had randomly rotated dendritic trees at different locations. The axons of L5 PNs were defined to curve beyond the boundary between GM and WM in the direction of the internal capsule (Supplementary Figure S1). Further adjustments of L5 PNs geometry were as follows (illustrated in Supplementary Figure S2): each dendritic tree was oriented such that its principal axis would align with the normal vector of its associated triangular surface element. The bending part of the axon was calculated according to the normal vector of the surface element. The arc length of the axon bend was set to 0.6 mm when the z-component of the normal vector was positive and otherwise the arc length was 0.3 mm (compare Supplementary Figure S2). Note that when we varied the angle and arc length of the axon bend, it usually did not alter the activation threshold. The axons of L3 PNs were defined to terminate in layer 5/6 within the GM. Altogether, a total of 10,888 L5 PNs and 10,888 L3 PNs was constructed. This process was implemented in MATLAB (MathWorks, Natick, MA, USA)."     Figure 3 and definition of coil orientations in the text: It would be easier to define the orientations according to the current direction induced in the brain. It is somewhat counterintuitive that that an arrow pointing anteriorly (i.e., in the base orientation) represents a current flowing in the opposite direction, as stated on page 2. Physiological experiments show that M1 has the lowest threshold for monophasic pulses when currents are flowing from posterior to anterior. The results in Figure 4 indicate that this should be the case for the base orientation. This, however, is in contrast to the statement in the methods that current direction was anterior-posterior in that case. Please clarify. A schematic plot of all coil (or better current) directions would be helpful to prevent confusion. In a related manner, please also add information on the meaning of blue and red colors for the lower half of Figure 3 (normal component). Do red colors indicate inflowing or outflowing currents? Thanks for the suggestion. In the bottom row in Figure 3, red color indicated current directed inwards and blue color means current directed outwards. Thus, in the base orientation, the currents were flowing from posterior to anterior. We corrected the corresponding sentence as follows:

Adding an external current source to model the TMS field (page 3): How well is this justified for non-cable like structures such as the soma?
"(Field calculations) The base coil orientation was defined relative to the direction of the central sulcus that the electric field induced was in the posterior to anterior direction (the yellow arrow in Figure 1(c))." "(Legend in Figure 3) For the orthogonal component of electric fields (bottom row), red color indicates current flowing in the direction from superficial to lower cortical layers and blue color represents currents flowing in the opposite direction." Figure 4: It would be good to be able to see the fields and excitabilities on the medial part of the hand knob. For example, for 270° (L5 PNs), the hand knob seems mostly non-excited. However, this impression might simply result from not being able to see the medial part. We made Supplementary Figure S3 (an animated gif) to show the medial part of the hand knob. Figures 4 and 5, and related results: I strongly suggest to differentiate between M1 and S1, and present separate results for both, with a focus on M1. From Figure 5, it seems that coil orientation 135° activates the largest number of neurons. However, Figure 4 indicates that those neurons are mostly in S1. For this reason, it would make more sense to present separate subfigures for M1 and S1 in Figure 5. We appreciate the reviewer's suggestion. We separated results for precentral and postcentral gyrus, as shown in Figure 5 (for Figure 4, we already showed spatial distributions of thresholds separately for the precentral and postcentral gyrus). We rephrased corresponding sentences in Results and Discussion sections, as follows: "(Results section) 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." Figure 6, and related results: It seems that axon bends of L5 PNs were never the site where Figure 6, and related results: It seems that axon bends of L5 PNs were never the site where stimulation occurred most easily. This is in contrast to prior hypotheses, and also (simplified)  induced by magnetic stimulation also argued that action potential initiation is in the axon initial segment. Thus, these recent modeling studies reported the axon initial segment as action potential initiation site. The different observation compared to the results using simplified models might be induced because of different morphology and electrical properties of incorporated PNs. As we mentioned in the Discussion section, further study of the impact of neuronal morphology and electrical properties will be helpful. We rephrased corresponding sentences as follows: "(Discussion section) 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 (2015) combined a realistic et al. 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 (2011) have also et al. 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 . 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." In addition, when we control the arc length from 0.6 mm over 0.3 mm to 0 mm we found that the threshold did not alter as long as the bend was smooth, while a bend with 0 mm arc length produced an increased threshold. In this work, 0.6 mm arc length was determined to reproduce the results from Wongwarnpigoon , 2012 that investigated the impact of epidural cortical et al. stimulation using same pyramidal neurons. We added text how the axon bending parts were made as follows: "(Method -Multi-compartment neuronal models) The bending part of the axon was calculated according to the normal vector of the surface element. The arc length of the axon bend was set to 0.6 mm when the z-component of the normal vector was positive and otherwise the arc length was 0.3 mm (compare Supplementary Figure S2). Note that when we varied the angle and arc length of the axon bend, it usually did not alter the activation threshold." Figure 8: It would be good to replace the results using the coil orientation which stimulates the hand knob most strongly, rather than S1. We appreciate the reviewer's suggestion. We changed Figure 8 to the case for the base coil orientation.
Comments on the Discussion section: This is one of only few studies that combine field calculations with morphologically realistic models of neurons. It shows very encouraging results. Given its pioneering character, however, I would appreciate if the uncertainties involved in the modeling process could be more systematically discussed.I would like to emphasize that this suggestion is not meant as a critique of quality of the study. The match between modeling and experimental results is promising, but should not be misunderstood as a proof that the model is "correct".
A critical limitation involved in the proposed modeling process is that straightly stretched axons of L5 PNs inside the WM are not realistic. Thus, as we mentioned in discussion section, further developments in tractography using DTI may improve detailed neuronal morphology and thus we expect that more realistic PNs can be constructed in the future. The other uncertainties of the neuronal models were related to electrical and morphological variation among and within PNs. These issues were addressed below according to the next comments.
In addition to limitations in modeling the neurons, uncertainties regarding the tissue conductivities (resulting in uncertainties regarding the field distribution) should be mentioned. The neuron models were taken from cat visual cortex. Even in humans, the histology of visual and motor cortex differs substantially. While the visual cortex has a thick layer 4, the motor cortex is dominated by layer 3 and 5 (including the large "Betz" cells). The study results demonstrate a strong impact of neural features on the results, raising the question on how well a transfer between species and areas is possible. Here, the size of the neurons was rescaled to fit to the target cortical area. I am curious on how much is known about other differences, such as differences in axon and soma diameters, or channel densities, that might systematically impact the results, but where not taken into account? In a related manner, the study of Salvador (2011) hints towards putative additional et al. structures that might be affected, such as short-and intermediate connections in WM (their a1 and a2 association fibres), and terminals of incoming fibre projections. It would be worth mentioning them as well.
We agree that neural features have a strong impact on the results. In this work, we only varied the morphology of PNs (L3 and L5 PNs). However, further variation regarding to morphological and electrical properties should be considered. Thus, we rephrased the corresponding sentences as follows: "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 (2011) revealed that et al. 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 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 that incorporated the same models of PNs." I am curious to understand which features of the L3 PNs resulted in their better excitability to TMS? One could intuitively assume that the larger and far-projecting L5 PNs might be the better targets. Discussing this might also help to better reveal the important neural features that dominate the cortical excitability to TMS.
Mainen and Sejnowski, 1996, [41], tested the influence of dendritic structure by constructing neurons that share a common distribution of ion channels and differ only in their dendritic geometry (Fig 1). In that study, we used Fig.1  ). It is plausible therefore that the morphology of L3 PNs might contribute to their better excitability to TMS. We rephrased the corresponding sentences as follows: "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 , and the higher stimulation intensities required to produce D-waves." No competing interests were disclosed. Competing Interests: