A Model of TMS-induced I-waves in Motor Cortex

Brain Stimulation 7 (2014) 401e414

Contents lists available at ScienceDirect

Brain Stimulation journal homepage: www.brainstimjrnl.com

A Model of TMS-induced I-waves in Motor Cortex
ta
lin V. Rusu a, c, e, Max Murakami a, d, Ulf Ziemann b, Jochen Triesch a, d, * Ca a

Frankfurt Institute for Advanced Studies, Frankfurt am Main, Germany Department of Neurology & Stroke, and Hertie-Institute for Clinical Brain Research, Eberhard Karls University, Tübingen, Germany c Center for Cognitive and Neural Studies (Coneural), Romanian Institute of Science and Technology, Cluj-Napoca, Romania d Department of Physics, Goethe University, Frankfurt am Main, Germany e Department of Computer Science, Babes¸-Bolyai University, Cluj-Napoca, Romania b

a r t i c l e i n f o

a b s t r a c t

Article history: Received 8 February 2013 Received in revised form 17 February 2014 Accepted 17 February 2014 Available online 27 March 2014

Background: Transcranial magnetic stimulation (TMS) allows to manipulate neural activity non-invasively, and much research is trying to exploit this ability in clinical and basic research settings. In a standard TMS paradigm, single-pulse stimulation over motor cortex produces repetitive responses in descending motor pathways called I-waves. However, the details of how TMS induces neural activity patterns in cortical circuits to produce these responses remain poorly understood. According to a traditional view, I-waves are due to repetitive synaptic inputs to pyramidal neurons in layer 5 (L5) of motor cortex, but the potential origin of such repetitive inputs is unclear. Objective/hypothesis: Here we aim to test the plausibility of an alternative mechanism behind D- and I-wave generation through computational modeling. This mechanism relies on the broad distribution of conduction delays of synaptic inputs arriving at different parts of L5 cells’ dendritic trees and their spike generation mechanism. Methods: Our model consists of a detailed L5 pyramidal cell and a population of layer 2 and 3 (L2/3) neurons projecting onto it with synapses exhibiting short-term depression. I-waves are simulated as superpositions of spike trains from a large population of L5 cells. Results: Our model successfully reproduces all basic characteristics of I-waves observed in epidural responses during in vivo recordings of conscious humans. In addition, it shows how the complex morphology of L5 neurons might play an important role in the generation of I-waves. In the model, later I-waves are formed due to inputs to distal synapses, while earlier ones are driven by synapses closer to the soma. Finally, the model offers an explanation for the inhibition and facilitation effects in pairedpulse stimulation protocols. Conclusions: In contrast to previous models, which required either neural oscillators or chains of inhibitory interneurons acting upon L5 cells, our model is fully feed-forward without lateral connections or loops. It parsimoniously explains findings from a range of experiments and should be considered as a viable alternative explanation of the generating mechanism of I-waves. Ó 2014 Elsevier Inc. All rights reserved.

Keywords: Transcranial stimulation Motor cortex I-waves Computational model Neuron

Introduction Non-invasive brain stimulation techniques such as transcranial magnetic stimulation (TMS) have gained much attention in recent years after promising results in treating several neurological

Supported by the LOEWE-Program “Neuronal Coordination Research Focus Frankfurt” (NeFF). Financial disclosures: The authors declare having no conflicts of interest or financial interests related to this study. * Corresponding author. Frankfurt Institute for Advanced Studies, Ruth-MoufangStr. 1, 60438 Frankfurt am Main, Germany. Tel.: þ49 (0) 69 798 47531. E-mail address: triesch@fias.uni-frankfurt.de (J. Triesch). 1935-861X/$ e see front matter Ó 2014 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.brs.2014.02.009

disorders such as depression or stroke [1,2]. The ability to influence brain activity non-invasively is very appealing, but it has been difficult to establish how exactly TMS activates different types of neurons in cortical circuits. In a standard single-pulse paradigm a TMS coil is placed over the motor cortex. The fluctuating magnetic field induces an electric field, which affects the excitability of central motor pathways and causes a strong depolarization of large neuronal populations. As a result of single-pulse TMS but also single-pulse transcranial electric stimulation (TES), high-frequency (w600 Hz) descending volleys of activity are observed at electrodes in the epidural space [3e6]. Early studies [7e9] using different stimulation types and comparing spike latencies have shown that weak electrical stimuli generated brief single waves. They

402

C.V. Rusu et al. / Brain Stimulation 7 (2014) 401e414

suggested that such waves originate as a result of direct stimulation of pyramidal tract neurons and therefore these waves were termed D-waves. With strong cortical stimuli the D-wave was usually followed by three or four subsequent waves considered to have an indirect origin and thought to be the result of action potentials (APs) from presynaptic fibers impinging on the dendritic tree of pyramidal tract neurons. They were therefore termed I-waves. In humans, the descending cortical volleys of activity were initially recorded in anaesthetized subjects from the surface of the spinal cord during surgery. More recently, they have been recorded in conscious subjects with no central nervous system (CNS) abnormality who had spinal cord electrodes implanted for the treatment of chronic pain. Such opportunities are uncommon and thus the effects of TMS are usually studied by using needle electromyography (EMG) recordings or by observing limb movements. However, such recordings give little insight into the effects of magnetic stimulation on populations of cortical neurons. Several theories regarding the mechanism responsible for the generation of I-waves have been proposed, but none of them has gained widespread acceptance [7,10,11]. In one of the proposed models, TMS causes the activation of a large pool of cortical cells that fire according to their intrinsic membrane properties and act as a resonating circuit that directly activates layer 5 (L5) cells. Such a model however has been argued to require a second, indirect activation of layer 2 and 3 (L2/3) cells. In another proposed model, TMS activates chains of inhibitory and excitatory neurons providing waves of activation and inhibition to L5. It is unclear, however, what the anatomical basis for such hypothesized chains should be. Here we aim to test the plausibility of an alternative mechanism behind D- and I-wave generation through computational modeling that relies on the broad distribution of conduction delays of synaptic inputs arriving at different parts of L5 cells’ dendritic trees and their spike generation mechanism. We constructed a model consisting of a morphologically detailed L5 cell stimulated by a pool of L2/3 excitatory and inhibitory cells in a 4:1 ratio [12] (Fig. 1). These cells project uniformly onto the basal and apical dendrites of the L5 cell. Ion channel kinetics were modeled using a modified HodgkineHuxley formalism [13], while neurotransmission is mediated by excitatory voltage-independent and voltage-dependent channels together with inhibitory channels (see Methods section). The model offers a

viable explanation of the generation of D- and I-waves and accounts for their frequency and timing. According to the model, the generation of I-waves is the product of intrinsic and extrinsic factors. Synchronous volleys of excitatory and inhibitory post-synaptic potentials (EPSPs and IPSPs) from L2/3 cells interact on the complex dendritic tree of the L5 cell. In response, the L5 cell’s spiking mechanism generates brief trains of action potentials at typical I-wave frequencies. We explore how the anatomical details of L5 cells play an important role in I-wave generation. Specifically, we show how by removing distal dendritic branches later I-waves are abolished. Thus, our model is consistent with the claim that individual I-waves have different generating mechanisms with early and late I-waves being supported by distinct pathways [9,14]. As shown in the following, our model reproduces findings from a range of experiments including pharmacological or behavioral modulation of I-waves and facilitation and inhibition effects observed in paired-pulse stimulation protocols. Overall, we find that the proposed model is a viable and parsimonious alternative to previous attempts to explain the mechanisms behind D- and I-wave generation. Methods The model All simulations were carried out in Python using the NEURON simulation environment [15,16]. The computational model presented here was simulated with a time step of 0.025 ms. Our model incorporates an L5 pyramidal cell and a total of 300 excitatory and inhibitory L2/3 cortical neurons (ratio 4:1) projecting onto it (Fig. 1). The morphology of the L5 cell was reconstructed from the experimental studies presented in Ref. [17] and modeled as a multicompartmental unit (188e201 dendritic branches, 11 segments for each compartment, average segment length of 10 mm; for a detailed morphological overview see original publication). The L2/3 cells The modeling of the L2/3 excitatory and inhibitory cortical neurons as regular and fast spiking single-compartment units was taken from Ref. [18]. The L2/3 neurons were modeled as cylinders with their radius equal to their height, measuring 96 mm and 67 mm for the excitatory and the inhibitory L2/3 neurons, respectively. The dynamics of the L2/3 single compartments were governed by:

Cm

dV ¼ g leak ðU  Eleak Þ  INa  IKd  IM ; dt

(1)

mF where U is the membrane potential, Cm ¼ 1 cm 2 the specific mS 4 membrane capacitance, g leak ¼ 10 the leak membrane 2 cm conductance, Eleak ¼ 70 mV its reversal potential and Ix the membrane current of ion type x. Specifically, INa and IKd are the sodium and delayed rectifier potassium currents responsible for action potential generation, and IM is a slow voltage-dependent potassium current responsible for spike-frequency adaptation. Details for each current are given in the Appendix.

The L5 cell For the dynamics of the L5 cell a modified version of the HodgkineHuxley formalism was used [19]: Figure 1. The used model including a reconstructed dendritic tree of an L5 pyramidal cell. A total of 300 excitatory and inhibitory L2/3 cells (ratio 4:1) project synapses onto the L5 cell.

Cm

dU ¼ g leak ðU  Eleak Þ  INa  IK ; dt

(2)

C.V. Rusu et al. / Brain Stimulation 7 (2014) 401e414

403

mF where U is the membrane potential, Cm ¼ 1 cm 2 the specific mS the leak membrane membrane capacitance, g leak ¼ 104 cm 2 conductance, Eleak ¼ 70 mV its reversal potential and Ix the membrane current of ion type x. Specifically, INa and IK are the sodium and potassium currents. Additionally, the L5 axon featured a delayed-rectifier potassium current IKDR which was adapted from Ref. [20]. Details for each current are given in the Appendix.

the soma of the L5 cell. Figure 2A shows a simulated voltage trace recorded at the L5 cell axon exhibiting fast APs in response to a somatic current injection of 3 nA for 50 ms. Figure 2B plots the mean firing rate over 50 ms against the intensity of the injected current, which was varied in steps of 0.1 nA. Indeed, with sufficient stimulation, the neuron generates firing rates of around 600 Hz, matching typical I-wave frequencies (6 spikes in 10 ms).

Frequency response

The synapses

During voluntary movement the firing rate of pyramidal tract neurons has been found to rarely exceed 100 Hz [21e23]. While this is much lower than the typical 600 Hz corresponding to I-waves, recordings from single corticospinal neurons and thalamocortical projection neurons have shown that these neurons can indeed fire up to I-wave frequencies [7,8]. To ensure that the simulated L5 cell is capable of firing spikes at I-wave frequencies, we slightly adjusted potassium channel conductances (from the original parameters [17]) to narrow individual action potentials (APs) [24]. To test this we simulated a direct injection of rectangular current pulses into

Each of the L2/3 neurons projected synapses onto the L5 cell with neurotransmission modeled by excitatory voltageindependent a-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA) channels and voltage-dependent N-methyl-D-aspartate (NMDA) channels together with inhibitory g-aminobutyric acid-A (GABAA) and g-aminobutyric acid-B (GABAB) channels. Channel conductance changes according to:


    t  sc t  sc  exp  ; gðtÞ ¼ gpeak exp 

s1

(3)

s2

where s1 is the rise time and s2 the decay time of the synapse. Here sc is a 1 ms conduction delay between an L2/3 cell generating a spike and this spike reaching its target location on the L5 cell’s dendritic tree. This conduction delay replicates the observed delays between D- and I-waves (typically between 0.75 and 1.34 ms [25]). Time constants, reversal potentials and peak conductances for each synaptic channel (see Refs. [26e31] for typical values) are given in Table 1. Synapses were distributed uniformly onto the dendritic tree. Figure 3 depicts the targets of the L2/3 excitatory projections. This placement, while simplistic, matches experimentally recorded data in the L5 of the cortex [32]. Synaptic strengths, modeled as a gain factor, were drawn from log-normal distributions. For inhibitory synapses, log WLI wN ð0:2; 0:5Þ. For excitatory synapses, log WLE wNð0:1; 0:5Þ. This essentially caused some synapses to be very strong while most were weak. As a consequence, depending on which L2/3 fibers were activated by the TMS pulse, the number of output spikes in the L5 neuron differed. Synaptic depression has been shown to be a prominent feature of transmission at neocortical synapses and is believed to be caused by depletion of the readily releasable pool of neurotransmitter vesicles [33e35]. In the present model, synapses from L2/3 to L5 exhibited a form of short-term synaptic depression. Specifically, synaptic response amplitudes were modulated using a depression factor D [36]. When the synapse was in the undepressed state, D ¼ 1. After each presynaptic spike, D was multiplied by a factor d ¼ 0.5 representing the amount of depression per presynaptic action potential:

D/D d:

(4)

Between stimuli, D recovered exponentially back toward one with time constant s (200 ms):

1D : D_ ¼

(5)

sD

Table 1 Time constants and reversal potentials for each synaptic channel as derived from the literature.

Figure 2. High-frequency firing of the L5 cell in response to direct current injection. (A) Repetitive axonal response evoked by the injection of a rectangular current of 3 nA (horizontal bar). (B) Mean firing rate as an average over time (50 ms after stimulation) vs. current amplitude.

Receptor

s1 (ms)

s2 (ms)

Erev (mV)

gpeak (mU1)

AMPA NMDA GABAA GABAB

0.2 2 0.3 45.2

1.7 26 2.5 175.16

0 0 70 90

0.1 0.03 0.5 0.5

404

C.V. Rusu et al. / Brain Stimulation 7 (2014) 401e414

majority of the induced current was present for less than a quarter of the total time after which it decayed. Figure 4 plots the simulated time course of the TMS-induced magnetic and electric field. As a result of TMS pulses directly activating the L5 cell, i was directly injected into the axon. To model different levels of excitability of pyramidal tract cells in response to TMS pulses, current intensities (ipeak) were drawn from a normal distribution: DN wNð2:7; 0:5Þ. Depending on this current, the L5 cell emitted or not a direct response as a result of the TMS pulse. The magnetic field induced by a TMS pulse directly activates both the pyramidal tract and L2/3 cortical fibers; however, it rapidly decays with distance from the surface of the skull. As a consequence, neurons located in deeper layers are likely to be less influenced by the magnetic field. In addition, the orientation of fiber termination relative to the magnetic field also determines how neurons are stimulated. To model this, we selected a percentage of L2/3 cells that were activated by the TMS pulse in addition to the L5 cell (Fig. 1). These cells generated a spike while the remaining others received a smaller depolarizing current (typically, half of the current threshold). During simulations additional shapes for the induced electric current (such as monophasic pulses, rectangular waves) were also considered. However, we found that the precise shape does not strongly impact the results, with the total current being the decisive factor. Simulating epidural recordings During experiments on conscious humans, the responses to TMS pulses have been recorded by placing electrodes in the epidural space [4,42,43]. Such recordings measure a compound response from axons of a large neuron population. To simulate comparable recordings, multiple instantiations of the model were simulated with L5 axonal voltage traces superimposed and convolved with a Gaussian filter (s ¼ 10) to account for different conduction delays from L5 cells to where the epidural recording is made.

Figure 3. Density of excitatory L2/3 to L5 projections.

Simulating TMS pulses Commercial TMS devices typically use coils with a diameter of 50e100 mm and induce a monophasic magnetic field of 1e2 T with a rise time of about 50e200 ms and a total duration of roughly 0.4 ms [37]. The shape of the coil has been shown to have an impact on the efficiency of magnetic stimulation with different shapes having been proposed to improve behavioral responses. According to Faraday’s law the fluctuating magnetic field induces an electric field across the scalp which generates currents across the cellular membrane of a neuron and causes a change in polarization. In the case of a monophasic pulse the first phase of the induced electric field roughly resembles a quarter cycle of a cosine wave followed by a relatively low electric field in the opposite direction [38]. Recently, TMS devices capable of generating controllable pulse widths and near-rectangular pulse shapes have been proposed and argued to be more efficient [39,40]. In order to match the result of [38] qualitatively, we modeled the induced current as an underdamped current of an RLC circuit [41]:

iðtÞ ¼ ipeak sin ðu tÞ exp

  t  ;

s

(6)

where s ¼ 0.08 ms is the decay time of the pulse and ipeak is the maximal amplitude in nA. This way, similar to TMS devices, the

Figure 4. Simulated electric field shape induced by a TMS pulse (arbitrary units).

C.V. Rusu et al. / Brain Stimulation 7 (2014) 401e414

L5 morphology generation To construct L5 morphologies, we applied the following algorithm: (i) create a list of all dendritic branches, from both the apical and basal dendrites; (ii) compute the individual path lengths from each individual branch to the soma; (iii) sort the list with respect to the path length; (iv) according to a predefined percentage P, keep the closest branches and remove the rest; (v) uniformly place all synapses (their number remains thus constant). The parameter P determines the size of the resulting dendritic tree, with P ¼ 100 referring to the full tree and P ¼ 0 to a point neuron, which includes just the soma and axon. TMS is thought to elicit responses from L5 cells through both direct and indirect activation. Direct activation is modeled as a brief direct current injection of variable strength into the L5 cell which sometimes causes the appearance of a single spike contributing to the D-wave. Indirect activation is caused by stimulation of presynaptic L2/3 neurons, which leads to a series of subsequent spikes in the L5 neuron contributing to I-waves. In vivo recordings of epidural responses typically show a D-wave followed by three or four larger I-waves at peak intervals of about 1.4 ms with the D-wave clearly separated from subsequent I-waves. Recordings from single neurons have shown that the delays between the spikes corresponding to D- and I-waves form a continuous distribution between 0.75 and 1.34 ms [25]. We model this delay by introducing a 1 ms conduction delay for connections from L2/3 to L5 (see Methods section). To simulate D- and I-waves we pool axonal voltage traces of different model instantiations with different random synaptic connections drawn from log-normal distributions, different synapse placement and different direct activation strengths drawn from a normal distribution (see Methods section). Figure 5A shows a simulated epidural response to a TMS pulse activating 100% of the excitatory and inhibitory L2/3 inputs obtained from 100 model instantiations. A smaller D-wave is followed by three strong I-waves and a small fourth one. An

Figure 5. Epidural response to a TMS pulse. (A) Simulated epidural recording consisting of a small D-wave followed by three larger I-waves and a small fourth. (B) Epidural recording of the response generated by TMS of the human motor cortex. The vertical line indicates the timing of the TMS pulse. Adapted from Ref. [42].

405

epidurally measured response adapted from Ref. [42] is shown for comparison (Fig. 5B). Epidural and EMG recordings have shown that as the strength of the stimulation is increased, the strength of individual I-waves and their overall number increase. We tested if I-waves generated by our model are also recruited in a dose dependent fashion and how their number depends on the balance between inhibition and excitation of the L5 cells due to L2/3 cells. Depending on the strength and orientation of the TMS pulse, different percentages of inhibitory and excitatory L2/3 fibers get activated. Inhibitory fastspiking (FS) cells have smaller somata compared to pyramidal regular-spiking (RS) cells and have been suggested to be more sensitive to TMS [44,45]. We model this by having a weak TMS pulse preferentially activating inhibitory cells (see Methods section). Figure 6 shows simulated epidural recordings for increasing fractions of active excitatory and inhibitory inputs. In line with experimental data [46,47], as the relative strength of excitation increases, more I-waves are generated. Cellular morphology dependency of I-waves Among cortical cells, pyramidal neurons are special because of their pyramidal-shaped somata and large dendritic trees that span several cortical layers. Because of the larger distance they need to travel, EPSPs originating from distal dendritic branches take longer to reach the soma. Hence, we hypothesize that the shape of L5 cells influences I-wave generation. Specifically, we studied how the size of the modeled L5 dendritic trees affects the amplitude and number of individual I-waves across epidural responses. To investigate this, we created L5 morphologies where different fractions P of dendritic branches were present during stimulation (see Methods section). The percentage P determined the size of the resulting dendritic tree, with P ¼ 100% referring to the full tree and P ¼ 0% to a neuron, which included just the soma and axon (Fig. 7). Importantly, we distributed the same number of synapses onto the different L5 morphologies to make results comparable. For each

Figure 6. Integrated epidural responses measured for increasing E:I ratio. Inhibition is fixed at 100% while excitation is increased in steps of 10%. I-waves are recruited in a dose-dependent manner: as more excitation is present, the number and size of I-waves increases.

406

C.V. Rusu et al. / Brain Stimulation 7 (2014) 401e414

Figure 7. L5 dendritic tree with different percentages of dendritic branches present. Starting from the full dendritic tree on the right, branches are sorted and removed based on their distance to the soma.

obtained dendritic tree, we recorded the number of spikes generated by the L5 cell for different placements of L2/3 synapses. To assure that results are not dependent on a specific dendritic tree, we additionally ran simulations using the cell from Ref. [48] with the same membrane and synaptic parameters (see Methods

section). Figure 8A shows the average number of recorded spikes (mean and standard error of 100 repetitions) for different values of P. As the percentage of removed dendrites increased, the number of L5 spikes decreased. Figure 8B shows the corresponding epidural recordings for different values of P. The reason for this effect is that it may take up to 10 ms for EPSP peaks originating in distal branches to reach the soma (Fig. 8C). Thus, by removing distal dendritic branches of the L5 cell, later I-waves are abolished. To further illustrate this point, we considered a “fan” neuron consisting of a soma (of a fixed diameter) connected to a number of identical cylindrical dendrites (Fig. 9A). A number of 300 synapses were placed uniformly on the dendrites. Membrane and synaptic parameters were not modified (see Methods section). By increasing their length for fixed numbers of dendrites (10) and recording the number of generated spikes (mean and standard error of 100 repetitions) we observed an increase in the number of generated spikes (Fig. 9B). EPSPs originating from distal portions caused a longer but smaller depolarization of the soma. As dendrites grew

Figure 8. Simulated epidural responses for different dendritic tree sizes. (A) Number of L5 spikes (average and standard error) out of 100 model instantiations with different percentages of dendrites removed. To assure results are not dependent on a specific morphology we also simulated the morphology presented in Ref. [48] (neuron 2) with the same synaptic and membrane channels (see Methods section). (B) Corresponding epidural responses from neuron 1 for different percentages of dendrites present (0, 60 and 100). (C) Probability density of EPSP peak timings recorded at the L5 soma for activations originating in different dendrites. The further an EPSP originates from the soma, the longer it will take to arrive.

C.V. Rusu et al. / Brain Stimulation 7 (2014) 401e414

407

Figure 9. Spike count versus length of dendrites for the fan neuron. (A) A sketch of the fan neuron morphology. (B) Number of L5 spikes (average and standard error) out of 100 model instantiations with different length of dendrites. As the length increases, more spikes are generated.

longer, such EPSPs delivered a long lasting excitation at the soma after a delay and enabled the generation of late I-waves. Additionally, by increasing the number of dendrites, we also observed an increase of the number of generated spikes (Fig. 10A). Figure 10B shows corresponding epidural responses for different numbers of attached dendrites. Inhibition close to the soma suppresses EPSP conduction more effectively than distal inhibition [49]. Thus EPSP transmission of individual dendrites may be suppressed if inhibitory L2/3 neurons project onto the somatic end of those dendrites. As the total number of dendrites grew, synapses were distributed more sparsely onto the dendritic tree. Due to the larger number of excitatory compared to inhibitory L2/3 neurons, a growing number of dendrites were not (or hardly) affected by inhibition, as the fraction of suppressed dendrites decreased. These effects led to a higher depolarization at the soma for an increasing number of dendrites. We have shown that both the number and length of dendrites affect the generation of I-waves; however, it is not evident which of these effects is more prominent. To investigate this, we distributed a fixed amount of dendritic length (1800 mm) among identical

dendrites (Fig. 11A). We note: (i) when the number of dendrites was low but their length large, we observed no I-waves; (ii) when the number of dendrites was larger but their overall length small, we observed a decrease in late I-waves; (iii) in an intermediate case when dendrites were sufficiently long and insufficient numbers, we observed an increase in late I-waves. Figure 11B shows corresponding epidural responses for different numbers of dendrites. Overall, our results suggest that the elongated dendritic trees of L5 pyramidal cells play an important role in the generation of I-waves: proximal excitation favors the growth of early I-waves, while distal excitation drives late I-wave production. I-wave modulation We reasoned that if our model of I-wave generation were essentially correct, then it should also be able to capture the effects of pharmacological and behavioral interventions on I-waves and results obtained with different stimulation protocols such as paired-pulse stimulations. Thus we set out to test our model’s ability to account for these forms of I-wave modulation.

Figure 10. Spike count versus number of dendrites for the fan neuron. (A) Number of L5 spikes (average and standard error) out of 100 model instantiations with different number of dendrites added. (B) Corresponding epidural recordings for different numbers of dendrites.

408

C.V. Rusu et al. / Brain Stimulation 7 (2014) 401e414

Figure 11. Spike count versus number of identical dendrites for the fan neuron with a fixed total dendritic length. (A) Number of L5 spikes (average and standard error) out of 100 model instantiations with different numbers of dendrites added. (B) Corresponding epidural recordings for different numbers of dendrites.

Pharmacological interventions The administration of CNS active drugs with known modes of actions has been shown to enable the use of TMS as a measure of cortical excitability (for a review see Refs. [50,51]). We set out to test if simulated changes to cortical excitability induced by pharmacological interventions produce alterations to I-waves similar to those observed during in vivo experiments. GABAergic neurons are the main source of inhibition in the cortex, with synaptic interaction mediated by either fast g-aminobutyric acid-A (GABAA) or slow g-aminobutyric acid-B (GABAB) receptors [52]. Benzodiazepines (such as lorazepam, diazepam or midazolam) bind to GABAA receptors causing hyperpolarization in the cellular membrane and affecting cortical excitability. Several studies have reported a decrease in amplitude of responses to TMS pulses after administration of GABA agonists [53,54]. These changes have been attributed to the enhancement of GABA-ergic activity, which caused a decrease in excitability at the level of the pyramidal tract. While such studies mainly relied on measuring evoked EMG responses, others have shown how the presence of GABA-ergic drugs affects the amplitude of epidurally recorded I-waves [5,55]. They reported that later I-waves lost amplitude while the D-wave and first I-wave remained unaffected. Regarding increased cortical excitability, while no direct pharmacological evidence of increasing amplitudes of the descending volleys exists, it was shown that when subjects increased cortical excitability by tightly clenching their fist, a TMS pulse which normally evoked three I-waves generated an additional one [5]. We tested whether our model shows similar changes to I-waves in response to increased inhibition or excitation. A decrease in excitability was simulated as an increase in GABA conductance while an increase in excitability was simulated as an increase in aamino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA) conductance. When the peak conductance of GABAA and GABAB channels was increased, the number of L5 spikes decreased (Fig. 12A). Interestingly, the GABAA channel has a more pronounced effect on I-waves because of its shorter time scale. The influence of adjusting AMPA conductance was also tested. By increasing the peak conductance of AMPA channels, the number of L5 spikes increased (Fig. 12B). Figure 12C shows the corresponding simulated epidural responses for different AMPA, GABAA and GABAB channel conductances. Each epidural recording was obtained from 100 model simulations (see Methods section). The size of the D-wave

was not altered by any of the interventions, since it is generated by direct pyramidal tract axon excitation and does not rely on synaptic interactions. Paired-pulse stimulation In a final set of simulations, we tested whether our model can also offer an explanation for findings from so-called paired-pulse stimulation protocols. In a well-established protocol, a sub-threshold TMS pulse over motor cortex that in itself is insufficient to produce an EMG response is delivered preceding a supra-threshold pulse. It has been observed that at short interstimulus intervals (ISIs) between the two pulses, the second pulse produced a weaker response compared to a control situation without the prior sub-threshold pulse. For longer ISIs, however, facilitation was observed [56]. To model this set of findings, we incorporated short-term synaptic depression into the synapses from L2/3 to L5 (see Methods section). It has been suggested that a sub-threshold TMS pulse mostly activates inhibitory neurons [44,45] because of their lower threshold for magnetic stimulation. To capture this we set a sub-threshold TMS pulse to activate 100% of inhibitory and 10% of excitatory L2/ 3 cells. The supra-threshold pulse was chosen to activate 100% of inhibitory and 100% of excitatory cells and was capable of generating three fully grown I-waves in a baseline condition without a prior sub-threshold pulse. For each individual I-wave in the integrated epidural response, the difference between its size as a result of paired- and single-pulse stimulation was measured for different ISIs. Figure 13A illustrates this change as a percentage of the baseline response resulting from single-pulse supra-threshold stimulation. For short ISIs of less than 5 ms, the amplitude of the second and third I-wave was strongly reduced. This is because the initial subthreshold-pulse activated mostly inhibitory L2/3 cells whose IPSPs still dampen the I-waves produced by the second (supra-threshold) pulse. In contrast, for ISIs greater than 5 ms, Iwaves increased in amplitude. This is because synapses from inhibitory L2/3 cells activated during the initial sub-threshold pulse were still in a depressed state and generate less inhibition in response to the second pulse. The cross-over from inhibition to facilitation depended on the time course of IPSPs and the time scale of synaptic depression. In a second stimulation protocol, a sub-threshold TMS pulse over motor cortex that in itself is insufficient to produce an EMG

C.V. Rusu et al. / Brain Stimulation 7 (2014) 401e414

409

Figure 12. Pharmacological interventions influence TMS responses. (A) For fixed AMPA peak conductance (0.2 mU1), increasing GABAA and GABAB conductances causes the number of L5 spikes to decrease. (B) For fixed GABAA and GABAB conductances (0.9 and 0.1 mU1 respectively) an increase in AMPA conductance causes the number of L5 spikes to increase. (C) Example epidural responses resulting from averaging 100 spike trains for different GABAA and GABAB conductances. Reduced excitability causes a decrease in size of late I-waves. Increased excitability causes an increase in size of late I-waves.

response was delivered following a supra-threshold pulse. In this case the stimulation led to facilitation in a wave-like pattern with peaks at ISIs of 1.2, 2.3e3.0 and 4.1e4.5 ms [57]. To model this set of findings, we chose a sub-threshold TMS pulse which activated 100% of inhibitory and 10% of excitatory L2/3 cells, and a supra-threshold pulse which activated 100% of inhibitory and 60% of excitatory cells and was capable of generating three I-waves. Compared to control, the change in integrated epidural responses obtained from pairings in the range of 1.2 and 4.5 ms followed a wave-like pattern similar to experimental observations (Fig. 13B) [57]. Although a single subthreshold pulse by itself was not strong enough to cause the

generation of L5 spikes as it activated few presynaptic elements, in conjunction with a preceding supra-threshold one at short ISIs it caused the generation of more L2/3 spikes (Fig. 13C). The amplitude of the facilitation is comparable to experimental data and other modeling work [11]. Discussion Magnetic stimulation has shown great potential for treating disorders such as stroke or depression [1,2]. Furthermore, pairedpulse stimulation paradigms, consisting of either paired TMS

410

C.V. Rusu et al. / Brain Stimulation 7 (2014) 401e414

Figure 13. Paired-pulse responses. (A) A sub-threshold pulse followed by a suprathreshold one was delivered at varying ISI values. For each ISI, the size of individual I-waves in the epidural response is compared to the size in the response to the suprathreshold stimulus alone. (B) A sub-threshold pulse preceded by a supra-threshold one was delivered at varying ISI values. For each ISI, the integrated epidural response is compared to the one in the response to the supra-threshold stimulus alone. (C) Voltage trace (black line) showing an L2/3 excitatory cell becoming depolarized but not spiking following a supra-threshold (I) TMS pulse (see Methods section). Subsequent stimulation by a sub-threshold (II) pulse causes the neuron to spike (green line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

pulses or pairs of low-frequency median nerve stimulation and TMS pulses, have been used to further explore the excitability of underlying cortical circuits and show the emergence of activitydependent forms of cortical plasticity [57e60]. However, for a successful and safe treatment the precise effects of TMS need to be well understood. This is particularly important considering the risks of repeated TMS [61] and studies which report adverse effects such as an impairment of working memory by single- and doublepulse TMS applied to the prefrontal cortex [62,63]. Simulations of detailed cortical models can help to better explain the effects of magnetic stimulation on cortical circuits. We investigated the effects of TMS on individual cells by constructing a model featuring a morphologically detailed L5 cell driven by a pool of L2/3 pyramidal cells and simulating TMS pulses as current injections. Our goal was to develop a parsimonious computational model explaining the TMS-induced generation of Dand I-waves measured in epidural recordings. In contrast to a traditional view, according to which I-waves are generated by repetitive synaptic inputs onto pyramidal tract neurons, we wanted to assess the viability of an alternative account relying on pyramidal tract neurons’ spike generation mechanism and synchronous synaptic inputs onto their extended dendritic trees. Our model demonstrates that this hypothesis can be implemented consistently with plausible choices of model parameters. Specifically, it captures different features of TMS: (i) realistic timings and frequency of Iwaves; (ii) dose-dependent recruitment of I-waves; (iii) pharmacological interventions on I-waves; (iv) results obtained during paired-pulse stimulation protocols. Previous modeling work either focused on constructing largescale models of the thalamocortical system [11], or on modeling the effects of TMS on simple pyramidal cells with arbitrary morphology but without explaining epidural responses [41,64,65]. Given this, there is still a debate on the exact mechanisms underlying the production of I-waves in the motor cortex, with several theoretical models proposed over the years [10,66]. A model in which TMS induces a large activation of L2/3 cells, which fire according to their intrinsic membrane properties and act as a resonating circuit directly activating the L5 neurons (Fig. 14A) [67], is especially relevant for our discussion. Such a model has been argued to require a second generation of spikes in L2/3 cells triggered indirectly by TMS. By using a detailed compartmental model with a rich dendritic tree we showed that such second spikes are not necessary. Our model is fully feed-forward featuring no lateral connections or loops. Depending on synaptic receptor kinetics and synapse positions, a spike generated in L2/3 affects the L5 cell with a certain delay. These dispersed dendritic inputs interact with the intrinsic membrane properties to cause the generation of periodic Iwaves. Another relevant model proposes that TMS activates a chain of inhibitory and excitatory neurons providing waves of activation and inhibition to L5 [14] (Fig. 14B). Again, our model shows that such chains are not required to explain the basic phenomena of Dand I-waves. Instead, we argue that I-waves are a product of both extrinsic and intrinsic factors. According to our model, a suprathreshold TMS pulse causes a synchronized volley of postsynaptic potentials from L2/3 cells to impinge onto different parts of the dendritic trees of L5 cells. The intrinsic membrane properties and spiking mechanism of the L5 cells are then responsible for generating trains of action potentials at the characteristic I-wave frequency. Our model shares some similarities with previous work presented by Esser et al. [11]. This model simulated a complex multi-layered thalamocortical region of the cortex, but it neglected the complex dendritic tree structure of L5 cells. Additionally, D-waves were not specifically modeled, and TMS pulses evenly

C.V. Rusu et al. / Brain Stimulation 7 (2014) 401e414

Figure 14. Hypothetical models of I-wave generation. (A) An L5 cell is stimulated by a pool of cortical neurons acting as a resonating circuit, providing rhythmic input to the L5 cell. (B) An L5 cell is stimulated by a chain of inhibitory and excitatory neurons.

activated both L2/3 and L5 cells, with projections from L2/3 to L5 being purely excitatory. While such a setup can capture how I-waves are formed (increased excitation coupled with a refractory mechanism), it does not reveal how distinct presynaptic populations targeting different parts of L5 dendritic trees may interact to generate different I-waves. We also studied the effects of paired-pulse stimulation protocols. A number of paired-pulse paradigms evoking responses that are either inhibited or facilitated depending on the strength, order and timing of the pulses exist. To model a first set of observations we integrated short-term synaptic depression of synaptic connections from L2/3 to L5. We showed that under these conditions a paired-pulse stimulation protocol, consisting of a sub-threshold TMS pulse followed by a supra-threshold activation, caused depression of epidural responses at short ISIs and facilitation at large ISIs in the model, which is in agreement with experiments [56,58,68]. The inhibitory effect for short ISIs is owed to the initial subthreshold-pulse, which is assumed to activate mostly inhibitory L2/3 cells (since they have been suggested to have a lower activation threshold [44,45]) and whose IPSPs will still dampen the I-waves produced by the second (supra-threshold) pulse. The facilitation for larger ISIs in the model is caused by the fact that synapses from inhibitory L2/3 cells activated during the initial subthreshold pulse will still be in a depressed state and generate less inhibition in response to the second pulse. In a second set of observations, a supra-threshold pulse followed by a sub-threshold one caused facilitation of epidural responses. In our model, the first pulse caused the generation of some L2/3 spikes, while the second one caused some of the remaining L2/3 cells to generate action

411

potentials, and thus more EPSPs impinged onto L5 cells. This led, compared to control, to facilitation of integrated epidural responses in the range of 1.2 and 4.5 ms in a wave-like pattern similar to experimental observations (Fig. 13B) [57]. Our model was also shown to reproduce the effects of pharmacological interventions on the size and number of I-waves. Such effects were modeled as changes in GABA and AMPA channel conductances. We have shown that, for an increase in GABA conductance, the amplitude of epidural responses decreased, while for an increase in AMPA conductance their amplitude increased. Most importantly, we have observed that later I-waves are more affected by such interventions than the D-wave and the first I-wave, consistent with experimental data [5]. Finally, we showed how the morphology of L5 cells affects simulated epidural responses. Using two distinct dendritic trees [17,48], from which we removed fixed percentages of dendritic branches (ordered by their distance to the soma), we showed how later I-waves are abolished when distal dendritic segments were cut off. These results suggest that later I-waves are formed due to distal synapses, while earlier ones are driven by synapses closer to the soma. The model also makes predictions about how the morphology of L5 cells influences epidural responses. Assuming an approximately uniform L2/3 synapse distribution on the L5 dendritic tree, we expect the following observations. First, ablating distal dendritic parts results in abolition of late I-waves. Early I-waves are barely affected. The same is expected for ablation of apical excitatory synapses or intensification of apical inhibition. Second, ablating basal dendrites results in diminishment of both amplitude and frequency of the early I-waves. The I-wave duration may increase. Third, intensification of distal excitation favors the growth of later I-waves and the extension of the I-wave train. Early I-waves are barely affected. The same is expected for ablation of apical inhibitory synapses. It has been observed that TES can produce similar results to TMS. It is thought that TES activates the corticospinal neurons primarily directly close to the axon hillock and, in addition at sufficiently high stimulus intensity, transsynaptically [9,14,47]. This is in contrast to TMS, where a D-wave due to direct excitation only occurs clearly above motor evoked potential (MEP) threshold intensity [4,43,47]. Our model can describe this situation quite easily, although it cannot explain it, since it does not explicitly describe how electric fields generated by TMS or TES produce neuronal firing. This has been the topic of previous research [69e71]. In our model, the activation of neurons due to electromagnetic fields is simply modeled as direct current injections. We can model TES stimulation by choosing current injections consistent with the effects of TES stimulation on individual neurons. Specifically, we can model weak or moderate intensity TES as current injections into L5 pyramidal cells that are sufficiently strong to make them fire, while modeling high intensity TES as additional current injections into L2/3 cells, which are sufficiently strong to make these fire. This way, the basic phenomenology of TES stimulation can easily be expressed by our model. In summary, by incorporating detailed L5 pyramidal cells we were able to develop a parsimonious model capable of reproducing D- and I-waves and offering an explanation of some of the biophysical mechanisms underlying magnetic stimulation. As inherent with all modeling work we made a number of simplifications. First, we have very much simplified cortical anatomy in our model. For example, we only modeled one type of fast spiking interneuron despite the large diversity with different functional roles observed in cortex [49,72,73]. We have also placed excitatory and inhibitory synapses uniformly on the L5 cell dendritic tree. A more careful placement based on anatomical data would be desirable.

412

C.V. Rusu et al. / Brain Stimulation 7 (2014) 401e414

Furthermore, we have completely neglected recurrent connectivity e either within L5 or between L2/3 and L5. Our model is sensitive to the chosen time constants. Our time constants are taken from actual experimental measures which however might not explicitly focus on pyramidal cells in motor cortex (see Methods section). Lastly, our model features only modified versions of Hodgkine Huxley currents and neglects other complex currents known to play an important role in cortical neuron dynamics [74]. We have tested the model with a more complex L5 cell featuring additional intrinsic currents capable of producing Ca2þ and NMDA spikes in dendrites [17], however we have found comparable results across simulations (not shown here). Adding such features to our model could make it more realistic. However, as we have shown, they are not required to explain the basic phenomena of D- and I-waves. Beyond these issues, our future work will focus on modeling the induction of long-term plasticity with TMS protocols.

References [1] O’Reardon J, Solvason H, Janicak P, et al. Efficacy and safety of transcranial magnetic stimulation in the acute treatment of major depression: a multisite randomized controlled trial. Biol Psychiatry 2007;62:1208e16. [2] Hsu W, Cheng C, Liao K, Lee I, Lin Y. Effects of repetitive transcranial magnetic stimulation on motor functions in patients with stroke: a meta-analysis. Stroke 2012;43:1849e57. [3] Boyd S, Rothwell J, Cowan J, et al. A method of monitoring note on motor conduction velocities. J Neurol Neurosurg Psychiatry 1986;49:251e7. [4] Nakamura H, Kitagawa H, Kawaguchi Y, Tsuji H. Intracortical facilitation and inhibition after transcranial magnetic stimulation in conscious humans. J Physiol (Lond) 1997;498:817e23. [5] Di Lazzaro V, Oliviero A, Meglioni M, et al. Direct demonstration of the effect of lorazepam on the excitability of the human motor cortex. Clin Neurophysiol 2000;111:794e9. [6] Di Lazzaro V, Oliviero A, Saturno E, et al. The effect on corticospinal volleys of reversing the direction of current induced in the motor cortex by transcranial magnetic stimulation. Exp Brain Res 2001;138:268e73. [7] Patton H, Amassian V. Singe- and multiple-unit analysis of cortical stage of pyramidal tract activation. J Neurophysiol 1954;17:345e63. [8] Kernell D, Wu CP. Response of the pyramidal tract to stimulation of the baboon’s motor cortex. J Physiol 1967;191:653e72. [9] Day B, Dressler D, Maertens de Noordhout A, Marsden C, Nakashima K, Thompson JRP. Electric and magnetic stimulation of human motor cortex: surface emg and single motor unit responses. J Physiol 1989;412:449e73. [10] Ziemann U, Rothwell J. I-waves in motor cortex. J Clin Neurophysiol 2000; 17:397e405. [11] Esser S, Hill S, Tononi G. Modeling the effects of transcranial magnetic stimulation on cortical circuits. J Neurophysiol 2005;94:622e39. [12] Beaulieu C, Colonnier M. A laminar analysis of the number of roundasymmetrical and flat-symmetrical synapses of spines, dendritic trunks and cell bodies in area 17 of the cat. J Comp Neurol 1985;231:180e9. [13] Hodgkin A, Huxley A. A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol 1952;117:500e44. [14] Amassian V, Stewart M, Quirk G, Rosenthal J. Physiological basis of motor effects of a transient stimulus to cerebral cortex. Neurosurgery 1987;20:74e93. [15] Hines M, Carnivale N. The neuron simulation environment. Neural Comput 1997;9:1179e209. [16] Hines M, Davidson A, Muller E. Neuron and python. Front Neuroinform 2009;3(1). [17] Larkum M, Nevian T, Sandler M, Polsky A, Schiller J. Synaptic integration in tuft dendrites of layer 5 pyramidal neurons: a new unifying principle. Science 2009;325:756e60. [18] Pospischil M, Toledo-Rodriguez M, Monier C, et al. Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 2008;99:427e41. [19] Bernander O, Koch C, Douglas R. Amplification and linearization of distal synaptic input to cortical pyramidal cells. J Neurophysiol 1994;72(6):2743e53. [20] Traub R, Buhl E, Gloveli T, Whittington M. Fast rhythmic bursting can be induced in layer 2/3 cortical neurons by enhancing persistent Naþ conductance or by blocking BK channels. J Neurophysiol 2003;89:909e21. [21] Evarts E. Relation of pyramidal tract activity to force exerted during voluntary movement. J Neurophysiol 1968;31:14e27. [22] Evarts E, Fromm C, Kroeller J, Jennings V. Motor cortex control of finely graded forces. J Neurophysiol 1983;49:1199e215. [23] Cheney P, Fetz E. Functional classes of primate corticomotoneuronal cells and their relation to active force. J Neurophysiol 1980;44:773e91. [24] Hasenstaub A, Otte S, Callaway E, Sejnowski T. Metabolic cost as a unifying principle governing neuronal biophysics. Proc Natl Acad Sci U S A 2010; 107(27):12329e34.

[25] Rosenthal J, Waller H, Amassian V. An analysis of the activation of motor cortical neurons by surface stimulation. J Neurophysiol 1967;30:844e58. [26] Otis T, Mody I. Differential activation of GABAA and GABAB receptors by spontaneously released transmitter. J Neurophysiol 1992;67:227e35. [27] Stern P, Edwards F, Sakmann B. Fast and slow components of unitary EPSCs on stellate cells elicited by focal stimulation in slices of rat visual cortex. J Physiol 1992;449:247e78. [28] Haeusser M, Roth A. Estimating the time course of the excitatory synaptic conductance in neocortical pyramidal cells using a novel voltage jump method. J Neurosci 1997;17(20):7606e25. [29] Feldmeyer D, Luebke J, Silver R, Sakmann B. Synaptic connections between layer 4 spiny neurone-layer 2/3 pyramidal cell pairs in juvenile rat barrel cortex: physiology and anatomy of interlaminar signalling within a cortical column. J Physiol 2002;538(3):803e22. [30] Bartos M, Vida I, Frotscher M, Geiger J, Jonas P. Rapid signaling at inhibitory synapses in a dentate gyrus interneuron network. J Neurosci 2001;21:2687e98. [31] Otis T, De Koninck Y, Mody I. Characterization of synaptically elicited gabab responses using patch-clamp recordings in rat hippocampal slices. J Physiol 1993;463:391e407. [32] Petreanu L, Mao T, Sternson S, Svoboda K. The subcellular organization of neocortical excitatory connections. Nature 2009;457:1142e5. [33] Markram H, Tsodyks M. Redistribution of synaptic efficacy between neocortical pyramidal neurons. Nature 1996;382:807e10. [34] Abbott L, Varela J, Sen K, Nelson S. Synaptic depression and cortical gain control. Science 1997;275:220e4. [35] Zucker R, Regehr W. Short-term synaptic plasticity. Annu Rev Physiol 2002; 64:355e405. [36] Varela J, Sen K, Gibson J, Fost J, Abbott L, Nelson S. A quantitative description of short-term plasticity at excitatory synapses in layer 2/3 of rat primary visual cortex. J Neurosci 1997;17(20):7926e40. [37] Weyh T, Siebner H. Hirnstimulation-technische Grundlagen. In: Siebner HR, Ziemann U, editors. Das TMS-Buch. Springer; 2007. pp. 17e26. [38] Rothkegel H, Sommer M, Paulus W, Lang N. Impact of pulse duration in single pulse TMS. Clin Neurophysiol 2010;121(11):1915e21. [39] Peterchev A, Jalinous R, Lisanby S. A transcranial magnetic stimulator inducing near-rectangular pulses with controllable pulse width (cTMS). IEEE Trans Biomed Eng 2008;55(1):257e66. [40] Peterchev A, Goetz S, Westin G, Luber B, Lisanby S. Pulse width dependence of motor threshold and input-output curve characterized with controllable pulse parameter transcranial magnetic stimulation. Clin Neurophysiol July 2013;124(7):1364e72. http://dx.doi.org/10.1016/j.clinph.2013.01.011. [41] Pashut T, Wolfus S, Friedman A, et al. Mechanisms of magnetic stimulation of central nervous system neurons. PLoS Comput Biol 2011;7(3):e1002022. [42] Di Lazzaro V, Restuccia D, Oliviero A, et al. Effects of voluntary contraction on descending volleys evoked by transcranial stimulation in conscious humans. J Physiol 1998;508:625e33. [43] Kaneko K, Kawai S, Fuchigami Y, Shiraishi G, Ito T. Effect of stimulus intensity and voluntary contraction on corticospinal potentials following transcranial magnetic stimulation. J Neurol Sci 1996;139:131e6. [44] Davey N, Romaiguere P, Maskill D, Ellaway P. Suppression of voluntary motor activity revealed using transcranial magnetic stimulation of the motor cortex in man. J Physiol (Lond) 1994;477:223e35. [45] Ilic T, Meintzschel F, Cleff U, Ruge D, Kessler K, Ziemann U. Short-interval paired-pulse inhibition and facilitation of human motor cortex: the dimension of stimulus intensity. J Physiol 2002;545(1):153e67. [46] Nakamura H, Kitagawa H, Kawaguchi Y, Tsuji H. Direct and indirect activation of human corticospinal neurons by transcranial magnetic and electrical stimulation. Neurosci Lett 1996;210:45e8. [47] Di Lazzaro V, Oliviero A, Profice P, et al. Comparison of descending volleys evoked by transcranial magnetic and electric stimulation in conscious humans. Electroencephalogr Clin Neurophysiol 1998;109:397e401. [48] Hallermann S, de Kock C, Stuart G, Kole M. State and location dependence of action potential metabolic cost in cortical pyramidal neurons. Nat Neurosci 2012;15:1007e14. [49] Jadi M, Polsky A, Schiller J, Mel B. Location-dependent effects of inhibition on local spiking in pyramidal neuron dendrites. PLoS Comput Biol 2012;8(6):e1002550. [50] Ziemann U. TMS and drugs. Clin Neurophysiol 2004;115(8):1717e29. [51] Paulus W, Classen J, Cohen L, et al. State of the art: pharmacologic effects on cortical excitability measures tested by transcranial magnetic stimulation. Brain Stimul 2008;1(3):151e63. [52] Kang Y, Kaneto T, Ohishi H, Endo K, Araki T. Spatiotemporally differential inhibition of pyramidal cells in the cat motor cortex. J Neurophysiol 1994; 71:280e93. [53] Schoenle P, Isenberg C, Crozier T, Dressler D, Machetanz J, Conrad B. Changes of transcranially evoked motor responses in man by midazolam, a short acting benzodiazepine. Neurosci Lett 1989;101(3):321e4. [54] Palmieri M, Iani C, Scalise A, et al. The effect of benzodiazepines and flumazenil on motor cortical excitability in the human brain. Brain Res 1999; 815(2):192e9. [55] Taniguchi M, Nadstawek J, Langenbach U, Bremer F, Schramm J. Effects of four intravenous anesthetic agents on motor evoked potentials elicited by magnetic transcranial stimulation. Neurosurgery 1993;33:407e15. [56] Ziemann U, Rothwell J, Ridding M. Interaction between intracortical inhibition and facilitation in human motor cortex. J Physiol 1996;496:873e81.

C.V. Rusu et al. / Brain Stimulation 7 (2014) 401e414 [57] Ziemann U, Tergau F, Wassermann E, Wischer S, Hildebrandt J, Paulus W. Demonstration of facilitatory I wave interaction in the human motor cortex by paired transcranial magnetic stimulation. J Physiol 1998;511:181e90. [58] Hanajima R, Ugawa Y, Terao Y, et al. Paired-pulse magnetic stimulation of the human motor cortex: differences among I waves. J Physiol 1998;509:607e18. [59] Stefan K, Kunesch E, Cohen LG, Benecke R, Classen J. Induction of plasticity in the human motor cortex by paired associative stimulation. Brain 2000; 123(3):572e84. [60] Wolters A, Sandbrink F, Schlottmann A, et al. A temporally asymmetric Hebbian rule governing plasticity in the human motor cortex. J Neurophysiol 2003;89(5):2339e45. [61] Rossi S, Hallett M, Rossini P, Pascual-Leone A. Safety, ethical considerations, and application guidelines for the use of transcranial magnetic stimulation in clinical practice and research. Clin Neurophysiol 2009;120(1):2008e39. [62] Mull B, Seyal M. Transcranial magnetic stimulation of left prefrontal cortex impairs working memory. Clin Neurophysiol 2001;112(9):1672e5. [63] Osaka N, Otsuka Y, Hirose N, et al. Transcranial magnetic stimulation (TMS) applied to left dorsolateral prefrontal cortex disrupts verbal working memory performance in humans. Neurosci Lett 2007;418(3):232e5. [64] Kamitani Y, Bhalodia V, Shimojo YKS. A model of magnetic stimulation of neocortical neurons. Neurocomputing 2001;38e40:697e703. [65] Silva S, Basser P, Miranda P. Elucidating the mechanisms and loci of neuronal excitation by transcranial magnetic stimulation using a finite element model of a cortical sulcus. Clin Neurophysiol 2008;119(10):2405e13. [66] Di Lazzaro V, Ziemann U, Lemon R. State of the art: physiology of transcranial motor cortex stimulation. Brain Stimul 2008;1:345e62. [67] Phillips C. Epicortical electrical mapping of motor areas in primates. In: Bock G, O’Connor M, Marsh J, editors. Motor Areas of Cerebral cortex. John Wiley; 1987. pp. 5e20. [68] Kujirai T, Caramia M, Rothwell J, et al. Corticocortical inhibition in human motor cortex. J Physiol 1993;471:501e19. [69] Ravazzani P, Ruohonen J, Grandori F, Tognola G. Magnetic stimulation of the nervous system: induced electric field in unbounded, semi-infinite, spherical, and cylindrical media. Ann Biomed Eng 1996;24:606e16. [70] Thielscher A, Kammer T. Linking physics with physiology in tms: a sphere field model to determine the cortical stimulation site in tms. Neuroimage 2002;17: 1117e30. [71] Manola L, Holsheimer J, Veltink P, Buitenweg J. Anodal vs cathodal stimulation of motor cortex: a modeling study. Clin Neurophysiol 2007;118:464e74. [72] Markram H, Toledo-Rodriguez M, Wang Y, Gupta A, Silberberg G, Wu C. Interneurons of the neocortical inhibitory system. Nat Rev Neurosci 2004;5(10): 793e807. [73] Wilson N, Runyan C, Wang F, Sur M. Division and subtraction by distinct cortical inhibitory networks in vivo. Nature 2012;488:343e8. [74] Hille B. Ionic Channels of Excitable Membranes. Sunderland, MA: Sinauer; 2001.

Sodium current

INa ¼ g Na m3 h ðU  ENa Þ dm ¼ am ðUÞ ð1  mÞ  m bm ðUÞ dt dh ¼ ah ðUÞ ð1  hÞ  h bh ðUÞ dt 0:32 ðU  UT  13 mVÞ 1 am ¼ exp½  ðU  UT  13 mVÞ=4 mV  1 ms

bm ¼

L2/3 Cell currents

0:28 ðU  UT  40 mVÞ 1 exp½ðU  UT  40 mVÞ=5 mV  1 ms

ah ¼ 0:128 exp½  ðU  UT  17 mVÞ=18 mV bh ¼

S cm2

UT ¼ 55 mV

IKd ¼ g Kd n4 ðU  EK Þ dn ¼ an ðUÞ ð1  nÞ  n bn ðUÞ dt 0:032 ðU  UT  15 mVÞ 1 an ¼ exp½  ðU  UT  15 mVÞ=5 mV  1 ms

bn ¼ 0:5 exp½  ðU  UT  10 mVÞ=40 mV

1 ms

UT ¼ 55 mV EK ¼ 100 mV 8 S > > for the excitatory L2=3 cells < 0:005 cm2 ¼ > > : 0:01 S for the inhibitory L2=3 cells cm2

Slow non-inactivating potassium current, RS cells only

IM ¼ g M p ðU  EK Þ  dp ¼ ðpN ðUÞ  pÞ sp ðUÞ dt

(7)

where Ex is the reversal potential of the ion involved and m and h are the gate activation and inactivation variables, respectively. The gating variables m and h obeyed first-order dynamics:

dx 1 ¼  ½x  xN ðUÞ; sx ðUÞ dt

g Na ¼ 0:05

Delayed-rectifier potassium current

Channel kinetics were modeled such that each voltagedependent ion current follows:

Ix ¼ g x mM hH ðU  Ex Þ;

1 ms

4 1 1 þ exp½  ðU  UT  40 mVÞ=5 mV ms

ENa ¼ 50 mV

g Kd Appendix

413

(8)

pN ðUÞ ¼

sp ðUÞ ¼

1 1 þ exp½  ðU þ 35 mVÞ=10 mV

smax

3:3 exp½ðUþ35 mVÞ=20 mVþexp½ðUþ35 mVÞ=20 mV

EK ¼ 100 mV

g M ¼ 7,105

S cm2

smax ¼ 1000 ms

with

sx ðUÞ ¼ ½ax ðUÞ þ bx ðUÞ1 :

(9)

The functions a and b are fitted to match experimental data. The corresponding steady-state variables xN depend on the membrane potential:

ax ðUÞ : xN ðUÞ ¼ ax ðUÞ þ bx ðUÞ Details of individual currents are given below.

(10)

L5 cell currents Channel kinetics were modeled such that each voltagedependent ion current follows:

Ix ¼ g x mM hH ðU  Ex Þ;

(11)

where Ex is the reversal potential of the ion involved and m and h are the gate activation and inactivation variables, respectively. The gating variables m and h obeyed first-order dynamics:

414

C.V. Rusu et al. / Brain Stimulation 7 (2014) 401e414

Potassium current IK

dx 1 ¼  ½x  xN ðUÞ: sx dt

(12)

In the case of the currents INa and IK, the activation and inactivation time constants do not depend on the membrane potential. The voltage dependence of the IKDR time constant is described in “Delayed-rectifier potassium current IKDR , axon only” further down.

Sodium current INa

INa ¼ g Na h s m3 ðU  ENa Þ 1   U þ 45 mV 1 þ exp 3 mV

hN ¼

sN ¼

1   U þ 44 mV 1 þ exp 3 mV

IK ¼ g K n3 ðU  EK Þ þ g K2 n32 ðU  EK Þ 1   U þ 40 mV 1 þ exp 3 mV

nN ¼

1   U þ 40 mV 1 þ exp 3 mV 8 0 for the dendrites > > > > > > > S < for the soma 0:06 gK ¼ cm2 > > > > > > S > : 1 for the axon cm2 8 S > > > > 0:0001 cm2 for the dendrites > > > < S gK2 ¼ > for the soma 0:3 > 2 > cm > > > > : 0 for the axon

n2N ¼

EK ¼ 80 mV

sn ¼ 0:5 ms sn2 ¼ 10 ms

Delayed-rectifier potassium current IKDR , axon only

1   U þ 40 mV 1 þ exp 3 mV 8 S > > for the 0:004 > > cm2 > > > > > < S ¼ for the 0:008 2 > cm > > > > > > > > : 1:9 S for the cm2 ¼ 40 mV sh ¼ 0:1 ms

IKDR ¼ g KDR m4 ðU  EK Þ

mN ¼

gNa

ENa

1   U þ 29:5 mV 1 þ exp 10 mV   8 U þ 10 mV > > 0:25 ms þ 4:35 exp ms; if U < 10 mV > > 10 mV <

mN ¼ dendrites soma axon

ss ¼ 50 ms sm ¼ 0:05 ms

sm ¼

  > > U  10 mV > > : 0:25 ms þ 4:35 exp ms; 10 mV

EK ¼ 80 mV

g KDR ¼ 10

S cm2

otherwise