The deceptively simple laminar structure of neocortex belies the complexity of intra- and interlaminar connectivity. We developed a computational model based primarily on a unified set of brain activity mapping studies of mouse M1. The simulation consisted of 775 spiking neurons of 10 cell types with detailed population-to-population connectivity. Static analysis of connectivity with graph-theoretic tools revealed that the corticostriatal population showed strong centrality, suggesting that would provide a network hub. Subsequent dynamical analysis confirmed this observation, in addition to revealing network dynamics that cannot be readily predicted through analysis of the wiring diagram alone. Activation thresholds depended on the stimulated layer. Low stimulation produced transient activation, while stronger activation produced sustained oscillations where the threshold for sustained responses varied by layer: 13% in layer 2/3, 54% in layer 5A, 25% in layer 5B, and 17% in layer 6. The frequency and phase of the resulting oscillation also depended on stimulation layer. By demonstrating the effectiveness of combined static and dynamic analysis, our results show how static brain maps can be related to the results of brain activity mapping.
The laminar structure of neocortex has been known for over a century, but cellular properties and the synaptic connectivity of the microcircuit are only now being elucidated. While both conventional electrophysiological recordings that characterize intrinsic properties of neurons (Chen & Fetz, 2005; Dembrow, Chitwood, & Johnston, 2010; Hattox & Nelson, 2007), and anatomical methods that reveal structural aspects of connectivity continue to be useful (Kameda et al., 2012; Tanaka et al., 2011), much progress has been made in elucidating the circuitry of microcircuits using a variety of high-resolution electrophysiological and optical techniques (for reviews, see Bastos et al., 2012; Douglas & Martin, 2004; Thomson & Bannister, 2003). Paired recordings have revealed the detailed laminar organization of the pyramidal neuron network in somatosensory (Lefort, Tomm, Floyd Sarria, & Petersen, 2009) and prefrontal (Morishima & Kawaguchi, 2006; Wang et al., 2006) cortices. Glutamate uncaging has been used with laser scanning to map excitatory connections in visual (Dantzker & Callaway, 2000), barrel (Bureau, Shepherd, & Svoboda, 2004, Schubert et al., 2001), and auditory (Oviedo, Bureau, Svoboda, & Zador, 2010) cortices. Optical and optogenetic methods have aslo emerged as a tool for mapping circuits involving interneurons (Kätzel, Zemelman, Buetfering, Wölfel, & Miesenböck, 2010; Packer & Yuste, 2011).
Assimilation of anatomical and physiological information through computational simulations is required for understanding functional connectivity. For example, quantitative approaches have drawn on morphological reconstructions to explore circuit structure (Binzegger, Douglas, & Martin, 2004; Stepanyants & Chklovskii, 2005; Stepanyants et al., 2008). A variety of network modeling efforts are beginning to incorporate aspects of the detailed cellular and synaptic-level circuit information to constrain simulations of the dynamic behavior of cortical circuits (Ainsworth et al., 2011; Lang, Dercksen, Sakmann, & Oberlaender, 2011; Markram, 2006; Roopun et al., 2008; Vierling-Claassen, Cardin, Moore, & Jones, 2010; Wang, 2010). The simulations were built on a data set obtained using a variety of techniques, including anatomical and cell recording, as well as laser-scanning photostimulation with glutamate uncaging and optogenetic stimulation.
Simulation is a valuable complement to brain activity mapping. Effects of underdefined parameters can be tested in the simulation to determine critical parameters, which can then be sought in experiment. Simulation also has the advantage of providing access to information that must otherwise be inferred, enabling systematic analysis of functional connectivity and network dynamics. In some cases, simulation permits compensation for experimentally imposed limitations. These features make simulation a useful complement to experiment and can potentially lead to predictions that can be verified or falsified through empirical testing.
Simulations can be performed at multiple scales. Mean field models, whose state variables represent the averages of action potential firing across a large population of cells (Wilson & Cowan, 1972), have the advantages of speed and analyzability at the cost of ignoring spike timing and unit synchrony. They also do not take into account detailed connectivity patterns, making them inadequate to study the effects of specific wiring patterns in the microcircuit. Simulations, however, can take into account electrical and chemical activity in individual neurons, with detailed characterization of subcellular compartments such as dendrites and synapses. While closely following empirical details from pharmacology and genomics, the computational cost of such simulations is extremely high. In addition, many parameters remain poorly constrained by experiment. The current simulation design is at an intermediate level of abstraction, utilizing a large network of modified integrate-and-fire neurons that can be tuned to grossly produce various dynamical features of different neuronal subtypes.
We have developed a spiking unit model of primary motor cortex (M1), based on experimental data from a single laboratory, on this single area, from a single organism (mouse), in a restricted age range. In this way, we attempted to avoid the chimerical nature of many cortical circuitry models, including our own prior models, which are based on data drawn from a variety of studies, species, and regions. However, in cases where specific measurements were not available, the model was parameterized based on models used in previous studies (see section 2).
We show two complementary analyses of the resulting model: static graph-theoretic analysis and dynamical stimulation-activation simulation. We demonstrate structure-dynamic relations: one population of layer 5 cells, the corticostriatal (STR) population, shows strong centrality in graph analysis and can also be shown to play a major role dynamically. Modifying network connectivity allowed us to identify critical cell populations and connections. Sustained oscillations emerged that differ depending on the laminar location of stimulation. This study provides a preliminary account of possible input-output transformations that can be performed by the M1 cortical microcircuit.
2.1. Spiking Neuron Model.
Individual neurons were event-driven, rule-based dynamical units used in prior models (Kerr, Neymotin et al., 2012; Kerr, van Albada et al., 2012; Neymotin, Lee, Park, Fenton, & Lytton, 2011; Neymotin, Jacobs, Fenton, & Lytton, 2011; Neymotin, Lazarewicz et al., 2011; Song, Kerr, Lytton, & Francis, 2013). We include details here for the convenience of readers. Units include key dynamical features such as adaptation, bursting, depolarization blockade, and voltage-sensitive NMDA conductance (Lytton, Neymotin, & Hines, 2008; Lytton & Omurtag, 2007; Lytton, Omurtag, Neymotin, & Hines, 2008; Lytton & Stewart, 2005, 2006). Event-driven processing provides a faster alternative to network integration: a presynaptic spike is an event that arrives after a delay at postsynaptic cells; this arrival is then a subsequent event that triggers further processing in the postsynaptic cells.
Each cell had a membrane voltage state variable (Vm), with a baseline value determined by a resting membrane potential parameter (VRMP). After synaptic input events (dendritic or somatic activations of AMPA, NMDA, or GABA A), if Vm crossed spiking threshold (Vth), the cell would fire an action potential and enter an absolute refractory period, lasting . After an action potential, an after-hyperpolarization voltage state variable (VAHP) was increased by a fixed amount WAHP, and then VAHP was subtracted from Vm. Then VAHP decayed exponentially (with time constant ) to 0. To simulate voltage blockade, a cell could not fire if Vm surpassed the blockade voltage (Vblock). Relative refractory period was simulated after an action potential by increasing the firing threshold Vth by , where WRR was a unitless weight parameter. Vth then decayed exponentially to its baseline value with time-constant .
The model used four excitatory and six inhibitory cell types, with four distinct parameterizations (see Table 1): two parameterization sets for excitatory cells (E) and two for inhibitory interneurons (I). L2/3 pyramidal cells (E2), L6 pyramidal cells (E6), and L5 STR cells all showed adaptation and were parameterized as in our prior models (Chadderdon, Neymotin, Kerr, & Lytton, 2012; Kerr, Neymotin et al., 2012; Neymotin, Kerr et al., 2011; Neymotin, Lee et al., 2011). Dynamics for these cells provided a rising firing threshold under continuous activation comparable to current clamp, resulting in spiking adaptation. L5 SPI cells had a different set of parameters providing a nonadapting pattern (Suter, Migliore, & Shepherd, 2012; Yu et al., 2008; Miller, Okaty, & Nelson, 2008) by omitting the rising firing threshold. The low-threshold spiking (LTS) inhibitory cells (designated with L) were provided with a lower threshold than were the other inhibitory cells.
|Cell Type .||VRMP .||Vth .||Vblock .||.||WRR .||.||WAHP .||.|
|E2, STR, E6||–65||–40||–25||5||0.75||8.0||1.0||400|
|I2, I5, I6||–63||–40||100||2.5||0.25||1.5||0.5||50|
|I2L, I5L, I6L||–65||–47||100||2.5||0.25||1.5||0.5||50|
|Cell Type .||VRMP .||Vth .||Vblock .||.||WRR .||.||WAHP .||.|
|E2, STR, E6||–65||–40||–25||5||0.75||8.0||1.0||400|
|I2, I5, I6||–63||–40||100||2.5||0.25||1.5||0.5||50|
|I2L, I5L, I6L||–65||–47||100||2.5||0.25||1.5||0.5||50|
Notes: VRMP = resting membrane potential (mV); Vth = threshold voltage (mV); Vblock = blockade voltage (mV); = absolute refractory time period (ms); WRR = relative refractory weight; = relative refractory time constant (ms); WAHP = after-hyperpolarization increment (mV); = after-hyperpolarization time constant (ms). Cell types: E2, layer 2/3 pyramidal neurons; STR, corticostriatal neurons; SPI, corticospinal neurons; E6, layer 6 pyramidal neurons; I2/I5/I6, fast spiking interneurons in layers 2/3, 5, and 6; I2L/I5L/I6L, low-threshold-spiking interneurons in layers 2/3, 5, and 6.
2.2. Network Connectivity.
The cell types were distributed uniformly within the depth range corresponding to their layer. Location 0 was at the surface (pia mater). The maximum location 1412 m represents the start of the white matter. The 775 cells were divided as follows: L2/3 cells distributed from 0 to 451.8 m: 150 E2, 25 I2, 25 I2L cells; L5A (451.8–663.6 m): 54 STR, 10 I5, 14 I5L cells; L5B (663.6–1059.0 m): 72 SPI, 113 STR, 30 I5, 26 I5L cells; and L6 (1059.0–1412.0 m): 192 E6, 32 I6, and 32 I6L cells. The relative numbers of cells were chosen according to estimated cell density measured as a function of depth from pia in vibrissal motor cortex (Hooks et al., 2011). Note that SPI cells are present only in L5B, whereas STR cells are present in both L5A and L5B.
Cells were connected probabilistically with connection densities and synaptic weights being assigned according to pre- and postsynaptic cell types (see Figure 1). For the I E and I I connections, we utilized densities from prior models (Kerr, Neymotin et al., 2012; Neymotin, Kerr et al., 2011; Neymotin, Lee et al., 2011), removing interlaminar IL E connections. For E E connections, and several E I connections, we calculated densities from experimental results (Apicella, Wickersham, Seung, & Shepherd, 2012; Kiritani, Wickersham, Seung, & Shepherd, 2012; Weiler, Wood, Yu, Solla, & Shepherd, 2008). The wiring matrix from Weiler et al. (2008) was used to estimate the strength of most of the E E connections. Kiritani et al. (2012) provided the important constraint that there were recursive STR STR and SPI SPI connections, and a significant STR SPI, connection, but no SPI STR connection. Apicella et al. (2012) justifies an interlaminar connection of E2 I5L and rules out E2 I5, STR I5L, and SPI I5L connections.
Synaptic weights Ws were tuned for cells of a particular projection type so as to provide reasonable firing rates postsynaptically (see Figure 1B). Weights did not change during simulation (no plasticity). These weights were assigned to the projections to dendritic AMPA synapse channels for E cells, and dendritic NMDA synapse channels were assigned a weight 1/10th these AMPA values. For fast-spiking I cells (I2, I5, I6), the weights in the figure were assigned to the somatic GABA synapse channels, whereas for the LTS I cells (I2L, I5L, I6L), the weights were assigned to the dendritic GABA synapse channels. Again, for the I E and I I connections, we borrowed from the topology of our previous models (Kerr, Neymotin et al., 2012; Neymotin, Kerr et al., 2011; Neymotin, Lee et al., 2011).
For E2 STR, E2 SPI, E2 I5L, weight was dependent on depth location, using data from Figure 3F of Apicella et al. (2012). For E2 STR connections, weights to L5A STR cells decreased with the cell depth of the E2 source cells, but weights to L5B cell depths were set to zero, as connections into L5B are not present in the data (see Figure 2 of Anderson, Sheets, Kiritani, & Shepherd, 2010). E2 SPI connections peak in weight at the middle of the L2/3 layer, and for both L5A and L5B targets, E2 I5L weights fall off with increasing depth of the source E2 units. Full connectivity details are available in the simulation code available on ModelDB (see below).
Our simulations corresponded to an in vitro slice preparation with no ongoing input activity driving the cells. The simulation was therefore quiescent until stimulated. For stimulation, a certain percentage of cells, both excitatory and inhibitory, in the selected layer were brought to firing threshold. The choice of which cells were activated was randomized; in each case, five different seeds were tested to make sure that a given pattern of activation was not dependent on a particular pattern of activation. Activity patterns were also tested for robustness with five different randomizations of detailed wiring (see Figure 5). For each of the four layers stimulated, 40 different stimulation percentages were tried. Thus, the primary simulation data set was on n=4000.
The model was implemented in NEURON 7.2 (Carnevale & Hines, 2006) for Linux. Full code is available on ModelDB (https://senselab.med.yale.edu/modeldb). One Second of simulation took approximately 30 s of CPU time on an Intel Core i5 CPU on a system with 2.8 GB of memory. Network analysis was performed using Networkx python package.
3.1. Static Network Analysis: Graph-Theory Measures.
We applied two graph-theoretic measures to identify central cell populations and connections: node betweenness centrality for cell types and edge betweenness centrality for connections between cell types (Freeman, 1977). For any two cell populations, there exists a shortest path along which activity propagates through the network. Nodes that are present in a larger proportion of shortest paths between all possible cell type pairs have a higher betweenness centrality value. That is, a population will have a high value if information between several other populations in the network is relayed through it. This measure is important because activity in the whole network might be modulated by controlling activity in a highly central population. This makes connections in and out of this population highly important. A similar measure, edge betweenness centrality, is used to assess the relative importance of connections in a network. A connection with a larger proportion of shortest paths passing through it has a higher edge betweenness centrality. It is not dependent on connection sign (excitatory or inhibitory) or strength as static analysis is used in this study primarily to identify potentially important populations and connections. Hence, a connection might be strong but not central, and vice versa. Subsequently, simulations are performed to analyze the dynamic instantiation of this network where the impact of connection strength and sign are taken into account.
Static analysis demonstrated a high betweenness centrality for the STR subpopulation of L5 (centrality 0.57; see Figure 2). STR sits solidly between two other centers with high centrality: E2 (0.41) and E6 (0.39). It is also well connected with a one-way conduit to the other L5 population, the SPI cells, which provide a major output downward to lower centers. STR therefore appears to serve as the critical relay of activity from L2/3 and L6 and a critical relay to downward-projecting outputs, as well as itself providing outputs to other areas of cortex and to striatum. STR thus occupies a highly central position in the wiring diagram owing to its numerous connections (several that are reciprocal) with other populations and more specifically due to its critical role acting as a bridge for excitatory information flow between superior and inferior layers. The centrality of E2 and E6 indicates that they will play important roles as signal relays in managing inputs to their respective layers.
The most central excitatory edges are STR E2 (0.25) and STR E6 (0.23), again emphasizing the centrality of STR and emphasizing this triad of centers as a likely backbone of signal flow. I5 STR is an inhibitory projection that also shows strong centrality (0.16), presumably helping to balance the excitatory influences from these two E populations. Associated with the SPI I5 projection (0.14), this pathway provides the only feedback from SPI to STR, given that the excitatory projection between these populations is unidirectional in the opposite direction. Other edges show far smaller centrality values in the 0.01 to 0.09 range.
Centrality analysis provided just one view of this complex network, limited in that it did not consider connection strengths. A certain connection might be central yet have less impact on network dynamics due to low strength. We therefore added another view, analyzing overall connection strength, defined as the product of weight and number of connections, normalized to account for variations in population size (see Figure 3). The STR population stood out clearly as a major hub in this analysis as well. The same major connections also still stood out: STR E2 (108.9) and STR E6 (114.7). Also of note is the moderately strong inhibitory backprojection from SPI to STR via I5, the one backprojection from SPI back into the microcircuit. Outside of this projection, SPI is entirely isolated as a pure output of the microcircuit.
Both E2 and E6 exhibit strong self-connections and strong reciprocal connections with the fast-spiking interneurons of their respective layers. Two features distinguish E2 from E6. First, E2 sends a much stronger (434.1) connection to SPI compared to E6 (115.2). Second, the strength of reciprocal connections to STR is less for E2 (E2 STR: 76.3 and STR E2: 108.9) when compared to E6 (E6 STR: 113.1 and STR E6: 114.7).
3.2. Dynamic Network Analysis: Activation Patterns Differ with Stimulus Location.
Subthreshold stimulation of a single layer produced firing activity that propagated to other layers but was not sufficient to trigger sustained oscillations (see Figure 4, left column). Activity died out typically within approximately 100 ms. With suprathreshold stimulation, activity spread to other layers and triggered sustained oscillations (see Figure 4, right column).
Layer 2/3 has powerful excitatory projections downward to the L5 populations (see Figure 3). We therefore anticipated that this layer would be particularly effective in activating the column. With 12% activation of L2/3, activity spread to layer 5 neurons but was not sufficient to cause sustained oscillations (see Figure 4A). With 13% activation, activity rapidly spread to layer 5 and from there to layers 2 and 6, subsequently triggering oscillations in the network (see Figure 4B). Spread of activity with subthreshold stimulation was as expected from the wiring. E2 activates SPI and STR populations at about the same time, but the response of the STR population is more sustained at suprathreshold stimulation, due largely to the reciprocal connectivity with E2. L6 is recruited only with suprathreshold stimulation. Interestingly, L2/3, the layer being stimulated, begins to participate in the sustained oscillatory response only after L5 is fully engaged. This delay in E2 activation also reflected the role of inhibitory cells. Initial activation of E2 activated local I2 inhibitory cells, which will delay recurrence of activity into L2/3, restricting any immediate response in E2 that might have otherwise emerged from E2E2 connections. (Although I2 cells were also directly activated by the stimulation itself, this activation was minimal and not significant compared to the activation via E2 projections causing feedback inhibition onto themselves.) By and large, E I projections are local to the layer. A demonstrated exception is the E2 I5L projection (Apicella et al., 2012).
Proceeding down the laminae, we next assessed stimulation to L5A (see Figures 4C and 4D), which contains STR (but not SPI) cells (Anderson et al., 2010). With 53% stimulation there was no sustained activity (see Figure 4C). Sensitivity to stimulation was substantially less than in L2/3: 54% stimulation was required for sustained activity (see Figure 4D). The initial response of the directly activated STR cells was paradoxically far less than seen with L2/3 stimulation due to the immediate powerful activation of L5 inhibitory cells in this setting (0–20 ms in Figure 4B versus Figure 4D). This inhibition also helps explain the high stimulation required to produce activation with L5A stimulation. However, a paradoxical corollary of this inhibition was the more prolonged activity in the transient stimulation case: SPI cells remained active for about 40 ms, about twice as long as seen with L2/3 stimulation (see Figure 4A). This sustained response reflects a brief interplay of excitation and inhibition involving SPI and I5 cells.
L5B consists of both SPI and STR neurons. With 24% L5B stimulation, sustained oscillations were not produced (see Figure 4E). This stimulation produced a primary SPI and I5 transient response with slightly less activation of STR than was seen with L5A stimulation (see Figure 4E versus Figure 4C). Sustained oscillations with 25% stimulation were also different, with activation of all layers happening after a greater delay (see Figure 4F versus Figure 4D). This difference can also be traced to the early activation and dominance of the SPI activity: SPI cells project exclusively outside the network and hence do not influence activity locally (except within the SPI population itself).
L6 stimulation produced a different activation pattern than seen with other stimulation sites. The transient response with 16% stimulation has less activity, with only minimal involvement of STR and SPI cells (see Figure 4G). This reflects similar E6 STR and E6 SPI weights (see Figure 1B) and no activation of E6 cells beyond the time of stimulation. A 17% stimulation showed sustained activity (see Figure 4H).
In order to assess the robustness of lamina-specific thresholds required to generate sustained oscillations, we performed simulations with different random seeds for the wiring and stimulus: five wiring seeds and five input seeds (see Figure 5). Median values were 15%, 50%, 30%, and 16% for L2/3, L5A, L5B, and L6 stimulations with similar mean values. The differences in layer sensitivity appeared to be due to several causes, including layer size and layer centrality.
Following the responses out to longer times (up to 120 s) demonstrated that the oscillatory responses at higher stimulation strengths were sustained indefinitely (see Figure 6). More coherent oscillations (i.e., the degree to which the firing phases are aligned across cells in a population) were produced by stimulation in L2/3 or L6. For example, the firing pattern of E2 population is more coherent in Figure 6A (L2/3 stimulation) when compared to the corresponding pattern in Figure 6B (L5A stimulation). Generally STR firing was less coherent than SPI firing, with the degree of coherence scaled according to that of the SPI population for a particular stimulation condition. SPI activity gradually lagged E2 and STR following L2/3 stimulation, but led E2 and E6 following L6 stimulation, where E2 and E6 were grossly synchronized. E6 and E2 were in phase following L6 stimulation and, less coherently, following L5A stimulation. E6 led E2 in the case of L5B stimulation.
In addition to differences in coherence, there were also differences in response phase and frequency with stimulation in the different layers (see Figure 6E). Measured population responses in the L5 population were about 60 Hz for most stimulation cases. However, in the case of L6 stimulation, frequency doubled to approximately 125 Hz (note that this localized effect reflected differences in initial conditions with no change in network parameters). This frequency doubling contrasts with the frequency doubling of the SPI population produced by L5B stimulation (see Figure 6C). The combined L5 activity masks the sharper population firing patterns of the SPI population alone. Phase of oscillations generated by L2/3 stimulation was earlier than those generated by L5B stimulation, with L5A stimulation showing the greatest relative phase lag.
3.3. PING: The Mechanism for Sustained Oscillations.
Using a subnetwork with STR and I5 neurons, we investigated whether sustained oscillations might be generated by a pyramidal-interneuron gamma or PING mechanism (Tiesinga & Sejnowski, 2009). This subnetwork consists of reciprocal pyramidal-interneuron connections as well as mutual connections within those populations (see Figure 3). With 60% activation, pyramidal neurons exhibited oscillations with a frequency of approximately 60 Hz (see Figure 7A). These oscillations in the gamma frequency range of 30 Hz to 100 Hz were sustained for up to 120 s (the longest duration for which simulations were run). The PING mechanism involves the interaction of excitatory and inhibitory populations, with the excitatory volley triggering an inhibitory volley. We tested this by cutting connections between STR and I5. Pyramidal neurons exhibited asynchronous activity (see Figure 7B), suggesting that oscillations were indeed generated through PING.
3.4. Impact of I5STR, E2I5L and SPISTR Connections on Network Dynamics.
We investigated the functional relevance of three specific connections in the network: I5STR (strongest connection in the network), E2I5L (only excitatoryinhibitory connection across layers), and the absent SPISTR connection (if present, would link the two major output neuron types).
After making specific modifications to the connectivity, we assessed oscillations generated by these modified networks in response to suprathreshold stimulation of L5A. Activity generated by an intact network was used as control (see Figure 8A) against which activity of modified networks was compared (see Figures 8B to 8F).
The impact of I5STR on activity was assessed. An increase in STR activity, relative to control, was expected with a halving of I5STR strength (see Figure 8B). STR activity not only increased but also became more dispersed. The frequency of STR oscillations increased from approximately 55 Hz to approximately 85 Hz. An unexpected and second-order effect was an increase in the E2 and E6 activity, especially approximately 150 ms after stimulation. While static analysis did not reveal a direct link between I5 and E2/E6 activity, simulations revealed that I5 indirectly controls activity in these populations. Completely cutting the I5STR connection counterintuitively resulted in highly suppressed STR activity (see Figure 8C). This suggests that increased STR firing resulted in its own inhibition through the STRE2I5LSTR path and, to a lesser extent, STRSPII5I5LSTR path. E2 and E6 activity became further dispersed. It was not possible to predict these results from the wiring diagram as they depend on precisely tuned activity of both excitatory and inhibitory populations. For a given stimulus and wiring diagram, simulations help identify specific connections and subnetworks that will be instantiated and become relevant to network activity.
We then investigated the impact of the E2I5L connection, which is unique by virtue of being the only interlayer excitatoryinhibitory connection in the network. Decreased excitatory drive to I5L resulted in decreased inhibitory drive to I5. This resulted in sharpened STR activity (see Figure 8D). In order to compensate for elevated I5 activity, we halved the strength of I5STR connection (see Figure 8E). The resulting dynamics exhibit characteristic properties of both modifications: an increase in STR oscillation frequency with more dispersed E2/E6 activity and sharpened STR activity. This suggests that I5STR is a critical determinant of STR oscillation frequency (and E2/E6 activity), while E2I5L connection determines the dispersion of STR activity.
Introduction of an excitatory SPISTR connection with the strength of E2STR connection caused little gross effect in steady-state network dynamics except for an increase in SPI cell firing (see Figure 8F). Similarly, cutting E2STR or E6STR had little gross effect on dynamics (data not shown).
We have built a network model of M1 cortex based directly on functional connectivity measures. Parameters were largely based on a consistent set of brain activity maps from a single area of a single species (mouse) over a narrow age window (young adult from approximately 1 to 2 months), with experiments done in a single laboratory. The consistent nature of the underlying experimental data allows for greater confidence in the robustness of simulation results. Realistic and lamina-specific oscillatory patterns observed in the current simulations were not observed in our previous models constructed using multiple data sets. Not only is this a methodological advancement, but it also makes the results and predictions of this study directly relevant to experimentalists who can look for predictions made for a particular brain area. It must be noted that due to the nature of our data set and model design, oscillations generated by the simulation correspond closely to an in vitro slice.
We made the following observations:
Lamina-specific activation is strongly thresholded and produces oscillations with distinct characteristics.
A decrease in the influence or firing rate of fast-spiking I5 activity increases the oscillation frequency of STR and shifts it toward higher frequencies within gamma. This is testable by using optogenetics to alter interneuron firing rates (Cardin et al., 2009).
The effects of modifying both E2I5L and I5STR can be observed simultaneously. It is possible to effectively tune dynamics of all excitatory populations by tuning the activation of I5 and I5L neurons.
Oscillations developed in L5 through a PING mechanism also have implications for oscillations generated in L2 and L6. The central location of STR, L5, and I5L populations makes them part of a critical subnetwork within the overall network.
The major conclusion of static graph-theoretic analysis of the circuit was that STR cells constitute a major hub in the M1 microcircuitry. This was demonstrated by both betweenness centrality and a simple connectivity-strength criteria. This finding was then confirmed by dynamical analysis. STR cells were generally second to be activated with activation from any layer (e.g., from L2/3 in Figure 4B). They also engage in the strong oscillatory responses that involve an interplay with the other populations (see Figures 6A to 6D). STR cells project to the SPI cells and also to L2/3 and L6, providing information flow both upward toward superficial layers and downward toward deep layers. They also pass information to contralateral striatum and form significant connections among each other. As such, they are positioned to receive a rich diversity of inputs and influence the activity of motor cortex and basal ganglia.
Our initial explorations of the dynamics of the model demonstrated complex responses to synchronous population stimulation that varied depending on the location and strength of the stimulation. These results are especially relevant in the context of experimental data, which suggests that long-range connections to the motor cortex are lamina specific (Hooks et al., 2013). For instance, L2/3 and 5A receive inputs predominantly from the sensory thalamus, L5B and L6 receive inputs predominantly from the secondary motor cortex, and L5A and L5B receive inputs mainly from the motor thalamus. The sensitivity of the model varied with laminar stimulus location, with L2/3 being most sensitive, followed by L6, L5B, then L5A. The L2/3 L5 connection was critical to the interpyramidal cell activation across cortical layers. Transient and sustained activity were supported by both recurrent interlayer activation loops and excitatory recurrent connections within the pyramidal cells of a particular layer. This is consistent with previous experimental work suggesting that different oscillatory patterns are generated by laminar-specific stimulation (Sakata & Harris, 2009; Beltramo et al., 2013). The underlying mechanism of the oscillation is pyramidal-interneuron network gamma (PING).
Previous studies have analyzed the static connectivity or modified static connectivity in order to explore functional consequences. Betweenness centrality measures have been used to understand epileptiform activity in both models (Morgan & Soltesz, 2008) and experimental studies (van Diessen et al., 2013; Wilke, Worrell, & He, 2011; Bernhardt, Chen, He, Evans, & Bernasconi, 2011). As our connectivity data were constrained by experiment, we used dynamic analysis not merely to explore consequences of static connectivity but as a way to gain greater insight into it. We did this by identifying dynamics that were not readily predicted by the wiring diagram. This allowed our simulations to effectively complement experimental studies by exploring the effect of network modifications that are not typically feasible in experiments. A novelty of this work is that we performed both static and dynamic analysis. Static analysis of the writing diagram was used to identify central populations and connections and predict interesting dynamics. Dynamic analysis was used to validate and extend these analyses. The utility of this combined study is illustrated by the following example. Diminishing the strength of I5STR connection expectedly increased STR activity. This confirmed results from static analysis. But the increased dispersion of E2 and E6 activity was revealed through simulation. Further, cutting this connection revealed entirely novel dynamics, where E2 and E6 activity become highly dispersed at the expense of much reduced STR activity. Hence, both static and dynamic analyses combined to provide a holistic view of network functioning.
I5 activity was a major determinant in STR oscillation frequency and dispersion of E2 and E6 activity. I5L activity driven by E2 was a major determinant in dispersion STR activity. These results are consistent with previous work that suggests a close relationship between excitatory and inhibitory cells (Isomura, Harukuni, Takekawa, Aizawa, & Fukai, 2009). We have also shown that the subnetwork with STR and I5 can generate oscillations using a PING mechanism. When SPISTR, E2STR or E6STR connections were cut, no gross effect was observed, suggesting that oscillatory activity generated in L5 critically determines oscillations in other layers but not vice versa.
This research supported by DARPA grant N66001-10-C-2008 (J.T.F.,W.W.L.), NIH/NINDS grant R01 NS061963 (G.M.G.S.); and NIH grant T32 NS041234 (B.A.S.); and Whitman Center, Marine Biological Laboratory (G.M.G.S., W.W.L.). We thank Larry Eberle (SUNY Downstate) for Neurosim lab support, Michael Hines (Yale) and Ted Carnevale (Yale) for NEURON support, and Tom Morse (Yale) for ModelDB support.