| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
1Neuroscience Training Program and 2Department of Psychiatry, University of Wisconsin, Madison, Wisconsin
Submitted 2 December 2004; accepted in final form 17 March 2005
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
; 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.
|
METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
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 x 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.
|
; 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 x 5 group of topographic elements from each layer (for a total of 225 cells per column). MP consisted of an 8 x 8 array of these columns interspersed with an 8 x 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.
INPUT AREAS.
MP received input from PM, SI, and T. PM and SI were modeled as 40 x 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).
THALAMUS.
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 x 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 x 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.
Connectivity
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.
|
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.
INTRATHALAMIC CONNECTIONS.
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.
THALAMOCORTICAL CONNECTIONS.
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 CONNECTIONS.
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.
Transmission delays
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.
Model neurons
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)]
 |
 |
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.
|
Synaptic input, Isyn, constituted the sum of all synaptic channel currents
 |
 |
1...m, respectively, and the time to peak tpeak as
 |
; Otis et al. 1993; Stern et al. 1992).
|
; 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
 |
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).
Simulation techniques
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.
|
RESULTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
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).
|
), 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).
|
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.
|
). 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).
|
). 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.
|
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 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).
|
|
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 spi
