Transcranial magnetic stimulation (TMS) is commonly used to activate or inactivate specific cortical areas in a noninvasive manner. Because of technical constraints, the precise effects of TMS on cortical circuits are difficult to assess experimentally. Here, this issue is investigated by constructing a detailed model of a portion of the thalamocortical system and examining the effects of the simulated delivery of a TMS pulse. The model, which incorporates a large number of physiological and anatomical constraints, includes 33,000 spiking neurons arranged in a 3-layered motor cortex and over 5 million intra- and interlayer synaptic connections. The model was validated by reproducing several results from the experimental literature. These include the frequency, timing, dose response, and pharmacological modulation of epidurally recorded responses to TMS (the so-called I-waves), as well as paired-pulse response curves consistent with data from several experimental studies. The modeled responses to simulated TMS pulses in different experimental paradigms provide a detailed, self-consistent account of the neural and synaptic activities evoked by TMS within prototypical cortical circuits.
Transcranial magnetic stimulation (TMS) is a technique frequently used for the noninvasive stimulation of the human cerebral cortex. Depending on stimulation parameters, TMS can disrupt or induce activity in a particular cortical region, thereby providing a valuable tool for the study of motor and sensory processes, attention, memory, language, and neuronal plasticity (Di Lazzaro et al. 2004; Lisanby et al. 2000; Terao and Ugawa 2002; Ziemann and Rothwell 2000). Moreover, TMS is being actively investigated in the clinic as a promising treatment modality for several psychiatric and neurological disorders (Kobayashi and Pascual-Leone 2003).
Despite its increasing use in both research and clinical settings, the mechanisms of TMS action on the cerebral cortex are still unclear. Studying such mechanisms is hampered by several technical constraints. Developing paradigms for animal experiments is complicated by difficulties in designing small TMS coils (Barker 2002). Additionally, the strong magnetic field induces severe recording artifacts, reducing the ability to evaluate its effects with sufficient spatial and temporal resolution (Virtanen et al. 1999).
Researchers have nevertheless been able to obtain some mechanistic insights by developing sophisticated approaches that are applicable to humans. Specifically, the direct output of motor cortical neurons (called I-waves) has been recorded by implanting epidural electrodes along the motor pathway (Di Lazzaro et al. 1998a,b, 2000, 2001b; Edgley et al. 1997). Furthermore, by applying closely spaced TMS pulses (paired-pulse paradigms) it has been possible to make inferences about the cortical circuitry involved in the response to TMS (Di Lazzaro et al. 1998c; Fisher et al. 2002; Hanajima et al. 1998, 2003; Roshan et al. 2003; Ziemann et al. 1998).
To gain a better understanding of the cortical response to TMS it would be useful to develop a sufficiently detailed model of the effects of TMS at the cellular and circuit level. For instance, does a TMS pulse exert preferential effects on different cortical layers and synaptic circuits? A way to address these questions is to explicitly model the effects of a TMS pulse on the cerebral cortex by using large-scale computer simulations that incorporate a considerable number of anatomical and physiological facts.
Building on previous work (Hill and Tononi 2005; Lumer et al. 1997), we constructed a large-scale model of a prototypical thalamocortical circuit and perturbed it with one or more simulated TMS pulses. The model contains over 33,000 spiking neurons, organized into a 3-layered cortical area and the associated thalamic and reticular regions. The circuitry includes an anatomically based set of intralayer, interlayer, thalamocortical, corticothalamic, and interareal connections. Each model neuron is endowed with a set of intrinsic as well as synaptic currents [α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA), N-methyl-d-aspartate (NMDA), γ-aminobutyric acid-A (GABAA), and γ-aminobutyric acid-B (GABAB)]. To facilitate comparison with the experimental literature, we chose to simulate a portion of motor cortex because it is the most frequently studied area in the TMS literature. Furthermore, the effects of TMS examined in this paper are believed to occur locally within the cortical circuit (Ziemann and Rothwell 2000).
To constrain the simulations, we aimed at replicating, within a single model and with the same set of parameters, several independent results from the experimental literature. Specifically, the model reproduces the frequency, timing, dose response, and behavioral/pharmacological modulation of I-waves. The same model also reproduces paired-pulse response curves and the response timing difference observed when the orientation of the coil is reversed. In doing so, the model overcomes the technical constraints inherent to TMS experimental methods to provide a detailed picture of the synaptic, neural, and circuit activity that leads to these responses. In this way, the simulations presented below offer a detailed self-consistent account of how cortical circuitry responds to TMS pulses.
A model of the motor thalamocortical system
Simulations were performed using a model of motor cortex that was based on previous large-scale models of the visual thalamocortical system (Hill and Tononi 2005; Lumer et al. 1997). For the present purposes, the model was adapted to reflect the physiology and anatomy of motor cortex. The level of detail chosen for the model was based on a desire to study the effects of TMS on cortical circuits in a realistic, but computationally feasible, manner.
The model (Fig. 1) contains over 33,000 spiking integrate-and-fire neurons and over 5 million connections. The focus of the model is a primary motor cortical area (MP) built to reflect a 5.5 × 5.5-mm portion of primary motor cortex (area 4). MP is subdivided into 3 distinct layers—L2/3, L5, and L6—corresponding to the superficial layers (layers 2 and 3), layer 5, and layer 6, respectively. Layer 1 is not modeled because of the sparseness of neurons there and layer 4 is not included because it is very thin in the motor cortex (Gatter et al. 1978). Each layer has a distinct set of intralayer and interlayer connections.
MP receives input from thalamus (T), specifically from motor thalamic subdivisions. MP also receives input from a somatosensory area (SI), corresponding to areas 1, 2, 3a, and 3b, and input from a premotor area (PM), corresponding to area 6. These inputs reflect the major sources of input to primary motor cortex (Ghosh et al. 1987; Stepniewska et al. 1994a,b). These areas are connected to MP with distinct corticothalamic, thalamocortical, and corticocortical connections.
Excitatory and inhibitory neurons are modeled as single-compartment integrate-and-fire units. Neurotransmission is modeled by excitatory voltage-independent AMPA channels and voltage-dependent NMDA channels as well as inhibitory GABAA and GABAB channels. The effect of a TMS pulse is modeled as the synchronous activation of cortical synapses. Below, we provide a detailed look at the modeled regions, connections, neurons, synaptic channels, and TMS activation.
Organization of model regions
PRIMARY MOTOR CORTICAL AREA.
Anatomically, cells in the primary motor cortex (area 4) form columns, about 300 μm in diameter, separated by 50-μm cell-sparse regions (Meyer 1987; Mountcastle 1998). Functionally, the firing rate of a given cell in motor cortex increases when movements are made in a specific direction. This direction is referred to as the cells' preferred direction. Neurons within a column typically share a single preferred direction and adjacent columns cover a range of preferred directions. Cells of opposite preferred direction are typically about 350 μm apart (Amirikian and Georgopoulos 2003), close enough to be in adjacent columns. In terms of connections, interlayer projections arising from small (<200 μm diameter) spots in motor cortex terminate in discrete patches with center-to-center distances ranging from 450 to 1,550 μm (Lund et al. 1993), close enough to range from 2 to several columns away. Columns of similar function are often connected in cortex (Mountcastle 1998). Together, this suggests that columns of similar preferred direction are typically a column's width apart from each another, separated by columns of differing (including opposite) preferred directions.
The construction of MP followed this organization. A topographic element was defined to contain 2 excitatory and one inhibitory cell. The neurons in a given topographic element were located adjacent in space and had the same preferred direction. Columns were formed from one 5 × 5 group of topographic elements from each layer (for a total of 225 cells per column). MP consisted of an 8 × 8 array of these columns interspersed with an 8 × 8 array of columns having the opposite preferred direction. Based on this organization and the above literature on columnar sizes, in MP the diameter of columns was 300 μm, the edge-to-edge distance between adjacent columns with opposite preferred direction was 50 μm, and the edge-to-edge distance between columns with the same preferred direction was 400 μm.
MP received input from PM, SI, and T. PM and SI were modeled as 40 × 40 grids of randomly firing neurons. The mean firing rate of neurons in SI and PM was adjusted to 1 Hz to match spontaneous firing rates found in intercortical cells in somatosensory and motor cortex (Swadlow 1994, 1990).
A number of thalamic subdivisions connect to motor cortex, including the ventroanterior–ventrolateral complex and the ventral posterolateral area. Connections between the thalamus and motor cortex have a basic topographical organization (Stepniewska et al. 1994b). Microstimulation of motor thalamus evokes movements focused around single joints (Vitek et al. 1996), suggesting that individual cells in motor thalamus may approximately correspond to specific preferred directions. Thalamic cells have a typical spontaneous firing rate of 10–20 Hz (Farley 1997; Raeva et al. 1999). The thalamus receives inhibitory input from the reticular nuclear and local inhibitory interneurons that help to modulate its activity (Ando et al. 1995).
For the purpose of the present study, we simulated a motor thalamus, T, which included two 40 × 40 grids of cells, one excitatory and one inhibitory, that corresponded topographically with MP. Half of the excitatory cells in T had a preferred direction of 0°, whereas the other half had a preferred direction of 180°. Random input, mimicking spontaneous activity from motor-related structures such as the cerebellum and other motor cortical areas, caused cells in T to fire with a rate of 10–20 Hz. A reticular nucleus (R) was modeled as a 40 × 40 grid of inhibitory cells that corresponded topographically to the grids in T. It should be noted that in the model, T was an important part of the neural circuit and helped to generate realistic spontaneous activity. However, the thalamus is generally believed to play a minor role in the generation of early responses to TMS (Ziemann and Rothwell 2000). Therefore the discussion of the thalamus in the remainder of the paper will be minimal.
The connectivity of the model is described in Table 1. The table details the layer of origin and cell type as well as the target layer, cell type, and synapse type of a given modeled connection. The table describes the peak height, width, and maximum radius of the Gaussian probability distribution used to generate the connections as well as the strength and conduction delays for the given connection. All of these parameters were constrained by data available in the literature. In what follows, we will briefly discuss the organization of the major connection pathways. Except where noted, excitatory connections are mediated by AMPA, whereas inhibitory connections are mediated by GABAA. Specific parameter values are available in the table.
HORIZONTAL INTRALAYER EXCITATORY CONNECTIONS.
Excitatory neurons in the superficial layers as well as layer 5 of motor cortex have intralayer horizontal projections that may extend with a radius of ≤2–3 mm, whereas layer 6 cells appear to give out few if any axon collaterals (Aroniadou and Keller 1993; Gatter et al. 1978; Ghosh and Porter 1988; Yamashita and Arikuni 2001). Over these distances, neurons that elicit responses from agonistic muscles have reciprocal excitatory connections, whereas neurons eliciting responses from antagonistic muscles do not (Capaday et al. 1998). Superficial neurons synapse onto other superficial neurons about twice as frequently as layer 5 neurons onto other layer 5 neurons (Thomson et al. 2002). Horizontal superficial and layer 5 connections have a strong NMDA-receptor–mediated component as well as an AMPA-receptor–mediated component (Aroniadou and Keller 1993).
To model this arrangement the intralayer excitatory connections of L6 were sparse and over a short radius, whereas intralayer excitatory connections in L2/3 and L5 were spread out with a radius corresponding to about 2 mm. Twice as many synapses were generated between L2/3 cells as between L5 cells. Both L2/3 and L5 intralayer connections had synapses mediated by AMPA and NMDA receptors. Excitatory cells with the same preferred direction were directly interconnected, whereas those with the opposite preferred direction were not.
VERTICAL INTERLAYER EXCITATORY CONNECTIONS.
Excitatory vertical connections in motor cortex project from one layer to another with a divergence of about 250 μm (Gatter et al. 1978; Kaneko et al. 1994, 2000). The projections from superficial layer neurons make connections onto layer 5 neurons 4 times more frequently than onto layer 6 neurons (Kaneko et al. 2000). Connections from layer 5 in motor cortex are primarily constrained to layer 5 and layer 6 (Ghosh and Porter 1988; Ghosh et al. 1988), although intracortical microstimulation of layer 5 is able to elicit responses in the superficial layers (Aroniadou and Keller 1993; Asanuma and Rosen 1973). Projections from layer 6 appear to be confined to deep cortical layers in motor cortex (Noda and Yamamoto 1984).
In MP, neurons sent vertical interlayer projections with a narrow columnar-like distribution (i.e., with a limited horizontal spread but extensive vertical range). The projections from L2/3 to L5 and L6 preserved the 4/1 ratio observed in motor cortex. The projections from L5 had relatively weak connections to L2/3 compared with projections to L5 and L6. L6 sent weak vertical projections to L5.
INTRACORTICAL INHIBITORY CONNECTIONS.
Inhibitory interneurons in cerebral cortex can have synaptic interactions mediated by fast GABAA receptors or slow GABAB receptors (Kang et al. 1994). GABAA-mediated inhibitory cells have primarily intralayer projections, whereas GABAB-mediated inhibitory cells are found primarily in superficial layers and project in a columnar fashion to most layers (Conde et al. 1994; Kawaguchi 1995; Porter et al. 2000). About 15% of cells in motor cortex are inhibitory cells, with deeper layers typically having fewer inhibitory cells than superficial layers (Beaulieu 1993). About 10–20% of synapses in cortex are inhibitory (Beaulieu and Colonnier 1985). Layers 5 and 6 of motor cortex have a GABAA receptor density that is about half of what it is in layers 2 and 3 (Schwark et al. 1994). Inhibitory interneurons typically form intralayer connections at the same to twice the frequency of excitatory intralayer connections (Thomson et al. 2002). The radius of the GABAA-mediated connections in motor cortex is typically about 1 mm, whereas GABAB-mediated connections have a columnar radius of about 450 μm (Kang et al. 1994).
In MP, GABAA-like inhibitory cells were modeled in each layer and had horizontal intralayer projections. For simplicity, GABAB-like connections were modeled by giving excitatory cells in L2/3 narrow GABAB-like projections to all 3 layers. Inhibitory connections formed 10–20% of all cortical synapses. Inhibitory connections were one to 2 times as common as intralayer excitatory connections. Twice as many GABAA-mediated connections were generated in L2/3 as in L5. GABAA connections had a radius of 1 mm and GABAB connections had a radius of 450 μm.
Input from other cortical areas
Motor cortical afferents from somatosensory cortex terminate primarily in layer 3 at about the same frequency as thalamocortical projections (Ichikawa et al. 1985; Kosar et al. 1985), whereas afferents from area 6 terminate strongly in layer 5 compared with their terminations in layer 3 (Tokuno and Nambu 2000). AMPA and NMDA receptors mediate the response of motor cortical afferents from both somatosensory cortex and area 6 (Shima and Tanji 1998). Area 6 has about 25% more cells than somatosensory cortex that project to motor cortex. Additionally, retrograde tracer studies have shown that cells outside of motor cortex account for half of all cells connecting to a given portion of motor cortex (Ghosh et al. 1987).
In the model, SI projected to L2/3 with the same frequency as the thalamocortical projections. PM sent strong projections to L5 compared with its projections to L2/3, and sent out 25% more projections than SI. Inputs from SI and PM were both mediated by AMPA and NMDA synapses. SI, PM, and thalamocortical projections accounted for about half of all connections to a given portion of MP.
Excitatory neurons in the motor thalamus send projections to the motor cortex and the reticular nucleus, but do not send local collaterals (Aumann et al. 1998; Liu et al. 1995). Inhibitory neurons send projections within their immediate vicinity (Scheibel and Scheibel 1966). The reticular nucleus sends local collaterals within the reticular nucleus as well as focal topographically organized projections to the related thalamic area (Liu et al. 1995).
In the model's thalamic area, the excitatory population of T sent focused projections to R, whereas the inhibitory population sent collaterals locally within T. R sent moderately focused projections to both populations of T.
Projections from the thalamus to motor cortex have a slight preference for superficial layers over layer 5 and terminate in a number of discrete patches, each typically 0.5 mm in diameter and spread over an area about 2–4 mm across (Aumann et al. 1998; Shinoda et al. 1993). A coherent functional map in the motor thalamus can be observed with microstimulation, which typically evokes a response involving a single joint (Vitek et al. 1996). This suggests that thalamic cells have preferential access to one preferred direction over another, which is feasible because of the patchy nature of connections formed by individual thalamocortical cells.
In the model, excitatory cells in T sent projections to populations in L2/3, L5, and L6 that had the same preferred direction as that of the thalamic cell. Connections to L5 were strong, to L2/3 were moderate, and to L6 were weak. Projections from a single cell covered a large diameter in MP.
Corticothalamic projections originate in layers 5 and 6 of motor cortex, although the projections from layer 6 are far more numerous. Layer 6 neurons send projections to topographically corresponding neurons in the thalamus and reticular nucleus (Kakei et al. 2001). Corticothalamic projections are reciprocal with thalamocortical connections, but have a wider radius (Ando et al. 1995; McFarland and Haber 2002). Corticothalamic projections have been shown to elicit monosynaptic responses from inhibitory neurons in the thalamus and reticular nucleus (Ando et al. 1995).
In the model, MP and T were reciprocally connected, but with a wider connection radius for connections from MP to T than vice versa. Because of the vastly larger number of projections originating from layer 6 than from layer 5, connections from MP to T were simplified as originating only from L6. L6 also sent projections to the corresponding topographical location in R and to the inhibitory population of T.
We define transmission delays as the time between the emergence of a presynaptic neuronal spike and the moment postsynaptic current integration begins, thus including axonal conduction time, synaptic delay, and dendritic conduction time. The transmission delay between layer 3 and layer 5 pyramidal cells is 1.4 ms (Thomson et al. 2002). Information on other connections in cortex is not as readily available but can be calculated indirectly from the following information. Synaptic delays are about 0.4 ms (Eccles 1964). Axonal conduction velocities of short-range excitatory projections in motor cortex are about 1 m/s (Swadlow 1994) and those of inhibitory interneurons are about 0.4 m/s (Kang et al. 1994). Finally, in human motor cortex, the distance between the centers of layers 5 and 6 is about 0.5 mm (Rockel et al. 1980). Transmission delays from thalamus to cortex average about 3 ms (Noda and Yamamoto 1984). Corticothalamic delays also average about 3 ms, whereas delays of cortical projections to the reticular nucleus average about 2.5 ms. Thalamic interneurons act over a very short time interval (<1 ms), as do connections between the reticular nucleus and thalamus (Ando et al. 1995).
In the model, all conduction times were generated using a random Gaussian distribution. Specific mean conduction times were determined from the direct experimental measurements given above, when available. If direct measures were not available, the mean conduction delays were calculated from the synaptic delay, conduction velocity, and cell-to-cell distances. Because PM and SI fire independently of network activity, accurate measures of connection times from interareal input areas were not needed and were therefore only approximated. See Table 1 for the details of this parameterization.
Neurons were modeled as in Hill and Tononi 2005, with minor modifications. Briefly, a single-compartment leaky integrate-and-fire (IandF) neuron including Hodgkin–Huxley-like currents was used to simulate excitatory and inhibitory neurons (Koch and Segev 1998). To simulate the subthreshold dynamics of voltage-gated sodium channels, the firing threshold (θ) varied based on a resting threshold (θeq), a constant C (0.85), voltage (V), and a threshold time constant (tq) through the equation [modified from Bradley (1961); Hill (1936); and MacGregor and Oliver (1974)] The membrane potential (V) changed for each neuron as where gNa(leak) is the sodium leak, ENa is the sodium reversal potential (30 mV), gK(leak) is the potassium leak, EK is the potassium reversal potential (−90 mV), Isyn is the synaptic current, Iint is the intrinsic current, τm is the membrane time constant, gspike is 1 during an action potential and 0 when a cell is not spiking, and τspike is the time constant during a spike. Spikes lasted tspike.
Table 2 provides the values used for these properties in the model. L5 cells are considered separately from other MP excitatory cells because they were modeled after the large pyramidal tract cells in layer 5. Values for θeq and τm in MP were based on in vivo studies in the cat, with values for τm used that reflect behavior after current applied for <50 ms (Baranyi et al. 1993b). The spiking characteristics of the modeled cells in MP, defined by tq, τspike, and tspike, were set so the modeled cells matched traces of cellular membrane potential during a spike (Baranyi et al. 1993b) as well as firing rate characteristics of the cells (Baranyi et al. 1993a; Chen et al. 1996). Values for gNa(leak) and gK(leak) were selected so that during spontaneous activity modeled cells reached observed resting membrane potentials (Baranyi et al. 1993b). Because of similar firing properties of thalamic cells and inhibitory interneurons in cortex (Baranyi et al. 1993a; Sawyer et al. 1994), these 2 cell types were modeled with similar properties.
Model synaptic channels
Synaptic input, Isyn, constituted the sum of all synaptic channel currents The amplitude and time course of a postsynaptic potential was determined by the conductance of each afferent i, on each channel j, as well as that channel's reversal potential Ej. Channel conductance g(t) changed in response to synaptic activity according to Rise and decay time constants are given by τ1…τm, respectively, and the time to peak tpeak as Conductances were adimensional because the volume of the modeled neurons was not specifically defined. Table 3 lists time constants and reversal potentials for each channel as derived from the literature (Otis and Mody 1992; Otis et al. 1993; Stern et al. 1992).
Single cell to single cell connections in cortex typically show primarily AMPA (accounting for 95% of response peak) mediated excitatory postsynaptic potentials of about 1.3 mV and GABAA-mediated inhibitory postsynaptic potentials of about 1.2 mV in layer 5 pyramidal cells having a resting membrane potential of −60 mV (Markram et al. 1997; Thomson and Bannister 1998; Thomson et al. 1996), with similar inhibitory postsynaptic potential values seen in the more superficial layers (Tamas et al. 1997). Accounting for differing reversal potentials of the various synaptic channels involved as well as the voltage sensitivity of the NMDA receptor, these values were used to establish the values of, and ratios between, the peak conductances involved in the model (Table 3). The peak conductance for GABAB was set so that GABAB had a weak modulatory role in the model. Research on rats has indicated that excitatory connections between excitatory and inhibitory cells are much stronger than between excitatory cells (Holmgren et al. 2003), which may be an important factor in creating balanced network activity. To ensure the proper ratio of excitation and inhibition in the model, the conductance of AMPA synapses onto inhibitory cells was raised to the needs of balancing the network. This resulted in an AMPA conductance that was 75% larger in inhibitory cells than in excitatory cells.
NMDA channel activation was expressed in the model as A dual exponential impulse response modeled ḡNMDA(t), whereas m(V) was modeled as the sum of a fast and slow exponential that modulates the voltage-dependent change in NMDA conductance. Whereas blocking of the NMDA channel occurs almost instantly (<0.06 μs), processes occurring on both a slow and a fast timescale mediate its unblocking (Vargas-Caballero and Robinson 2003). This was used as the basis for modeling an unblocking process with a fast unblock time (about 1 ms) and a second process with a slow unblock time (about 20 ms).
All synapses were designed to saturate to a certain percentage from a single pulse, rendering further responses in a short time span less effective. A single NMDA pulse saturated 54% of a single receptor and an AMPA pulse saturates 38% of a single receptor, reflecting modeling work done on these synapses (Franks et al. 2002). A single pulse on a GABAA channel saturates 60% of a single receptor (Grabauskas and Bradley 2003).
The role of short-term synaptic dynamics in the cortical response to TMS is currently unclear in the experimental literature. However, we found that the inclusion of synaptic depression in our model synapses did not significantly affect our modeled TMS results. For simplicity, we therefore do not include synaptic depression in the simulations presented in this paper.
Intrinsic ion channels
The model neurons were variously equipped with noninactivating hyperpolarization-activated cation current, low-threshold calcium current, persistent sodium current, and depolarization-activated potassium current as in Hill and Tononi (2005). These intrinsic channels are important for the transition between sleeping and waking states. We tested the role of these channels in generating all of our simulated results and found that they did not play a significant role in our simulations. As such, they will not be discussed further.
Simulating TMS pulses
A TMS pulse is a transient magnetic field generated by running current through a wrapped wire in the stimulator coil. The magnetic field from TMS decays rapidly with distance from the cortex, so it is believed to affect cortical but not thalamic targets (Bohning et al. 1997). The effect of TMS on cortex is typically 10 to 20 mm in diameter, with the peak effect nearly constant within 10 mm, allowing TMS to target neural groups in motor cortex that control specific muscles (Jalinous 1995; Ravazzani et al. 1996; Thielscher and Kammer 2002). Experimental studies show that peripheral nerves are most easily stimulated at terminations (Maccabee et al. 1993). Nerve modeling studies support this experimental work by suggesting that a neuron is most likely to be activated by TMS where geometric discontinuities occur, such as fiber terminations (e.g., synapses or dendrites). Modeling studies also suggest that the bends of cortical axons may be sites for TMS activation, whereas dendritic activation is not likely to occur at the stimulus strengths commonly used (Abdeen and Stuchly 1994; Nagarajan et al. 1993). The orientation of fibers or terminations relative to the magnetic field is important in determining response threshold (Abdeen and Stuchly 1994; Maccabee et al. 1993). However, we assumed that the orientations of fiber terminals in all populations vary randomly so that fiber orientation is not a factor in distinguishing the direct response of specific populations to TMS. The majority of current induced in tissue by TMS is present for <100 μs (Jalinous 1995), whereas the fibers responding to TMS are in a strong refractory state 1 ms after TMS, which fades away by 2 ms after TMS (Fisher et al. 2002; Hanajima et al. 2003; Roshan et al. 2003).
Based on the preceding information about how TMS effects cortex, we modeled TMS as a homogeneous activation of fiber terminals within a region of MP 5.5 mm in diameter (i.e., the entire modeled cortical area), but not affecting thalamic areas, PM or SI, because they are outside of the range of TMS effect. For simplicity, because there would be little difference between TMS activation at axonal bends and at fiber terminals in motor cortex except for a small disparity in response latencies, we modeled a TMS pulse as synchronously activating a chosen fraction of fiber terminals throughout MP. For a given cell, a single TMS pulse was therefore modeled as transiently changing the conductance g for a specified fraction of that cell's synaptic channels. The dynamics of this change follow the equations and parameters specified in the section Model synaptic channels. Because the orientation of a given fiber affects its activation threshold, the natural variation of fiber terminal orientations in cortex means that a stronger TMS pulse will activate more fiber terminals than a weak pulse. For this reason, we modeled stronger or weaker TMS pulses by increasing or decreasing the percentage of fiber terminals activated. Based on the short duration of a TMS pulse, we assumed a TMS pulse occurs during a single time step. After this activation, the fibers entered an absolute refractory period during which they could not be activated again for 1 ms, followed by a relative refractory period during which the number of fibers that could be activated increases in a linear fashion from 0 to 100% between 1 and 2 ms after the TMS pulse. This simulation of TMS was intended to reflect the use of a focal coil that induces posterior to anterior (PA) current in the brain. Other coil orientations can induce different effects in the brain, as will be discussed later in the paper.
TMS is believed to elicit firing from pyramidal tract fibers through both direct and indirect activation. Direct activation is the result of the TMS pulse exciting the pyramidal tract fibers themselves, causing spiking. Indirect activation is the result of TMS causing cortical fibers to fire, which causes a perturbation in cortical circuits that leads to the spiking of pyramidal tract fibers. Because the focus of this paper is to explore how TMS affects cortical circuitry, and because TMS stimulation frequently does not induce direct activation, we did not explicitly model direct activation. For reference purposes, the moment of direct activation would correspond to the moment when the TMS pulse occurs (fiber activation).
All simulations were performed using a general-purpose object-oriented interactive neural simulator called Synthesis, written by S. Hill (www.infinitedegrees.info). Differential equations were numerically integrated using the Runge–Kutta 4th-order method (Press 1992) with a step size of 0.1 ms. The simulations were carried out on a dual-processor 1.0-GHz Apple PowerPC G4 machine running Mac OS 10.2 and equipped with 1.5 GB RAM. Each simulation of the model used nearly 1.7 GB of memory. Computational performance varied based on average firing rates, but was typically about 1,000 ms/h during spontaneous activity.
The model, incorporating basic aspects of the anatomy and physiology of motor cortex, showed patterns of spontaneous activity resembling those recorded from awake animals
We chose to model motor cortex because of the large amount of available data from TMS experiments in this area. Figure 1 illustrates the basic anatomical features of the motor cortical system that are incorporated in the model to the extent needed to accurately reproduce experimental results. These features include a complex array of interconnections, a columnar structure, and excitatory and inhibitory model cells with properties based on the literature of cortical neurons of the corresponding type and layer. Further details are available in the methods section.
The model neurons exhibited firing comparable to spontaneous activity in experimental recordings. The general firing pattern of the model is demonstrated in a membrane potential plot diagram (Fig. 2A) showing the activity of 50 representative, adjacent cells chosen from each layer (about 1% of total cells for cortex and about 3% for the thalamus). The membrane potential plot reveals that the spontaneous activity in the model was even, except in L6 where activity was very low. The mean firing rate of excitatory neurons was 1.3 ± 1.5 Hz in L2/3 (mean ± SD), 5.6 ± 3.7 Hz in L5, 0.5 ± 1.0 Hz in L6, and 11.5 ± 9.1 Hz in T (Fig. 2B). The membrane potentials of the modeled neurons were also examined over this period. The mean membrane potential was −66.1 ± 2.2 mV in L2/3, −64.6 ± 2.9 mV in L5, −68.9 ± 2.0 mV in L6, and −66.3 ± 4.9 mV in T (Fig. 2C).
The model firing rates were consistent with values recoded in vivo of about 0.4 Hz in layer 6 cells, 6 Hz in layer 5 cells (Beloozerova et al. 2003), and 10–20 Hz in the thalamus (Farley 1997; Raeva et al. 1999). Simulated membrane potentials were consistent with observations of in vivo motor cortical cellular membrane potentials of 64–66 mV (SD of 2–5 mV) (Baranyi et al. 1993a).
Although the model contains neurons with opposing preferred directions that produce differential responses, we did not explore these properties further in this paper.
The application of a simulated TMS pulse produced a sequence of I-waves that resemble those observed during in vivo experiments
The motor cortical system is known to respond to stimulation with waves of activity that have been recorded from the human epidural space. The electrodes used in these recordings are typically positioned in the epidural space such that a compound action potential is recorded from the axons descending from motor cortex, which originate in layer 5 (Di Lazzaro et al. 1998a, 2001b; Edgley et al. 1997; Fujiki et al. 1996). Comparisons of single-unit recordings made from animals with multiunit epidural recordings confirmed that the epidural response reflects axonal activity (Edgley et al. 1997). To simulate an epidural recording we convolved a vector, whose elements contained the total number of L5 spikes in a given time step, with the shape of an axonal spike derived from layer 5 axonal recordings (Stuart et al. 1997). This convolution creates a measure of the electrical activity that would be recorded in the epidural space after spiking in motor cortex. The convolution did not alter the basic characteristics (timing and relative amplitude) of L5 spiking activity.
The waves of the in vivo epidural response have peak-to-peak intervals of around 1.4 ms and commonly occur in groups of 3 or 4 (Di Lazzaro et al. 1998a, 2001b; Edgley et al. 1997; Fujiki et al. 1996). As a first measure of how our model system responds to stimulation, we simulated a TMS pulse sufficient to induce a strong in vivo–like TMS response. This pulse was simulated by activating 25% of all cortical fiber terminals (we will call this a stimulus of strength 25%).
The simulated epidural response shows 3 large waves of activity separated by 1.3 ms, followed by a smaller fourth wave (Fig. 3), similar to responses recorded in humans. Simulated responses of single-cell membrane potentials, average layer membrane potentials, and membrane potential plots show that responses in all populations occurred in a very synchronous manner. Excitatory neurons in L2/3 and L5 spiked immediately after the stimulus, whereas those in L6 and T spiked at successively later times. The spiking in T was relatively small, which is to be expected, given that the thalamus is believed to play a minor role in the cortical response to TMS (Ziemann and Rothwell 2000). This firing pattern suggests that a TMS pulse elicits a response initially from layers 2, 3, and 5, which then propagates to layer 6 and then to the thalamus. The firing pattern also suggests that individual layer 5 cells spike multiple times, whereas the remainder of cortical excitatory cells spike only once. This is consistent with observations made in monkeys, where pyramidal tract cells spike multiple times at a frequency similar to that observed in epidural recordings (Edgley et al. 1997).
After the simulated TMS pulse, inhibitory neurons were generally more active than excitatory cells throughout the circuit and fired before, during, or after the excitatory cells in the population spiked. The role of excitatory and inhibitory cells is further examined below.
I-waves were recruited in a dose-dependent fashion
I-waves in motor cortex can be elicited in a dose-dependent fashion. We created a dose–response curve by running simulations across a wide range of stimulus strengths. The integrated epidural response was measured over the 10 ms after TMS as a means of judging total activation. We noted the stimulus strength at which the number of elicited I-waves changed (Fig. 4A). Four or more I-waves were elicited in the range of stimulus strengths examined (0–50%), with successive I-waves appearing with increasingly stronger stimuli (I-waves just began to emerge at 7, 14, 18, and 25% stimulus strength for I1–I4, respectively). Note how the dose–response curve demonstrates that each individual I-wave initially appeared small in size and grew with a stronger stimulus before the next I-wave was added to repeat this process.
TMS evokes a similar dose–response curve for individual and summed I-waves in human motor cortex. Epidural and EMG responses in humans have shown that as the strength of stimulation increases more I-waves are recruited and the strength of each individual I-wave increases (Nakamura et al. 1996).
To shed light on the differential neuronal responses after TMS, a dose–response curve was generated for excitatory and inhibitory cells in L2/3 and L5 by summing total spikes in the 10-ms period after a simulated TMS pulse (Fig. 4B). Note how a greater number of inhibitory than excitatory neurons spiked at low stimulus intensities, reflecting the fact that inhibitory neurons have a lower response threshold to TMS (an effect particularly prevalent in L2/3). This is to be expected from the properties of the model, including synaptic conductances and intrinsic properties that make inhibitory cells more excitable than excitatory cells.
Modifying cortical excitability produced changes in the I-waves similar to those observed in vivo
Modifying the excitability of cortex through pharmacological or behavioral means is a useful way to explore the relative role of excitatory and inhibitory elements. We simulated decreasing cortical excitability by adjusting GABAA levels. The control stimulus we used (35% strength) was capable of generating 4 I-waves (Fig. 5A). When we increased the peak conductance of GABAA by 25% the fourth I-wave was lost completely and the third I-wave was significantly reduced in size (Fig. 5B). Next, we increased cortical excitability by adjusting AMPA conductance. We simulated a TMS pulse (22.5% strength) capable of inducing 3 I-waves as our control (Fig. 5C). When we increased the peak conductance of AMPA by 15% a fourth I-wave appeared (Fig. 5D).
These results are comparable with those reported in the literature. When the GABAA agonist lorazepam is used to increase cortical inhibition, a stimulus that normally evokes 4 I-waves will evoke only 3 I-waves (Di Lazzaro et al. 2000). When subjects increase cortical excitability by forming a tightly clenched fist, a stimulus normally capable of evoking 3 I-waves will evoke an additional I-wave (Di Lazzaro et al. 1998b).
Applying TMS in paired-pulse paradigms produced patterns of inhibition or facilitation similar to those observed during in vivo experiments
Delivering subthreshold and/or suprathreshold TMS pulses in rapid succession in motor cortex elicits distinct patterns of response and can be used to explore excitation and inhibition. To implement various paired-pulse paradigms, we first determined, by comparison with Di Lazzaro et al. (2001b), that our model could produce a thresholdlike response at 7.5% stimulus strength. This was used as the threshold to study paired-pulse responses in a subthreshold/suprathreshold paradigm and a suprathreshold/subthreshold paradigm.
In humans, the subthreshold/suprathreshold paradigm has been examined by comparing EMG responses. These responses show inhibition at early interstimulus intervals (ISIs) up to about 7 ms and facilitation at later ISIs (Kujirai et al. 1993; Ziemann et al. 1996).
Although these studies measured only EMG responses, Hanajima et al. (1998) suggest that these experiments used suprathreshold pulses that evoke more than one I-wave. Therefore we simulated this paradigm by modeling a subthreshold pulse (6.5%) followed by a suprathreshold pulse (15%) normally capable of eliciting 2 I-waves. The simulation was run from the same initial conditions for ISIs varying between 1 and 12 ms. We measured the integrated epidural response for 7 ms after the TMS pulse. Compared with the response to the suprathreshold stimulus alone, this response was inhibited at early ISIs (1–8 ms) and facilitated at later ISIs (9–12 ms), as shown in Fig. 6A.
In humans, a variant of the subthreshold/suprathreshold paradigm has been used by Di Lazzaro et al. (1998c) to study individual I-waves at early ISIs with epidural recordings. They found that I1 is unaffected, I2 is somewhat inhibited, and I3 is strongly inhibited at the early time points in this paradigm. Two inhibitory mechanisms are believed to be at work in this response. At the earliest time intervals, the fibers excited by TMS are thought to be in a refractory state, whereas at later ISIs GABAergic inhibitory circuits come into play (Fisher et al. 2002; Hanajima et al. 2003).
To simulate comparable results, we modeled a subthreshold (5%) pulse followed by a suprathreshold pulse (20%) normally capable of eliciting 3 I-waves. The simulation was run from the same initial conditions for ISIs varying between 1 and 8 ms. We compared the peak size of each I-wave to the size of that I-wave in response to the suprathreshold stimulus alone (Fig. 6B). We found I1 was almost unaffected, I2 was moderately inhibited, and I3 was strongly inhibited. Furthermore, our results support the importance of a refractory effect of TMS in this paradigm: When the simulations were run without the refractory effect of TMS on the synapses, no inhibition was observed at an ISI of 1 ms.
In humans, the delivery of a suprathreshold pulse followed by a subthreshold pulse at short intervals has been explored by Ziemann et al. (1998). They found that this paradigm leads to facilitation in a wavelike pattern with peaks at ISIs of 1.2, 2.3–3.0, and 4.1–4.5 ms.
We modeled this paradigm using a suprathreshold pulse (22.5% of presynaptic fibers) capable of eliciting 3 I-waves followed by a subthreshold pulse (6.5% of presynaptic fibers). Comparing the paired-pulse response to control (Fig. 6B) over a range of ISIs from 1 to 5 ms shows that the suprathreshold response was enhanced compared with controls. The precise amplitude of this facilitation is not comparable with our simulated data because the human recordings were made from muscle responses, whereas ours were made directly from simulated L5 cortical activity. However, this enhancement does follow a wavelike pattern with peaks at ISIs of 1.0, 3.2, and 4.6 ms as was observed in humans.
Simulating activation of afferents from the somatosensory area produced a response similar to that obtained with anterior–posterior coil orientation in vivo
The sequence of I-wave induction and the paired-pulse curve change when TMS is used to induce currents in an anterior–posterior (AP) orientation, as opposed to the more common posterior–anterior (PA) orientation. Specifically, the I-wave that appears just above the motor threshold has a latency corresponding to that of I3 [3.6–4.5 ms after the direct wave (Di Lazzaro et al. 2001b; Hanajima et al. 1998), compared with 1.0–1.4 ms for I1 (Di Lazzaro et al. 1998a, 2001a; Nakamura et al. 1996)].
It has been hypothesized that these differences may be attributable to different populations of fibers being excited by AP versus PA currents (Di Lazzaro et al. 2001b). Although fiber terminals are the most likely site for activation by TMS, it is also possible that large afferent axons from premotor and somatosensory areas, which constitute the main cortical input to motor cortex (Braitenberg and Schüz 1991; DeFelipe et al. 1986; Sutor et al. 2000), may be especially sensitive to AP currents. These afferents bend into motor cortex 0.5 mm from layer 5 and 1.5 mm from layer 3 (Rockel et al. 1980), and it is known from modeling and peripheral nerve stimulation studies that bends in large fibers have a low threshold for TMS activation (Abdeen and Stuchly 1994; Maccabee et al. 1993).
To explore this possibility, we simulated the stimulation of PM or SI fibers at a point equivalent to where they bend into cortex (Fig. 7A), using a 1 m/s conduction velocity (Swadlow 1994) to calculate the appropriate delay from the bend to the fiber terminals. As shown in Fig. 7B, stimulation of PM fibers resulted in a synchronous response 2.1 ms after the TMS pulse. By contrast, stimulation of SI fibers resulted in a synchronous response 4 ms after the TMS pulse, which is within the time range of when I3 would be expected [3.6 to 4.5 ms (Di Lazzaro et al. 2001b; Hanajima et al. 1998)]. This longer latency is the result of SI afferents projecting exclusively to L2/3, requiring the TMS-induced perturbation to travel to L2/3 before reaching L5, whereas PM afferents project directly to L5 (as can be seen in Fig. 7A). For comparison, stimulation of fiber terminals within MP, which is how PA stimulation was simulated in other sections of the paper, resulted in synchronous firing 1.4 ms after TMS, as expected for I1. Altogether, these results suggest that the changes in I-wave latency observed after AP stimulation are consistent with preferential activation of somatosensory afferents.
The effects of TMS on cortical circuits: a modeling perspective
The simulations presented above demonstrate that, by taking into account the anatomical and physiological organization of the cortical motor system, the model was able to reproduce a variety of data from in vivo experiments. Having thus validated the model, it is now appropriate to use it to obtain a detailed view of cortical activity after one or 2 TMS pulses.
Effects of a single TMS pulse
To obtain a mechanistic understanding of how motor cortical circuits might respond to TMS, we simulated the delivery of one or 2 TMS pulses and examined membrane potentials, synaptic currents, and population firing rates in each cortical layer (Figs. 8 and 9). In the present analysis, synaptic currents were categorized as originating from synapses activated directly by TMS, intralayer excitatory input, interlayer excitatory input, or local inhibitory input. These synaptic currents were further separated according to the particular wave of spiking from which they originated. Excitatory currents constituted the combined response of AMPA and NMDA channels, whereas inhibitory currents constituted the response of GABAA channels. GABAB channels did not change considerably over the time course examined, which is expected because GABAB is believed to play a role primarily in TMS responses of longer latency than those considered here. For brevity, only currents generated by afferents to excitatory cells are shown. Excitatory and inhibitory cells in a given layer receive inputs from the same sources, so the timing of the synaptic input they receive is the same. Thalamic layers are not shown because they do not provide significant input to the cortex during the time period examined and, furthermore, because experimental evidence has demonstrated that the thalamus is not important to the generation of I-waves (Amassian et al. 1987; Di Lazzaro et al. 2003).
Figure 8 shows the model's response to a single simulated TMS pulse that activates 20% of fiber terminals in MP (a stimulus capable of eliciting 3 I-waves). Briefly, the course of events is as follows.
1) The immediate effect of the TMS pulse is the synchronous activation of 20% of all populations of cortical fiber terminals. These currents are strongly depolarizing because excitatory AMPA channels are more numerous and open more rapidly than inhibitory GABAA channels. This net depolarization is sufficient to cause spiking in excitatory and inhibitory neurons in L2/3 and L5 but causes spiking in inhibitory neurons only in L6. The spiking of the L5 excitatory cells is recorded as I1 in the epidural response.
2) Intralayer excitatory, interlayer excitatory (from L5), and local inhibitory inputs from the cells that spiked in step 1 activate synapses onto L2/3 cells. All currents, once activated, require many milliseconds to completely decay. Total current described here is the summation of all currents that are still active. At this point, the total current in L2/3 is depolarizing. However, most cells in L2/3 are in a refractory state from firing in the previous step, therefore minimal spiking occurs in L2/3. Meanwhile, intralayer excitatory, interlayer excitatory (from L2/3), and local inhibitory inputs activate synapses onto L5 cells, leading to a depolarizing total current. L5 cells receive less inhibitory input and have a shorter refractory period than that of L2/3 cells (see methods for details), and thus L5 excitatory cells spike about 1.3 ms after I1 (which is recorded as I2 in the epidural response) and inhibitory cells, because of their marginally longer refractory period than that of L5 excitatory cells, spike shortly thereafter. Meanwhile, interlayer excitatory (primarily from L5 and less from L2/3) and local inhibitory inputs activate synapses onto L6 cells, leading to a depolarizing total current. This net depolarization is sufficient to cause the spiking of cells in L6.
3) Interlayer excitatory input from L5 cells spiking in step 2 activates synapses onto neurons in L2/3. Because excitatory currents decay more rapidly than inhibitory currents, the net current from synapses previously activated becomes hyperpolarizing, and therefore the total current becomes nearly zero. As a result, L2/3 excitatory cells, although no longer refractory, do not spike again. However, inhibitory cells have shorter refractory periods and are more sensitive to input than excitatory cells, and thus L2/3 inhibitory cells spike. Meanwhile, intralayer excitatory, weak interlayer excitatory (from L6), and local inhibitory inputs activate synapses onto L5 cells, resulting in a depolarizing total current (inhibition is weaker in L5 than in L2/3; see methods for details). This net depolarization is sufficient to cause the synchronous spiking of L5 cells about 1.3 ms after I2 (which is recorded as I3 in the epidural response) and inhibitory cells shortly thereafter because of their slightly longer refractory period. Meanwhile, intralayer excitatory, interlayer excitatory (from L5), and local inhibitory inputs activate synapses onto L6 cells, resulting in a depolarizing total current. This net depolarization is insufficient to overcome the refractory period of L6 excitatory cells, but is sufficient to cause L6 inhibitory cells to spike.
4) Intralayer excitatory and local inhibitory input from cells spiking in step 3 activate synapses onto L5 cells, resulting in nearly zero total current. As a result, L5 cells do not spike again. Meanwhile, interlayer excitatory (from L5) and local inhibitory inputs activate synapses onto L5, resulting in nearly zero total current. As a result, L6 cells do not spike again.
Effects of a TMS paired pulse
Figure 9 shows the model's response to a simulated subthreshold TMS pulse that activated 5% of fiber terminals followed 4 ms later by a suprathreshold stimulus that activated 20% of fiber terminals. In the model, Fig. 4B shows that this paradigm resulted in inhibition of the I2 and I3 normally induced by the suprathreshold stimulus alone. The course of events can be described as follows.
1) The immediate effect of the subthreshold TMS pulse is the synchronous activation of a small percentage of all cortical presynaptic fiber terminal populations. These currents are mildly depolarizing. This net depolarization is sufficient to cause the spiking of inhibitory neurons in all layers and a small number of excitatory cells in L5. The majority of excitatory cells, which are less sensitive to input than inhibitory cells, do not spike.
2) The immediate effect of the suprathreshold TMS pulse is the synchronous activation of a large percentage of all cortical presynaptic fiber terminal populations. This TMS pulse as well as interlayer excitatory (from L5) and local inhibitory inputs from cells spiking in step 1 activate synapses onto L2/3 cells. All currents, once activated, require many milliseconds to completely decay. Total current described here is the summation of all currents that are still active. At this point, total current in L2/3 is mildly depolarizing. This net depolarization is sufficient to induce a small number of excitatory and most inhibitory cells in L2/3 to spike. The suprathreshold TMS pulse as well as intralayer excitatory and local inhibitory inputs activate synapses onto L5 cells, resulting in a depolarizing total current. This net depolarization is sufficient to cause the spiking of L5 cells (which is recorded as I1 in the epidural response). The suprathreshold TMS pulse as well as interlayer excitatory (from L5) and local inhibitory inputs activate synapses onto L6 cells, resulting in a mildly depolarizing total current. This net depolarization is sufficient to cause the spiking of L6 inhibitory cells.
3) Intralayer excitatory, interlayer excitatory (from L5), and local inhibitory inputs from cells spiking in step 2 activate synapses onto L2/3 cells, resulting in nearly zero total current. As a result, no further spiking occurs in L2/3. Intralayer excitatory, weak interlayer excitatory (from L2/3), and local inhibitory inputs activate synapses onto L5 cells, resulting in nearly zero total current. This current is sufficient to induce a small amount of spiking in L5 (which is recorded as a weak I2 in the epidural response). Intralayer excitatory, interlayer excitatory (from L5), and local inhibitory inputs activate synapses onto L6 cells, resulting in a depolarizing total current. This net depolarization is sufficient to cause the spiking of L6 cells.
4) Intralayer excitatory and local inhibitory inputs activate synapses onto L5 cells, resulting in nearly zero total current. As a result, L5 cells do not spike again. Intralayer excitatory and local inhibitory inputs activate synapses onto L6 cells, resulting in a slightly depolarizing total current. This net depolarization is sufficient to cause the spiking of L6 inhibitory cells.
In all layers, most currents after the suprathreshold pulse were inhibited by the subthreshold pulse (Fig. 9, right). Currents to L5 from interlayer sources (predominantly L2/3) were particularly inhibited.
It should be noted that this response description is applicable only for short ISIs. At longer ISIs (>8 ms), the synaptic channels opened by the subthreshold pulse have time to close and do not directly contribute to the circuit's response to a suprathreshold pulse. However, over this time the currents do affect the membrane potential of neurons throughout cortex, shifting membrane potentials so that neurons are more likely to fire, which results in the facilitation observed at later I-waves (data not shown).
In this paper, we investigated the effects of TMS on cortical circuits by constructing a large-scale model of a portion of the thalamocortical system and simulating the delivery of TMS pulses. To date, simulations of TMS have been limited to exploring the distribution of induced electric fields in cortical structures (Nadeem et al. 2003; Ravazzani et al. 1996) and their effects on single cells (Nagarajan et al. 1993). Therefore the model presented here is the first to simulate the effects of TMS by accounting for the complexities of cortical circuits. Simulations using this model show that large-scale modeling can complement the experimental investigation of the cortical response to TMS and propose a comprehensive picture of TMS-induced activity. Specifically, using a consistent set of parameters, the model is able to 1) reproduce spontaneous activity patterns similar to those observed in vivo; 2) respond to a simulated TMS pulse with I-waves of the appropriate timing; 3) generate a realistic dose–response curve; 4) replicate experiments performed under conditions of modified cortical excitability; 5) produce realistic paired-pulse results; and 6) mimic the different responses observed during AP stimulation.
Having been validated with respect to a large number of experimental results, we took advantage of the model to obtain a detailed picture of the effects triggered by the delivery of a single TMS pulse to a cortical area (Fig. 8). In short, the emerging picture is as follows: First, a suprathreshold TMS pulse directly activates a large fraction of cortical fiber terminals, thereby inducing spiking activity in both excitatory and inhibitory neurons in all cortical layers. These spiking volleys lead to excitatory and inhibitory currents, which are especially prominent at the comparatively strong synapses made by neurons from L2/3 and L5. Net depolarizing currents are strong enough to induce further firing in excitatory neurons in L5 because of their short refractory period and high ratio of excitatory to inhibitory synapses, but not in other layers. L5 neurons respond to this net depolarization with one to 3 more spikes, which are timed by their intrinsic neuronal properties. Finally, as a significant number of synaptic channels opened by earlier inputs close, synaptic currents become insufficient to induce further firing in L5 neurons, and the local effects of the TMS pulse fade.
By carefully evaluating the available experimental evidence, Ziemann and Rothwell (2000) recently considered as many as 5 different scenarios of how a TMS pulse may lead to the generation of I-waves, 2 of which are especially relevant in light of our simulated results. According to one scenario, TMS induces a large, primary input to cortical cells that causes them to spike with a timing dictated by their intrinsic properties. The present simulations confirm a role for intrinsic cellular properties in the timing of I-waves. However, they also indicate that secondary firing within cortical circuits, triggered indirectly by TMS, is required for the generation of later I-waves. According to another scenario, TMS activates a chain of excitatory and inhibitory neurons that provide alternating excitatory and inhibitory input to layer 5 cells. The simulations confirm the contribution of secondary waves of activation. However, they indicate that the secondary activation evolves as a single depolarizing envelope rather than as a series of depolarizing and hyperpolarizing waves. Thus the actual course of events underlying the generation of I-waves is best approximated as a combination of intrinsic neuronal properties and circuit interactions.
The effects of TMS on cortical circuits have also been investigated by recording muscle and epidural responses to paired TMS pulses to motor cortex. A number of paired-pulse paradigms exist that evoke responses that are inhibited or facilitated compared with a single pulse, depending on the strength of each pulse and the ISI used. In this paper, we focused on the reduced response to a suprathreshold TMS pulse when it is delivered 4 ms after a subthreshold pulse. Based on the results presented here, the response of cortical circuits to this paired-pulse paradigm is predicted to unfold as follows: First, the subthreshold TMS pulse directly activates a small fraction of cortical fiber terminals in all cortical layers. This small direct activation induces spiking activity primarily in inhibitory neurons as a result of their greater sensitivity to depolarizing inputs. As a consequence, hyperpolarizing synaptic currents are induced throughout all cortical layers. When the suprathreshold TMS pulse is delivered at a short ISI (4 ms), these hyperpolarizing currents are sufficient to prevent the firing of L2/3 excitatory cells that normally would spike in response to the suprathreshold TMS pulse. On the other hand, L5 neurons do fire, thanks to a comparatively high ratio of excitatory to inhibitory synapses. However, without subsequent input from L2/3, most L5 neurons do not spike more than once. At a longer ISI (>9 ms), the inhibitory effect induced by the subthreshold pulse fades as GABAA synaptic channels close.
This predicted scenario is in line with the observations made using epidural recordings in humans, which show that a subthreshold pulse causes the inhibition of I2 and I3, but not of I1 (Di Lazzaro et al. 2000, 2004; Hanajima et al. 1998). The present simulations confirm that, as was originally hypothesized, the subthreshold pulse inhibits later I-waves by activating local inhibitory neurons. Moreover, the model demonstrates that the main effect of the activation of inhibitory neurons is to prevent the firing of L2/3 neurons. As a consequence, L2/3 neurons do not contribute excitatory input to L5 neurons and later I-waves cannot occur.
In summary, by incorporating the detailed anatomy and physiology of motor cortex, the present model confirms and extends experimental observations and predicts in considerable detail the sequence of cortical events that are triggered by one or more TMS pulses. Other predictions of the model include a reentrant thalamocortical response to TMS stimuli, and the preferential activation of fibers terminating in superficial layers by AP stimulation. This preferential fiber activation does not suggest that a TMS pulse with one orientation exclusively activates fiber tracts, whereas a TMS pulse with the opposite orientation exclusively activates fiber terminals. Rather, we speculate that TMS, at certain orientations and intensities, may activate large fiber tracts in addition to fiber terminals. We further speculate that an AP-oriented TMS pulse activates somatosensory fibers at a lower intensity than a similarly oriented pulse activates fiber terminals. Therefore a TMS pulse near threshold activates primarily motor cortical afferents from somatosensory cortex, resulting in an I-wave with the timing of I3, whereas a stronger TMS pulse activates fiber terminals as well as fiber tracts, resulting in the generation of early and late I-waves.
As with any model (cf. Traub et al. 2004), these predictions should be qualified because of inevitable simplifications and limitations. For example, although the model was created by accounting for a large number of details about the anatomy and physiology of the motor cortical system, relevant parameters were not always available or were loosely constrained. For these situations, we performed extensive parameter searches to ensure that the chosen values resulted in a model that generated spontaneous activity and the responses to TMS stimulation observed in vivo. Furthermore these parameter searches have made it clear that the model produces the expected results for a wide range of possible parameter values. In choosing a level of modeling that was sufficient for understanding the effects of TMS on cortical circuits yet could be computationally feasible, we made several design choices. We did not explicitly model magnetic fields, although we took into account experimental and theoretical literature on how TMS affects neural fibers. We modeled a layered motor cortex containing explicitly modeled minicolumns of excitatory and inhibitory neurons (>30,000) endowed with AMPA, NMDA, and GABA receptors. Moreover, these neurons were linked by a network of intra- and interlayer connections, as well as by corticothalamocortical and interareal connections (>5,000,000). Although the scale of the model precluded computationally intensive features such as multiple-compartment neurons with explicit dendritic and axonal arborization, the model was nevertheless able to replicate a wide array of experimental results, suggesting that it captured most of the relevant features. We restricted the model to a limited section of motor cortex, corresponding to the area of peak TMS intensity. However, I-waves and paired-pulse responses occur so rapidly that their generation should involve primarily local circuits. We also did not attempt to model the generation of actual movements or EMG-like responses by TMS. However, we did model epidural-like responses, which have been measured in several human studies. Finally, although we were able to account for the effects of single- and paired-pulse paradigms, we did not address certain other effects of TMS on cortical circuits. Nevertheless, the current model presents a useful starting point for exploring the cumulative effects of fast and slow TMS trains, the induction of plastic changes in cortical circuits, and the long-range effects of TMS-induced volleys on connected cortical regions.
TMS has a growing number of clinical and research applications with the potential to investigate and treat disorders such as depression, schizophrenia, or stroke. For these reasons, it is important that the mechanisms of action of TMS on cortical circuits are thoroughly understood. As shown here, large-scale modeling can overcome experimental difficulties to provide a comprehensive picture of the synaptic, neural, and circuit responses triggered by TMS pulses.
This work was supported in part by National Alliance for Research on Schizophrenia and Depression as well as National Institute of General Medical Services Grant GM-07507.
We thank M. Massimini, F. Ferrarelli, and R. Huber for helpful conversations.
The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “advertisement” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
- Copyright © 2005 by the American Physiological Society