| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
Department of Psychiatry, University of Wisconsin—Madison, Madison, Wisconsin
Submitted 1 September 2004; accepted in final form 26 October 2004
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
). The slow oscillation is composed of a hyperpolarized phase or "down-state" during which neurons are deeply hyperpolarized and remain silent for a few hundred milliseconds. The down-state is followed by a depolarized phase or "up-state," also lasting several hundred milliseconds, during which neurons fire at rates that are even higher than in quiet wakefulness. These fluctuations in membrane potential are reflected in the cortical EEG as slow waves of high-voltage activity.
The slow oscillation is the fundamental cellular phenomenon that groups and organizes sleep rhythms such as slow-wave activity and sleep spindles ( Steriade 2003). After its discovery in anesthetized cats ( Steriade et al. 1993), the slow oscillation has been investigated during natural sleep in vivo ( Achermann and Borbely 1997; Steriade et al. 2001), in cortical slabs ( Timofeev et al. 2000), in vitro in cortical slice preparations ( Mao et al. 2001; Sanchez-Vives and McCormick 2000), and in computo ( Bazhenov et al. 2002; Compte et al. 2003). These studies have revealed that both intrinsic currents and various kinds of synaptic interactions are involved in initiating, maintaining, and terminating the slow oscillation, and that corticocortical circuits are involved in synchronizing it ( Amzica and Steriade 1995b).
To integrate the information gathered from these different experimental approaches, we have constructed a large-scale model of the thalamocortical system that aims to provide a coherent account of the transition from wakefulness to sleep and the generation of the slow oscillation at several different levels—from ion channel kinetics to global EEG phenomena. The model incorporates key aspects of the neuroanatomical organization of the thalamocortical system, including two visual cortical areas subdivided into multiple layers, corresponding thalamic and reticular sectors, and several millions of intra- and interareal connections linking >65,000 spiking neurons. Moreover, the model incorporates several types of intrinsic conductances (mediating the hyperpolarization-activated cation current Ih, low-threshold calcium current IT, persistent sodium current INa(p), potassium leak current IKL, depolarization-dependent potassium current IDK—representing Ca2+ and Na+-dependent K+ currents) and synaptic currents [
-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA), N-methyl-D-aspartate (NMDA), 
-aminobutyric acid-A (GABAA), -aminobutyric acid-B (GABAB)].
Because of these properties, the model is able to reproduce experimental data ranging from intracellular traces and multiunit rasters to optical imaging-like voltage patterns and EEG-like field potentials. Moreover, by simulating changes in intrinsic currents arising from the reduced release of neuromodulators, the model can switch from a waking to a sleep mode of activity. Specifically, in the waking mode the model reproduces spontaneous activity patterns as well as selective responses to visual stimuli that are seen in vivo. After transitioning to the sleep mode, the model engages in slow oscillations that closely resemble those observed experimentally. By providing a comprehensive view of all system variables and by permitting idealized "experimental" manipulations, the model provides a self-consistent account of the mechanisms responsible for the initiation, maintenance, and termination of the slow oscillation and of its synchronization within and across thalamocortical circuits.
|
METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
Regional organization
PRIMARY CORTICAL AREA. The model (Fig. 1) is organized in regions and pathways consisting of a primary and a secondary area of visual cortex, two corresponding regions of the dorsal thalamus, and two regions of the reticular thalamic nucleus. The primary visual area (Vp) represents a restricted portion of cat striate cortex (area 17) and it contains units with small receptive fields that are selective for oriented segments. The simulated cortex is divided into 3 layers with different patterns of afferent, efferent, and local connectivity corresponding to supragranular layers (L2–3), infragranular layers (L5–6), and layer 4 (L4).
|
; Kisvarday et al. 1997; Lowel et al. 1987) and in terms of the average separation of neighboring axonal termination patches in supragranular layers ( Kisvarday and Eysel 1992). To constrain our model, we assume that an area of roughly 1.0 mm2 forms a basic macrounit in the striate cortical mosaic array. Vp is scaled to span 64 (8 x 8) such macrounits. Therefore Vp corresponds to approximately 0.8 cm2 of striate cortical surface and spans a monocular patch of 8 x 8° in the parafoveal visual field. Each macrounit in the model contains discrete groups of orientation selective cells. In the present model, 2 groups of cells selective for vertical and horizontal oriented input represent a simplified version of the X pathway for one eye.
Figure 2 shows the detailed orientation-selective, feedforward, and feedback circuitry for one horizontally selective and one vertically selective macrounit. Each topographic location (topographic element) in the model cortex is considered to correspond to a cortical column, which is represented by 9 model neurons (2 excitatory and 1 inhibitory for each of the 3 layers). All topographic elements in Vp are organized in maps of 40 x 40 elements for each of the 2 modeled orientation selectivities (horizontal and vertical). Orientation selectivity is achieved by the convergence of afferents from an oriented rectangular region in Tp onto individual cortical cells in L4 and L5–6. The subdivision of the modeled cortical areas in elements spanning all layers reflects the developmental, anatomical, and physiological evidence for a basic columnar organization of neocortex ( Gilbert 1993; Mountcastle 1997, 1957; Rakic 1995).
|
), each topographic element corresponds to a cortical column with a surface area of 1,454 µm2 and containing approximately 94 neurons. Because we explicitly model only 9 cells for each topographic element, each modeled cell represents the activity of approximately 10 cortical neurons, making the topographical elements in the model comparable to the basic cortical modules proposed by Peters and Payne (1993). Although the ratio of modeled excitatory/inhibitory cells (66%/33%) is not exactly the same as that observed in vivo (about 80%/20%) because of computational considerations, the observed ratio of excitatory/inhibitory synapses is maintained (see following text).
SECONDARY CORTICAL AREA.
The secondary visual area (Vs) corresponds to an extrastriate area located along the ventral occipitotemporal pathway. Although Vs does not represent in detail any particular region of visual cortex, we use area 21 in the cat as a reference, which is the presumed homolog of cortical area V4 in the monkey ( Payne 1993). Vs is assumed to be about half the size of Vp [in the monkey, V1 is 1,120 mm2 and V4 is 540 mm2 ( Felleman and Van Essen 1991)]. In the model, Versus is based on some general properties associated with extrastriate areas (e.g., an enlargement of receptive fields) and with termination patterns of "forward" and "backward" corticocortical projections ( Felleman and Van Essen 1991; Van Essen et al. 1992). Vs contains neurons that are selective for either vertical lines, horizontal lines, or line crossings, organized in a coarse topographic map. For each of its 3 selectivities, Vs has a map of 30 x 30 elements (for a total of 24,300 model neurons) as compared with the 40 x 40 (totaling 28,800 model neurons) elements in Vp.
THALAMIC SECTORS.
According to Peters and Payne (1993), there is a rough correspondence between the number of X-cells in the lateral geniculate nucleus (LGN) and the number of basic cortical modules in area 17. We therefore model a geniculate map (Tp) composed of the same number of elements (40 x 40) as Vp. Each element of Tp contains 2 modeled neurons that correspond respectively to an X-relay cell and to an inhibitory interneuron. For simplicity of implementation, only the ON-portion of thalamic receptive fields is modeled. The secondary thalamic map (Ts) has 30 x 30 elements and its visuotopic arrangement has a much lower spatial resolution than that of Tp. Two sectors of the reticular nucleus, primary perigeniculate (Rp) and secondary higher-order (Rs), are modeled respectively as a 40 x 40 and a 30 x 30 map of inhibitory neurons.
Connectivity
In constructing the model, special emphasis was placed on the incorporation of realistic network properties, such as the spread and relative proportions of the various sets of connections composing the intra- and interregional thalamocortical circuitry. Specific patterns of arborization are classified as either focused or diffuse, on the basis of anatomical data. The focused connection pattern diverges for single arbors over a topographically registered region with a diameter of 5 target elements. Diffuse projections typically cover an area with a diameter of 25 elements for a single arbor. Contacts from individual arbors in the target area are made probabilistically according to Gaussian spatial density profiles. The proportion of synapses from different sources was used as a constraint in the parameterization of the various density profiles (Table 1). Two books, by Sherman and Guillery (2001) and White and Keller (1989), were particularly helpful in the development of this model.
|
; Mountcastle 1997): from layer IV to supragranular layers, from supragranular to infragranular layers, and from infragranular back to layer IV and to supragranular layers ( Callaway and Wiser 1996; Wiser and Callaway 1996, 1997). All these projections are made in a focused manner in the model. As a further constraint, we consider the proportion of synapses from different sources in each layer. As an example, for each simulated excitatory cell of layer IV, there are on average 40 interlaminar connections from infragranular layer and 23 intralaminar connections from layer IV, in close agreement with the 45% ratio reported in the cat striate cortex ( Ahmed et al. 1994). A similar connectivity was established in each map of Versus.
HORIZONTAL INTRALAMINAR CONNECTIONS.
Individual excitatory neurons in the supragranular layers of striate cortex have intralaminar horizontal projections that tend to be organized in patches of 200–400 mm in diameter. These patches typically interconnect neurons that have similar orientation preference ( Kisvarday et al. 1997). Patches originating from a single location extend over a region of roughly 2–4 mm ( Gilbert 1993). In the model, horizontal connections in the supragranular layer of Vp are made diffusely onto isoorientation cells, with an equivalent spread of 5.5 x 5.5 mm2. Intrinsic connections in the infragranular layer of Vp have a similar organization. Intralaminar connections in layer IV extend over a more limited area with a diameter of 15 elements. This reduced projective field reflects the more compact arborization in layer IV ( Douglas and Martin 2003).
INTRACORTICAL INHIBITORY CONNECTIONS.
The cerebral cortex contains many different types of GABAergic inhibitory interneurons ( Douglas and Martin 2003; Jones 1993). Among these, basket cells are ubiquitous and project mostly to the same layer where their soma is located. Double-bouquet cells are concentrated in supragranular layers ( Conde et al. 1994; Kawaguchi 1995; Kawaguchi and Kubota 1997; Peters and Sethares 1997) and their projections are organized in a restricted columnar arrangement that extends to most layers. There are indications that basket cells and other inhibitory interneurons act through fast GABAA-receptors, whereas double-bouquet cells may preferentially activate GABAB receptors ( Kang et al. 1994). In the model, basketlike cells provide a fast (GABAA-like) inhibition within each cortical layer to all cell types; double-bouquet analogs located in supragranular layers provide a slow (GABAB-like) inhibitory control of a narrow cylinder extended to all 3 layers.
The relationship between inhibition and orientation selectivity in the visual cortex is complex. However, some recent studies suggest that a single basket cell in the cat visual cortex provides input to surrounding regions representing the whole range of orientations, including iso- and cross-orientations to that basket cell soma ( Kisvarday and Eysel 1993; Kisvarday et al. 1994). In the model, we assume that lateral inhibition (GABAA) is provided equally to both of the modeled orientation selectivities in Vp. In Vs, half of the terminals of individual basketlike cells in Vs provide input to cells with similar selectivity to that of the parent soma, whereas the remaining half are split evenly between cells of other selectivities. The density profile of inhibitory connections was adjusted such that the relative proportions of inhibitory connections per layer are comparable to the values reported in the literature (i.e., about 10–20% of all synapses) ( Beaulieu and Colonnier 1985; Beaulieu et al. 1992).
FORWARD AND BACKWARD INTERAREAL CONNECTIONS.
According to several studies, backward connections are considerably more divergent than forward connections. This has been documented for projections from area MT to V1 and V2 of primates ( Krubitzer and Kaas 1989; Rockland and Knutson 2000; Shipp and Zeki 1989; Zeki and Shipp 1989) and from V2 to V1 ( Henry et al. 1991; Rockland and Van Hoesen 1994; Rockland and Virga 1989). Reconstructions of single axons indicate that forward projections from V1 and V2 ( Rockland 1992; Rockland and Knutson 2000; Rockland and Virga 1989) to V4 have discrete terminal clusters (2–4 clusters per axon, 250 mm wide), which are distributed over 2 to 2.5 mm. Conversely, individual axons from V4 to V1 diverge 
5 mm ( Rockland et al. 1994). These values should be compared with values around 2 to 5 mm for horizontal connections in V1. According to a classic description, forward projections tend to originate in superficial layers and to terminate in layer 4, whereas backward connections tend to originate from infragranular as well as supragranular layers and to terminate outside layer 4 ( Rockland and Pandya 1979). This basic scheme has since become considerably more complicated ( Felleman and Van Essen 1991; Maunsell and van Essen 1983). However, in the model, forward connections originating from the supragranular layer of Vp defined the 3 feature-specific responses in Vs. These selectivities resulted from biased convergent projections onto individual L4 neurons of Versus from either a 19 x 3 region of vertical selective cells, a 3 x 19 region of horizontal selective cells, or from both selectivities of Vp (9 x 3 and 3 x 9 regions, respectively). Backward projections from vertical and horizontal selective cells of Vs originate in the infragranular layer and terminate diffusely in the supragranular layer of Vp, targeting cells of similar orientation specificity. In contrast, backward projections from cross-selective cells extend to both selectivities in Vp. The laminar specificity of projections between Vp and Vs was consistent with that observed between cat areas 17 and 21 ( Rosenquist 1985).
THALAMIC CONNECTIONS.
Relay cells in the LGN have strong, driving connections to the cortex, and form collaterals only with the reticular nucleus (RT), whereas interneurons in the LGN inhibit other interneurons as well as relay cells with a focused connectivity pattern. RT neurons make diffuse connections within the RT nucleus and to thalamic relay nuclei ( Dubin and Cleland 1977). In the model, local interneurons produce (fast) GABAA-mediated inhibitory postsynaptic potentials (IPSPs) in thalamocortical relay cells. Thalamic to reticular projections (i.e., from Tp to Rp and from Ts to Rs) are made according to the focused connection scheme. RT projections target their corresponding relay sectors of the thalamus in a diffuse manner. Both GABAA and GABAB IPSPs mediate the RT inhibition of thalamic relay and interneurons, in accord with the inhibitory effect observed when RT cells fire tonically ( Kim and McCormick 1998; Kim et al. 1997; Pinault and Deschenes 1992).
THALAMOCORTICAL AND CORTICOTHALAMIC CONNECTIONS.
X-cells in cat laminae A and A1 of the LGN send axons that terminate mainly in layer IV and VI of area 17 ( Freund et al. 1989; LeVay and Gilbert 1976; Leventhal 1979). In the model, each simulated cell in L4 of Vp received connections selected from an 8 x 2 region of the thalamic map for the vertical selectivity (2 x 8 for the horizontal selectivity). Infragranular cells received about half as many connections from the same geniculate regions. The convergence of projections from these horizontal or vertical patches within the thalamus promoted orientation-specific responses in the cortex. Note that, in the model, the same X-cell targets both horizontal and vertical cortical cells, such that its arbor extends over at least half of an orientation cycle or 0.55 mm. This dimension is consistent with anatomical evidence that X axonal terminals form a single elongated clump in area 17, about 1 mm long x 0.6–0.8 mm wide ( Freund et al. 1985). Versus cells in L4 and the infragranular layer receive thalamocortical projections converging from a region of Ts with a diameter of 4 elements. Thalamocortical projections account for about 8% of all connections received by layer IV neurons, consistent with anatomical estimates ( Ahmed et al. 1994; Latawiec et al. 2000; Peters and Payne 1993). Corticothalamic axons descend from the infragranular excitatory cells into their corresponding thalamic relay sectors, contacting all cell types present in these structures, consistent with anatomical data ( Montero 1991; Robson 1983; Weber et al. 1983). En route, such fibers send collaterals to the RT nucleus. Consistent with experimental observations ( Golshani et al. 2001), corticoreticular projections are substantially stronger (2.5x) than corticothalamic projections. The topography of corticothalamic connectivity matches that of the thalamocortical connectivity ( Jones 2002).
Transmission delays
Transmission of signals within and across cortical areas occurs through several successive stages, including axonal conduction, synaptic delays, and postsynaptic potential (PSP) generation. Each of these stages is associated with delays in the transmission of a signal. Measured latencies between the firing of successive visual cortical areas in the cat have been estimated to lie between 5 and 15 ms ( Dinse and Kruger 1994) in the forward direction. Geniculocortical latencies may be even shorter ( Bullier and Henry 1979). Backward connections may be slower conducting. For instance, latencies from areas 18 and 19 to area 17 are 6 and 10 ms, respectively ( Bullier et al. 1988).
Because of the fact that the simulated cortices contain only 3 layers (instead of 6), we account for experimentally measured latencies along polysynaptic pathways by assuming comparatively longer transmission delays along certain pathways. Transmission delays for individual connections are sampled from Gaussian distributions with a SD of 1 ms. Each set of connections in the model is associated with a specific mean delay. Mean conduction delays are set to 2 ms for intralaminar connections and for most interlaminar connections. Infragranular to layer 4 connections are delayed on average by 7 ms, to account for disynaptic transmission through layers 5 and 6. Thalamocortical connections and forward connections from Vp to Vs have a mean delay of 3 ms, whereas corticothalamic connections and backward connections from Vs to Vp have a mean delay of 8 ms, again taking into account a disynaptic pathway through layers 5 and 6.
Model neurons
Both excitatory and inhibitory neurons are modeled as single-compartment spiking neurons incorporating Hodgkin-Huxley style currents. To model the contributions of key intrinsic currents, while preserving the computational efficiency of integrate-and-fire neurons that is necessary when computing a large-scale network, we devised a simplification of the fast-spiking currents (INa and IK). Model neurons thus behave like a hybrid between traditional integrate-and-fire neurons and full-fledged Hodgkin-Huxley neurons.
A dynamic threshold (
) is defined for each cell that determines at which membrane potential the cell should fire
 |
eq) determines the equilibrium threshold potential for both excitatory and inhibitory neurons. The threshold time constant 
determines the time to return to the equilibrium threshold. The specific values were chosen to match absolute refractory periods for different neuron types (Table 2).
|
 |
When the membrane potential V exceeds the threshold , a spike is generated by setting both V and instantaneously to the sodium reversal potential (ENa = 30 mV), modeling the contribution of the fast-spiking INa current. The activation of a fast potassium current during a spike is represented by a brief pulse (duration tspike, Table 2) with an amplitude of gspike = 1, thereby driving the membrane potential toward the potassium reversal potential (EK = –90 mV), while continuing to integrate intrinsic and synaptic currents. The integration of the fast hyperpolarizing current occurs faster than the membrane potential and is therefore governed by a "spike" time constant (spike m). The cell is unable to fire until V. With these three parameters, , spike, and tspike, we model key characteristics of spike generation including action potential width, afterhyperpolarization, and relative refractory period (Table 2).
The membrane time constants m are consistent with experimental data ( Baranyi et al. 1993; Connors et al. 1982; Kim and Connors 1993; Mason et al. 1991).
Two main categories of input currents contribute to the membrane potential, synaptic input (Isyn) and intrinsic currents (Iint), which are described below.
Synaptic channels
The synaptic input Isyn is the sum of all synaptic channel currents, 
Simulated synaptic channels provide voltage-dependent (NMDA-like) and voltage-independent (AMPA-like) excitation, as well as fast (GABAA-like) and slow (GABAB-like) inhibition. The conductance for each afferent i, on each channel j, specifies the amplitude and time course of the PSPs. The reversal potential for each channel Ej determines whether a current is inhibitory or excitatory. Electrical couplings between cortical inhibitory populations have been observed experimentally ( Galarreta and Hestrin 1999) but are not modeled here.
Synaptic activation is expressed as the change of a channel conductance, g(t), according to a dual-exponential response to single spike events, given by
 |
1 and 2 are the parameterizing the rise and decay time constants, respectively, and tpeak is the time to peak
 |
; Otis et al. 1993; Stern et al. 1992). The peak conductances gpeak were chosen to conform to a few simple constraints that led to regular network behavior. These constraints consisted of having: 1) peak excitatory postsynaptic potentials (EPSPs) of 1 mV in AMPA-like channels; 2) matched integrated excitatory postsynaptic currents (EPSCs) through AMPA and unblocked NMDA channels; 3) matched integrated inhibitory postsynaptic currents (IPSCs) through GABAA and GABAB channels. The activation of NMDA-like channels was expressed as
 |
). Based on this work, we model a 2-stage unblocking process with one component that unblocks quickly (about 1 ms) and a second slow component that unblocks in about 20 ms. Specific parameter settings for the different type of synaptic channels are listed in Table 3.
|
; Rosier et al. 1993); in addition, the effectiveness of these connections seems in part contingent on a concomitant depolarization of target cells that are not in visuotopic register with the sources of afferentation ( Bullier et al. 1988; Hirsch and Gilbert 1991; Salin and Bullier 1995).
Inhibition in the thalamus was mediated by fast (GABAA-like) synapses. Because of nucleus specific differences in the chloride reversal potential, reticulothalamic GABAA (TC) channels had a reversal potential more negative for thalamic relay cells (about –80 mV) than for cells found in the reticular nucleus (about –70 mV) ( Ulrich and Huguenard 1997).
Synaptic depression
There is substantial evidence that the rapid plasticity of excitatory and inhibitory synaptic responses is dominated by short-term depression and caused by the depletion of presynaptic pools of readily releasable neurotransmitter vesicles ( Zucker and Regehr 2002). In the model, short-term depression of both excitatory and inhibitory connections was based on a simple vesicle pool model ( Abbott et al. 1997; Galarreta and Hestrin 1998; Tsodyks and Markram 1997). Synaptic depression was modeled by scaling the peak conductance of a given synaptic channel by the size of the corresponding presynaptic pool of synaptic "vesicles." The dynamics of this pool was governed by the simple first-order equation dP/dt = –spike · 
P · P + (Ppeak – P)/P. The pool P decreases by the fraction P for each spike = 1. The pool recovers its peak value Ppeak according to the time constant P.
Intrinsic ion channel properties of thalamic and cortical neurons
Ion channel currents that influence intrinsic firing properties of thalamic and cortical neurons were modeled according to the Hodgkin-Huxley formalism Iint = gpeakmNh(V – Eint), where gpeak is the maximal conductance for the channel, m and h determine the activation and inactivation respectively (see following text), and Eint is the reversal potential for the given channel. The factor N allows the activation to occur on a different order than inactivation. The gating of activation and inactivation follows the same first-order kinetics equation: dx/dt = [x
(V) – x]/x(V) where x is the steady-state activation/inactivation value for the channel.
PACEMAKER CURRENT IH.
Ih is a noninactivating hyperpolarization-activated cation current that is believed to underlie a depolarizing "pacemaker" potential observed in many cells throughout the brain, including the thalamus and the cortex ( Huguenard and McCormick 1992; McCormick and Bal 1997; Robinson and Siegelbaum 2003). The activation variable mh for Ih is modeled by mh = 1/{1 + exp[(V – Vthreshold)/5.5]}, with Vthreshold = –75.0. The rate m of activation and deactivation also follows Huguenard and McCormick (1992): m = 1/[exp(–14.59 – 0.086V) + exp(–1.87 + 0.0701V)]. Only thalamic (Tp and Ts) and intrinsically bursting cells (IB: 30% of excitatory cells in L5–6) were endowed with Ih channels. The highest density of Ih channels in cortex is expressed in the dendrites of layer V neurons and was therefore included in L5–6 neurons ( Robinson and Siegelbaum 2003).
LOW-THRESHOLD CALCIUM CURRENT IT.
IT is a low-threshold fast-activating calcium current that underlies the generation of bursts in the thalamus and reticular nucleus ( Huguenard and Prince 1992; McCormick and Bal 1997). We use the formulation of IT from previous modeling work ( Destexhe et al. 1996a; Huguenard and McCormick 1992). Using the steady-state activation formulae, the activation variable is given by m = 1/{1 + exp[–(V + 59.0)/6.2]}, with the voltage-dependent time constant m = {0.22/exp[–(V + 132.0)/16.7]} + exp[(V + 16.8)/18.2] + 0.13. Inactivation of IT is defined as h = 1/{1 + exp[(V + 83.0)/4.0]} with the inactivation time constant h = 
8.2 + {56.6 + 0.27 exp[(V + 115.2)/5.0]}
/{1.0 + exp[(V + 86.0)/3.2]}. Only thalamic (Tp and Ts) and reticular (Rp and Rs) cells incorporated IT channels. This current combined with Ih (described above) endowed thalamic relay cells with intrinsic bursting properties. Although some cortical neurons contain T-type currents ( Paré and Lang 1998), we did not include them in model cortical neurons for the purpose of the present simulations. A slower T-current is known to exist in reticular neurons ( Destexhe et al. 1996b), but was not modeled here, although it is not expected that this current would have a significant impact on the present results.
PERSISTENT SODIUM CURRENT INA(P).
This sodium current is found in virtually all cortical neurons ( French et al. 1990; Kay et al. 1998; Mittmann and Alzheimer 1998; Stafstrom et al. 1984). It activates quickly near the resting potential and is considered persistent because it inactivates very slowly (on the order of seconds). We borrowed the formulation for INa(p) from previous work ( Compte et al. 2003; Fleidervish et al. 1996). The steady-state values for activation are used because they are considered to be instantaneous given their rapid time course (<1 ms). The steady-state activation is m(V) = 1/[1 + exp(–V + 55.7)/7.7]. We did not model the inactivation of this current given its very slow time course. All cells in the model contained INa(p).
DEPOLARIZATION-ACTIVATED POTASSIUM CURRENT IDK.
A Na+- or Ca2+-activated K+ current appears to play an important role in the termination of the depolarized phase of the slow oscillation ( Sanchez-Vives and McCormick 2000; Steriade et al. 2001). Both Na+ and Ca2+ currents are activated by the influx of ions that build up during periods of depolarization or spiking. To reduce the computational burden, we did not explicitly model either Ca2+ influx during spiking or the intracellular Na+ concentration. These concentrations increase most when the membrane potential is elevated. Therefore we chose to model this current as a generic activity-dependent potassium current with some activation properties taken from models of IKNa ( Wang and Lambert 2003). To characterize this depolarization-dependent influx, we use a sigmoid threshold function to determine how much the measure of depolarization D should increase for the current membrane potential. The factor D accumulates with depolarization and decays to the internal equilibrium concentration according to dD/dt = Dinflux – D · (1 – Deq)/D where Dinflux = 1/{1 + exp[–(V – D)/
D]}. The voltage-dependent influx Dinflux is determined by a sigmoid function with a threshold D = –10 mV and slope D = 5.0. The equilibrium level for the depolarization-dependent value is Deq = 0.001, and D = 1.25 s is the time constant that determines the return to Deq. The depolarization-dependent activation of the current IDK is given by m = 1/1 + (d1/2D)3.5. The parameter d1/2 = 0.25 determines the level of D necessary for half activation. All cortical (L5–6, L4, and L2–3; excitatory and inhibitory) cells contained IDK channels.
INFLUENCE OF DIFFUSE NEUROMODULATORY SYSTEMS.
Under physiological conditions, ascending neuromodulatory projections modulate the mode of firing in the thalamocortical system. Ascending neuromodulatory projections from several brain stem nuclei and the basal forebrain activate muscarinic, noradrenergic, serotoninergic, histaminergic, and glutamate metabotropic receptors, which modulate various cellular conductances that influence the overall level of depolarization on which the sleep-wake cycle critically depends ( McCormick 1992). As will be specified in the RESULTS, the various actions of neuromodulators are modeled as simultaneous changes of the conductances for IKL, Ih, IDK, INa(p), and AMPA synapses of the cortex, thalamus, and reticular nucleus.
SPONTANEOUS ACTIVITY: OPTIC NERVE FIRING AND MINIS.
The primary source of noise in the model was random spontaneous optic tract firing (45 spikes/s) modeled as 1,600 separate Poisson processes, which were independent of the behavioral state ( Mukhametov et al. 1970). The simulated optic nerve cells connect to Tp by diffuse projections with independent Gaussian latencies on each connection. This decreases the degree of synchronicity of spontaneous input and produces the slightly overlapping receptive fields observed in the LGN ( Kara and Reid 2003). This random activity percolated throughout the network, producing irregular spontaneous activity in all layers of the model.
In addition, we modeled "minis," the spontaneous release of neurotransmitter quanta ( Vautrin and Barker 2003), as low-amplitude PSPs (mean = 0.5 ± 0.25 mV) consistent with experimental observations ( Timofeev et al. 2000). The mean frequency of these Poisson distributed synaptic minis was set to 2 Hz (total for an individual cell).
Data recording
All state variables of the network (Vm, intrinsic and synaptic conductances and currents) were recorded during each simulation. These recordings were then used to visualize the model activity in a way that allowed comparison with experimental data—including local field potentials, optical dye recordings, and intrinsic and synaptic channel conductances.
The local field potential (LFP) as recorded in vivo is thought to be primarily a reflection of the net synaptic activity (i.e., not the mean firing rate) within a local region (several millimeters) around the measuring electrode ( Logothetis et al. 2001). A signal meant to resemble the local field potential was calculated as a spatial average of membrane potentials, 
, for all cells i within a given radius (unless otherwise specified, the radius = 20 units = the area of an entire model layer).
Voltage-sensitive optical dye recordings provide a technique to visualize spatiotemporal activity patterns in large-scale populations on a rapid timescale (about 10 ms) ( Fitzpatrick 2000). Accordingly, we visualized large-scale activity patterns by displaying average membrane potential (10 ms) while preserving topographic relationships between neurons.
Conductances of individual channels can be measured experimentally using patch-clamp techniques. In the model, all conductances are explicitly calculated, and therefore easily recorded under all conditions.
Simulation techniques
SYNTHESIS. All simulations were performed using a general-purpose object-oriented interactive neural simulator called Synthesis written by S. Hill (www.infinitedegrees.info). Synthesis provides a complete simulation environment, including: a computation server capable of parallel and distributed computation, a graphical user-interface for interactive visualization and manipulation, a scripting language for automating parameter searches and experiments, an interactive command-line environment for controlling the simulation, data agents for data gathering and analysis, and a library of standard neuron, synapse, and connection pattern models. The simulation server is multithreaded for multiprocessor computation and capable of distributing a simulation across multiple networked computers. Synthesis allows the user to interact fully with a simulation, visualizing and recording all variables—at different levels and timescales—while observing all network interactions and controlling all parameters in real time throughout the course of a simulation.
NUMERICAL METHODS.
Differential equations were numerically integrated using the Runge-Kutta 4th-order method ( Press et al. 1992) with a step size of 0.25 ms. The model was tested at smaller time steps with no significant differences observed. The generation of optic nerve activity and probabilistic connectivity patterns were based on standard pseudorandom number generator routines ( Press et al. 1992). Analysis of the simulation data was carried out using standard toolboxes in MATLAB (The MathWorks, Natick, MA).
COMPUTATIONAL TIME. The simulations were carried out on dual-processor 2.0 Ghz Power Macintosh G5 machines running Mac OS 10.3 and equipped with 3.5 GB RAM (Apple, Cupertino, CA). Each simulation of the full model used approximately 1.7 GB RAM. Computational performance varied with the mean firing rate and ranged from 500 to 700 ms/h for the full model. A single simulation run of 3 s in duration required over 5 h to compute.
|
RESULTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
). The full model contains 65,400 spiking neurons with 4,860,450 connections organized in regions and pathways, consisting of a primary and a secondary area of visual cortex (Vp and Vs), 2 corresponding regions of the dorsal thalamus (Tp and Ts), and 2 regions of the reticular thalamic nucleus (Rp and Rs). Vp represents a restricted portion of area 17 in the cat (i.e., about 1 cm2) and it contains units with small receptive fields that are selective for oriented segments. Versus represents a corresponding part of an extrastriate area with coarser topography, containing units with larger receptive fields that are selective for oriented lines and for line crossings. Tp and Ts correspond to a portion of the lateral geniculate nucleus (LGN) and to a portion of the lateral posterior pulvinar complex (LP), respectively. The thalamocortical model builds on a previously published architecture that was used to explore properties of synchrony and spike timing in response to visual stimuli and that successfully reproduced several experimental results ( Lumer et al. 1997a, b). A central feature of the model is the subdivision of the simulated cortex into 3 layers with different patterns of afferent, efferent, and local connectivity corresponding to supragranular layers, infragranular layers, and layer 4 (Figs. 1 and 2). Another key feature is the detailed simulation of horizontal intralaminar connections, vertical interlaminar connections within Vp and Vs, forward and backward connections between Vp and Vs, thalamocortical and corticothalamic connections, and connections from thalamic relay nuclei and cortex to the nucleus reticularis. We believe all these features constitute the minimum necessary components of a prototypical thalamocortical system.
Individual cortical and thalamic neurons, both excitatory and inhibitory, were modeled as single-compartment integrate-and-fire units using cellular constants from regular- and fast-spiking neurons, respectively ( Connors et al. 1982). Intrinsic currents were modeled using details of intracellular channel dynamics including those regulating hyperpolarizarion-activated cations Ih, low-threshold calcium IT, persistent sodium INa(p), potassium leak IKL, and a depolarization-dependent potassium current IDK. Synaptic interactions occurred through simulated channels that provided voltage-dependent (NMDA-like) and voltage-independent (AMPA-like) excitation, as well as fast (GABAA-like) and slow (GABAB-like) inhibition. All connections were endowed with conduction delays. Finally, neuromodulatory influences, such as those attributed to acetylcholine, were modeled as diffuse changes of intrinsic and synaptic conductances.
In what follows, we first show that our large-scale model of the thalamocortical system reproduces various aspects of spontaneous activity during wakefulness including low-voltage fast activity in the EEG, irregular firing, correlated subthreshold activity that reflects the functional architecture, and selective response to stimuli including gamma frequency synchronization in the evoked response. We then show that the model transitions to slow-wave sleep primarily arise from an increase in the potassium leak (IKL) current. The slow oscillation that subsequently emerges has many properties that are consistent with experimental results including a bimodal membrane potential distribution, a disfacilitated and silent down-state and a depolarized, high-conductance up-state that exhibits gamma frequency synchronization and mean firing rates in the range of 10–20 Hz. By performing several manipulations on the key model parameters, we investigate which intrinsic and synaptic currents underlie the initiation, maintenance, and termination of the slow oscillation. Finally, we demonstrate that the synchronization of the slow oscillation in the model is dependent on corticocortical connections, consistent with experimental observations.
The robustness of the results presented here was tested by replicating the experiments described below while systematically varying the relevant parameters as well as by starting from different initial conditions (the precise location of connections are generated probabilistically; see METHODS). For all simulation results described below, small perturbations in parameter values did not significantly alter the results (actual results of robustness tests are not shown because of the extremely large size of parameter space explored, amounting to more than 3 years of CPU time).
In the waking mode, the model exhibits spontaneous activity throughout the cortex and shows selective responses to visual stimuli
In the waking mode, spontaneous activity originating in retinothalamic afferents spreads throughout the thalamocortical system (Fig. 3, A–D, left). When a visual stimulus is applied, a selective response emerges from the spontaneous activity and high-frequency firing occurs for the duration of the stimulus (Fig. 3, A–D, right).
|
A moving visual stimulus, consisting of a 2-dimensional vertically oriented grating, was presented by firing a patterned subset of neurons in the optic tract for a duration of 250 ms. During the stimulus presentation, the vertically selective cells display a dramatically increased firing rate (10x) and high-frequency synchronized oscillations are clearly observable in the membrane potential (Fig. 3A, right). Poststimulus offset responses consist of marked inhibition and decreased firing, while the network reconfigures and resumes generating spontaneous activity after a few dozen milliseconds.
Intracellular unit recordings (Fig. 3B) depict spontaneous and evoked activity comparable to intracellular recording in vivo. Specifically, cortical excitatory and inhibitory neurons fire irregularly at low rates (2–10 Hz) during spontaneous activity, while firing at elevated rates (20–100 Hz) during stimulation, consistent with experimental data ( Azouz and Gray 1999). In these traces, the evoked activity appears to end prematurely because the stimulus is a slowly drifting grating and therefore does not rest within the receptive field of the individual cells shown for the entire time period.
The LFP (see METHODS), reflecting spontaneous and evoked activity, shows low-voltage fast-activity patterns in the absence of stimuli (Fig. 3C, left) and clear evoked responses with high-frequency oscillations during stimulus presentation (Fig. 3C, right), consistent with LFPs recorded in vivo during stimulation ( Gray and Singer 1989).
The topographic displays (Fig. 3D) show time-averaged (10 ms) membrane potentials in L2–3. Spontaneous activity in the model has a subthreshold correlational structure that reflects the orientation preference of the neural population (Fig. 3D, left). Specifically, subthreshold vertical and horizontal bars of depolarization are visible in the average membrane potential for each selective population even when no stimulus is applied. The spontaneous emergence of correlated activity reflects the underlying orientation-selective mechanisms, consistent with optical dye recordings in the monkey showing that cellular membrane potential fluctuations during spontaneous activity reflect the functional architecture and selective response mechanisms of visual cortex ( Kenet et al. 2003; Tsodyks et al. 1999). When a vertical grating is presented, vertically selective cells respond preferentially, whereas cells selective for horizontal features remain virtually silent, as is evident in Fig. 3D (right). A movie of the topographic view of several layers throughout the model network during both spontaneous activity and during stimulus presentation is available (see supplementary video #11 ).
An increase in potassium leak conductance triggers the transition from the waking mode to the sleep mode
In the model, the transition from the waking mode to the sleep mode (Fig. 4) comes about primarily by increasing the potassium leak conductance gKL (from 1.0 to 1.85), which is present in all model neurons. This corresponds to the unblocking of background potassium leak channels arising from the reduced actions during sleep of neuromodulators such as acetylcholine ( McCormick 1992). Without this single change of gKL, the network does not enter the hyperpolarized, silent state necessary for the sleep mode. A number of other network parameters are modulated during the transition from wakefulness to sleep. Several studies suggest that muscarinic receptor activation can significantly inhibit the INa(p) conductance ( Mittmann and Alzheimer 1998). The removal of acetylcholine would therefore cause an effective increase in the conductance of INa(p) throughout the network. We model this by increasing the conductance gNa(p) for the sleep mode (from 0.5 to 1.25). Acetylcholine has also been shown to shift the activation curve of Ih to more depolarized levels ( McCormick et al. 1993). We model this by increasing the conductance of Ih during the sleep mode (from 1.0 to 2.0). IT is also known to be inhibited by activation of muscarinic receptors ( McCormick 1992). We model this by increasing the peak conductance of IT from wakefulness to sleep (from 1.0 to 1.25). Neither of the changes to Ih or IT played a significant role in the transition to sleep. The removal of acetylcholine and norepinephrine unblocks slow potassium currents ( McCormick 1992), which are represented in the model by the depolarization-activated potassium current IDK. We increased gDK (from 0.5 to 1.25) in parallel with gKL. Finally, muscarinic receptor activation reduces intracortical EPSPs ( Gil et al. 1997), suggesting that—during wakefulness—excitatory synapses are depressed relative to slow-wave sleep. We model this change by increasing the amplitude of AMPA EPSPs (by 50%) from the waking mode to the sleep mode.
|
Figure 4B shows the same transition from the irregular firing of the waking mode to the sleep mode in 2 neighboring excitatory and inhibitory cells in L2–3 as well as single cells in the model thalamus and reticular nucleus. It is evident that the cortical cells undergo up- and down-states almost simultaneously. The reflection of slow oscillation activity is also seen in cells located in both Tp and Rp. The intracellular traces of these cells are comparable to those recorded in vivo ( Fuentealba et al. 2004; Steriade et al. 2001; Timofeev and Steriade 1996).
As revealed by the LFP (Fig. 4C), when the network transitions from the waking mode to the sleep mode, the population activity changes from irregular activity (low-voltage fast activity) to a synchronized oscillation that continues indefinitely (high-voltage slow activity).
Figure 4D depicts a topographic map of average membrane potentials (10 ms) from L2–3 during periods of the waking mode and during up- and down-states during the sleep mode. The topographic map during waking (Fig. 4D, left) is similar to that shown in Fig. 3D. The 2 topographical maps taken from the sleep mode illustrate the dramatic difference between up- and down-states. The up-state is clearly more depolarized than the down-state and neurons are intensely active. During the down state, neurons are hyperpolarized and virtually silent (Fig. 4D, right). The rapid development of this synchronization in the network is consistent with observations of the transition from wakefulness to natural sleep in cats ( Steriade et al. 2001). A movie of the topographic view of several layers throughout the model starting from the waking mode and showing the transition to the sleep mode is available (see supplementary video #2).
Additional simulations (not shown) reveal that when the network is transitioned more gradually from the waking mode to the sleep mode, and when the parameters influenced by neuromodulators are at an intermediate level, it exhibits sleep spindles. Spindle activity will be the object of a future publication and will not be discussed here.
In the sleep mode, the model displays a stable slow oscillation consisting of an up-state and a down-state
In the sleep mode, the network oscillates between a depolarized and hyperpolarized phase in a synchronized fashion (Fig. 5). The depolarized "up-state" lasts about 300–600 ms and is characterized by an average membrane potential of around –58 mV, elevated firing rates, and elevated intrinsic and synaptic currents. The hyperpolarized "down-state" lasts between 400 and 600 ms and is characterized by an average membrane potential of around –75 mV, an absence of firing, and very low intrinsic and synaptic currents. The LFP (Fig. 5, top) depicts the population level synchronization. Membrane potential rasters depict activity for each layer showing the synchronization of the up- and down-states across all layers and all neurons of the model cortex (Fig. 5, rasters). The network produces regular alternations between up- and down-states at about 1 Hz. Intracellular traces reflect the irregularity and variety of up- and down-states typical of the slow oscillation (Fig. 5, A–C, intracellular traces).
|
