An explicit goal of our modeling approach was to include as much biological realism as possible, in particular at the level of neural connectivity, but also at the level of the individual neuron. The overall computer simulation consisted of 400 pyramidal cells and 52 interneurons. We wished to include a number of neurons that was large enough to reveal the emergent properties of the hippocampal tissue, but within the constraints of our computing hardware. We used the 64 compartment pyramidal cell model described by Traub et al [81
], and the 46 compartment interneuron model implemented by Menschik and Finkel, which is an adaptation of the Traub and Miles [83
] model. Both models feature realistic dendritic arborizations, and incorporate Na+
, and K+A
channels distributed along the somato-dendritic axis, as well as an explicit representation of internal Ca++
concentration. Axon initial segments were explicitly modeled as compartments, but axons themselves were operationalized simply as delays. The number of neurons in the overall computer simulation is of an order of magnitude similar to those of previously published computational models. The standard for adequate network size, however, is not an arbitrary number of cells—rather, an important criterion is that it exhibits network level behaviors seen in actual tissue; evidence of this was given in the Results section.
This model was explicitly conceptualized as a three-dimensional section of hippocampal tissue, extending 154 microns in the septo-temporal direction, 154 microns in the transverse direction, and 634 microns in the dimension orthogonal to these, extending from the stratum lacunosum-moleculare to the alveus. Each soma has an x, y, and z location in the 3-space defined by these axes ().
Figure 1 Schematic diagram of overall model. Grey spheres indicate location of cell bodies, black cubes show placement of electrodes. Boundaries of simulated tissue block are shown with black lines. Dense layer of cells is s. pyramidale, with s. radiatum and s. (more ...)
The model was constructed, and will be described, as a two-step process: first, calculating the steric relationship of the various cells of the model, and second, representing their patterns of connectivity.
Densities of interneuron subtypes by calcium binding proteins in the various hippocampal strata (stratum oriens, pyramidale, radiatum, lacunosum-moleculare) were calculated from data given in Freund and Buzsaki [29
] and Jino and Kosaka [45
]. Freund and Buzsaki’s study (their Figure 23) indicates areal densities for 60 micron hippocampal slices; from this, we calculated spatial densities (see ). The breakdown of model interneurons is as follows: 16 parvalbumin (PV) cells, 6 calbindin (CB) cells, 19 calretinin (CR) cells, 7 somatostatin (SOM) cells, and 4 cholecystokinin (CCK) cells. For all interneuron subtypes, we use the same dendritic morphology and dendritic ion channel distribution. While it is certain that these subtypes differ in terms of dendritic morphology and neurophysiology, we felt that this was a reasonable approximation, particularly given the large number of neurobiological unknowns involved for each of the subgroups. In the model, the calcium-binding subtypes are distinguished from one another on the basis of axonal projection characteristics, as described below.
Table 1 Densities of hippocampal cells in model. Figures in the Table are areal densities in soma/mm2; data are for 60 micron thick slices from rat hippocampus taken perpendicular to the septo-temporal axis. Interneuron data from Freund and Buzsaki , pyramidal (more ...)
Density of pyramidal cell bodies in stratum pyramidale of CA1 was taken from the neuroanatomic studies of Boss et al [11
] and Hosseini-Sharifabad and Nyengaard [42
], which provide similar numerical estimates. For each cell, x, y, and z coordinates within the relevant stratum were randomly generated.
A good deal about neuronal connectivity in hippocampus remains unknown. To take maximum advantage of what is known at a neuroanatomic level, we constructed axon-to-dendrite synaptic connections according to the following algorithm: (i) Map cell categories based on calcium-binding proteins onto categories based on morphology, to the extent possible. (ii) For each morphological class, define axonal projection patterns based on synaptic targets, in terms of stratum, target cell type (pyramidal cell vs. interneuron), and synaptic target area on cell (initial segment vs. soma vs. dendrites). (iii) Calculate spatial bouton densities for the axon projection clouds of the various subtypes, and apportion synapses according to percentages of (ii) above.
Interneurons Classification by Calcium Binding Protein vs. Morphology
At a summary level, interneuron connectivity was based on the information shown in , which represents the distillation by Freund and Buzsaki [29
] of a large number of neuroanatomic studies. The precise manner in which the interneuron categorization based on calcium binding proteins used here, and in much of the current literature, maps onto traditional categories based on morphology (e.g., basket cells, chandelier cells, etc) is only beginning to become clear. However, to the extent that relationships are known, we attempted to include them in the model. For example, as indicated in , PV interneurons project to stratum pyramidale and the proximal area of stratum oriens. It is thought that interneurons that stain for PV are almost entirely basket cells or chandelier cells, though the precise percentage of each is not known. Nonetheless, we felt that it would be more realistic to assume two distinct subgroups of PV interneurons—one projecting densely to pyramidal cell somata, the other to pyramidal cell initial segments—than to assume a homogeneous category of soma/IS-projecting interneurons.
Figure 2 Summary diagram of axonal and dendritic arborizations of interneurons, by subgroup. The dark circles indicate the location of cell bodies for each of the interneuron subtypes; the dark lines emanating from them show the orientation and laminar distribution (more ...)
The manner in which the calcium binding categories included in our model—which are largely non-overlapping—were allocated to one or more morphological categories is described below. All of the following is summarized in .
Table 2 Connectivity parameters. a. Morphological categories of interneuron subtypes. For cells of the given calcium-binding category in the given stratum, breakdown by morphological class is given. For example, the first row indicates that 45% of PV cells located (more ...) Parvalbumin
Many researchers have found hippocampal PV cells to be either basket cells or chandelier cells [e.g., 9
]. One group [66
] found a small number of a sample of 36 PV-staining interneurons in CA1 to have a bistratified axon projection pattern. It was therefore assumed that 10% of all PV cells are bistratified, with the remainder evenly divided between basket and chandelier cells.
Researchers have found that CB cells tend to innervate the proximal and distal dendritic trees of interneurons in CA1 [56
]. Freund and Buzsaki [29
] divide CB cells into three subtypes: (a) Cells with a bistratified axonal arborization, projecting to SR (stratum radiatum) and SO (stratum oriens), but respecting the SP (stratum pyramidale); (b) Radiatum-projecting cells. We assume these to be the “type I” cells described by Gulyas and Freund [31
]. Somata of these cells are said to lie predominantly in the SR, and only occasionally in SP and SL (stratum lacunosum-moleculare). (c) Horizontal cells at the SR-SL border. The axonal arborization of these cells is poorly characterized. Therefore, of the CB cells whose soma lie in SR, 50% were assumed to be bistratified and 50% were assumed to be radiatum-projecting. CB cells elsewhere were assumed to be bistratified, so that our categories would be entirely non-overlapping. Cells that stained for both CB and SOM were classified as SOM cells. Approximately 87% of CB-staining cells in SO (and 0% of cells in other strata) fall into this category [29
These interneurons primarily target the dendritic processes of other interneurons [25
]. Therefore, CR cells, regardless of location, are classified as interneuron-projecting cells (see below).
In CA1, SOM cells are said to correspond to oriens-lacunosum moleculare (o-lm) cells [48
]. This was taken as a one-to-one correspondence in our model.
There is an extensive literature suggesting that CCK-immunoreactive cells are generally basket cells [29
]. One labeling study [66
] found a minority of 27 CCK immunoreactive cells to be bistratified cells, Shaffer collateral associated cells, or of other non-basket cell morphologies. Given the preponderance of data, we assumed that all CCK cells in the model are basket cells.
We define interneurons’ morphological categories in terms of axon projection pattern. (As mentioned above, we use a single canonical model for interneuron physiology and dendritic morphology.) These data are generally taken from studies in which an axon of a particular interneuron type is labeled with an anterograde tracer, allowing visualization of the entire axon with all of its ramifications. For example, Katona et al [48
] found that o-lm cells project only to SL and that of the total number of synaptic connections they form there, 86% are on dendrites of pyramidal cells, 3% are on pyramidal cell somata, and 11% are on interneuron somata; quantitatively similar results were obtained by Matyas et al [56
]. Analogous data were derived from comparable experiments on the other interneuron classes: basket cells [33
], bistratified cells [31
], oriens-lacunosum-moleculare cells [48
], and pyramidal cells [79
In terms of computational implementation, we chose a data structure that was consistent with the level of detail available in the neuroanatomic literature. Thus, as shown in , for each neuron type we are able to specify the stratum to which it projects, whether the target cell is a PC or interneuron (but not which subtype
of interneuron), and the region on the target cell on which synapses fall (IS, cell body, or dendrites). In rare cases, however, data at a higher degree of resolution are available. For example, Gulyas et al [32
] found that CR interneurons synapse with dendrites of other interneurons, with the following caveats: (a) CR cells avoid PV interneurons; and (b) when CR interneurons contact other CR cells, they do so at both dendrites AND somata. Given the data structure of the model, it was not possible to include these facts in the current simulation.
Spatial Bouton Density
Quantitative estimates of spatial synaptic density for each of the interneuron classes were derived from the work of Sik et al [76
]. These researchers anterogradely stained axonal projections of each of a number of different classes of interneurons. Making several tissue sections, they calculated a total (linear) axonal distance per volume, as well as a bouton frequency, in boutons per length of axon. From these data, a spatial bouton (synapse) density can be calculated. In our model, we made the following assumptions: (i) for a given cell, bouton density depends on the particular type of interneuron, as defined by the type of calcium binding protein expressed by the cell. Breakdown of synaptic targets—or, cast another way, probability that a bouton of a given cell in a given stratum will contact a particular target—is given by the figures of , as described above. (ii) Spatial bouton density does not vary with septo-temporal or transverse distance from the parent cell. While Sik et al show that densities do decrease as distance from the soma increases to, for example, a mm or more, given the small scale of our model, this is likely a safe assumption. (iii) For a given interneuron type bouton density does not vary from stratum to stratum. The values used in the model, in units of synapses per micron3
are as follows: PV 2.12 × 105
, CB 1.13 × 104
, CR 6.36 × 103
, SOM 4.21 × 105
, CCK 2.12 × 105
Neuroanatomic literature [14
] suggests that CA1 axons course through striatum oriens toward alveus, giving off occasional collaterals. In our model, we assume that pyramidal to pyramidal connections occur in stratum oriens (that is, on cells’ basal dendrites). This should be contrasted with area CA3, where it is well known that recurrent excitatory connectivity is considerably denser.
It should be noted that postsynaptic targets, as percentages, are not known with a high degree of certainty for pyramidal cells [77
]. Somogyi et al [79
] state that CA1 pyramidal cells synapse on other pyramidal cells as well as a basket cells and SOM cells; synapses onto axo-axonic, bistratified, and other interneuron types are thought to be possible, but not yet proven. We therefore assume projection percentages (PCs vs. interneurons) in rough proportion to their presence in SO.
Using a methodology similar to that described for interneurons, spatial density for varicosities in pyramidal axons was calculated from data given in Perez et al [67
] and Knowles and Schwarzkroin [49
]. This was calculated to be 9.09 × 104
synapses per micrometer. It is well known that pyramidal to pyramidal connectivity is considerably lower in CA1 than in CA3, though few quantitative estimates of this have been made.
While we aimed to include a high level of biological realism in our model, particularly in the area of connectivity, certain simplifying assumptions were made. All computational models necessarily make abstractions. Which details to include and not include must be based on the purpose of the computational model in question; it is best to explicitly acknowledge the assumptions made.
We used a pyramidal cell model designed as a CA3 cell, and there are some differences between the electrophysiological response properties of CA1 and CA3 PCs. For example, in contrast to CA1 cells, CA3 cells show action potentials that decrease in amplitude as repetitive spiking occurs. Also, it is thought that CA1 cells show greater SFA than CA3 cells. Thus, while this may introduce a degree of error in the model, given the direction of the error, it is tolerable. That is, our model may understate the degree of SFA at play in CA1, and it is possible that the SFA-mediated effect we describe could in actuality be greater.
Also, our model does not include gap junctions. Traub et al [84
] created a network model incorporating gap junctions between dendrites. They show that these junctions underlie oscillatory activity, particularly in the gamma range. However, the primary focus of our study was the elucidation of system wide stable states and their intermittent transitions over several seconds. Gamma range (40 cycles per second) activity operates at a markedly different time scale, and we felt safe not including this, given our behavior of interest.
Synaptic conductances were assumed to obey a dual exponential function, as follows:
where A is a normalization constant chosen such that gsyn
reaches a maximum of gmax
We assumed all interneurons form GABAergic synapses on their target cells. We assumed that pyramidal cells form excitatory synapses on their target cells, either AMPA or NMDA, with a 0.5 probability for each. Channel characteristics that were used in our model, in terms of values given in Equation 1
, are presented in .
Table 3 Synaptic channel parameters. Column headings: E equilibrium potential; τ1, τ2, and gmax as given in Equation 1. Subscript p indicates postsynaptic cell is pyramidal, subscript i signifies interneuron.
Because we were interested in understanding behaviors that were properties of the tissue itself, our model is intentionally self-contained: all synaptic stimulation to cells (aside from the initial stimulation) arises from the other cells of the model. A by-product of the fact that our model represents a very small piece of tissue is that the amount of innervation received by a given cell is considerably lower than that received by a hippocampal cells in vivo
: for example, actual CA1 PCs receive about 30,000 excitatory and 1,700 inhibitory inputs [59
]; the number of synaptic inputs received by PCs in our model were orders of magnitude lower than this. To compensate, we multiplied gmax
for all synaptic conductances by a constant scaling or “weight” factor. We arrived at this multiplier by gradually increasing it and re-running the model until pyramidal interneuron transmembrane voltages showed realistic-appearing spike waveforms; the factor we used was 300. We were reassured in our choice of weight factor by subsequent statistical analysis (see Results) that indicated reasonable agreement between basic model behaviors and actual tissue.
Model Inputs and Outputs
In contrast to some modeling efforts, our intention was not to impress upon the network patterned stimuli implicitly containing information, or the output of a single neuroanatomic track. Rather, we wished to present the network with an input stimulus that was biologically realistic yet “neutral”, with the emphasis on understanding how this tissue processed and transformed the information. Thus, we placed a “dummy” AMPA channel in each cell soma of the model (D. Beeman, personal communication, May 2002), and activated each randomly, with a mean activation rate of 15 Hz; this stimulation was applied for one second, and the model was allowed to run for an additional five seconds of brain time. Transmembrane potential from the soma of each neuron was recorded.
The question of initial conditions is a difficult one, as any set of assumptions about starting values for transmembrane potentials, values of Hodgkin-Huxley channel gates, etc., introduces bias into the functioning of the model; ideally, one would select a starting state that is stochastic, but biologically plausible. Therefore, we allowed the system to start in a state in which all variables were set to 0; we then applied the random stimulation as described above. This produced unnatural behaviors briefly which damped out within about 0.4 second, and spike waveforms became normal-appearing. We took this 0.4 second mark as the starting point of the simulation, given that it would make for stochastic, but still neurophysiologically realistic, initial conditions. Data from the preceding “transient” period was excluded from our analysis.
We also placed “virtual electrodes” in 16 locations throughout the simulated tissue block to record local field potential (). We calculated LFPs using the following equation, given by Nunez [64
], and adapted for use in neurocomputational modeling [12
], as follows:
Where Φ is the field potential in volts, s is conductivity of the medium surrounding the neurons, in 1/Ωm, Imi is the transmembrane current in amperes across the ith neural compartment, and Ri is the distance from the ith neural compartment to the recording electrode. The sum is taken over every computational compartment of every neuron in the simulation.
The simulation was implemented using the GENESIS neural modeling language [12
] and run under LINUX on a dual-processor PC. The individual pyramidal cell and interneuron models were ported to GENESIS by Sampat and Huerta (http://www.genesis-sim.org/BABEL/babeldirs/cells
) and Menschik and Finkel [61
], respectively; both are available at the aforementioned URL. Programs to specify cell placement and connectivity were written in the C programming language by the author. These were designed to be maximally flexible: parameter files to the programs, in the form of and , can be updated as additional neuroanatomic information becomes available.
We used the Crank-Nicholson numerical integration method [19
], with an integration time step of 0.125 msec. This was among the fastest numerical integration methods available using the GENESIS software and, with the time step used, did not produce significant error [12
]. A simulation of 2.5 seconds of “neuronal time” required about 24 hours of computer time. For questions on programming details, readers are invited to email the author at psiekmeier/at/mclean.harvard.edu