Neurogrid simulates cortical cell-types, active dendrites, and top-down attention

A central challenge for systems neuroscience and artificial intelligence is to understand how cognitive behaviors arise from large, highly interconnected networks of neurons [1]. Bridging this gap in our understanding calls for large-scale yet biophysically realistic neural simulations to span levels of experimental investigation (biophysical, dendritic, neuronal, columnar, and area) [2]. Digital simulation bridges these levels, but the expense in time and electricity is great [3].

A Titan RTX GPU card consumes electricity at a rate of 250 W to simulate a cortex model with 4 million neurons connected by 24 billion synapses at a speed 510-fold slower than real-time (i.e. each biological second takes 510 s) [4]. These simulations are challenging because each neuron distributes its output to thousands of neurons and, in turn, aggregates inputs from thousands of neurons.

In the digital paradigm, replicating signals for distribution as well as summing signals for aggregation occurs entirely within computing cores (facilitated by programming constructs such as MapReduce). The communication network interconnecting these cores simply routes messages by relaying them across links that connect a core to its immediate neighbors. Sharing these core-to-core wires instead of routing dedicated point-to-point wires, as the cortex does, ameliorates the difficulty of mapping a three-dimensional brain onto a two-dimensional chip [5].

With shared wires communicating addresses that signal arrival of an action potential, or spike, at particular synapse, traffic scales as the total number of synaptic connections times the presynaptic population's average spike-rate [6]. This address-event bus has been successfully used to build networks with thousands of neurons with a few hundred synaptic connections each [7]. It has not scaled beyond millions of synaptic connections, the point at which traffic saturates the bus' signaling rate.

To break this communication bottleneck, Neurogrid adopts a hybrid analog–digital approach that follows design principles the neocortex uses to minimize its wiring [8]. Neurogrid uses fast digital routers, operating serially, to replicate signals for distribution and uses slow analog circuits, operating in parallel, to sum signals for aggregation. It follows the neocortex's design principles by performing both distribution and aggregation hierarchically. This neuromorphic architecture scales to billions of synaptic connections, between up to a million neurons.

Here, we demonstrate that Neurogrid simulates cortical models spanning five levels of experimental investigation: biophysical, dendritic, neuronal, columnar, and area (sections 5–7). Bridging these five levels with Neurogrid revealed a novel way active dendrites could mediate top-down attention. We begin with a brief review of Neurogrid's neuromorphic architecture (section 3) and simulation environment (section 4); detailed descriptions are available elsewhere [8].

2.1. Dimensionless neuron models

Neurogrid's models of soma, dendrite, gating variables and synapses are dimensionless. A leaky integrate-and-fire model of a neuron [9] is described by

$$\begin{equation}{C}_{\mathrm{s}}\frac{\mathrm{d}{V}_{\mathrm{m}}}{\mathrm{d}t}={G}_{\text{leak}}\left({E}_{\text{leak}}-{V}_{\mathrm{m}}\right)+{I}_{\mathrm{sin}}+{I}_{\text{Na}},\end{equation} \tag{ 1 }$$

where V_m is the soma membrane potential, C_s is the soma membrane capacitance, G_leak is the soma leak conductance and E_leak its reversal potential, I_sin is the injected current and I_Na is the voltage-dependent sodium current that generates spikes. Rewriting the above equation with E_leak as the reference voltage:

$$\begin{equation}{C}_{\mathrm{s}}\frac{\mathrm{d}{V}_{\mathrm{s}}}{\mathrm{d}t}=-{G}_{\text{leak}}{V}_{\mathrm{s}}+{I}_{\mathrm{sin}}+{I}_{\text{Na}},\end{equation} \tag{ 2 }$$

where V_s = V_m − E_leak. We modeled ${I}_{\text{Na}}\left({V}_{\mathrm{s}}\right)$ as

$$\begin{equation}{I}_{\text{Na}}\left({V}_{\mathrm{s}}\right)=\frac{1}{2}\frac{{V}_{\mathrm{s}}^{2}}{{V}_{\text{th}}}{G}_{\text{leak}},\end{equation} \tag{ 3 }$$

where V_th is the threshold voltage. By definition, it is the voltage at which the sodium conductance (dI_Na/dV_s) is equal to the leak conductance. At this voltage, the membrane potential stops decelerating and starts accelerating (inflection point), the onset of spiking. Substituting equation (3) into equation (2) and dividing by G_leak V_th gives the dimensionless model for a quadratic leaky integrate-and-fire (QLIF) neuron [10, 11]:

$$\begin{equation}{\tau }_{\mathrm{s}}\dot {{v}_{\mathrm{s}}}=-{v}_{\text{s}}+{i}_{\mathrm{sin}}+\frac{{v}_{\text{s}}^{2}}{2},\end{equation} \tag{ 4 }$$

where τ_s = C_s/g_leak (membrane time constant), v_s = V_s/V_th and ${i}_{\mathrm{sin}}={I}_{\mathrm{sin}}/\left({g}_{\text{leak}}{V}_{\text{th}}\right)$. Thus, models of the form equation (1) can be expressed as dimensionless models of the form equation (4) by changing the reference voltage to E_leak and normalizing all voltages, conductances and currents by V_th, G_leak and G_leak V_th, respectively.

The soma compartment's dimensionless model is

$$\begin{equation}\begin{aligned}{\tau }_{\text{s}}\dot {{v}_{\text{s}}}=& -{v}_{\text{s}}+{i}_{\mathrm{sin}}+\frac{{v}_{\text{s}}^{2}}{2}-{g}_{\text{K}}{v}_{\text{s}}-{g}_{\text{res}}{v}_{\text{s}}{p}_{\text{res}}(t)+{v}_{\mathrm{d}}+\sum {g}_{{\text{syn}}_{i}}\left({e}_{{\text{syn}}_{i}}-{v}_{\text{s}}\right)\end{aligned},\end{equation} \tag{ 5 }$$

where v_s is the membrane potential, τ_s the membrane time constant, i_sin the input current, and ${g}_{{\text{syn}}_{i}}$ and ${e}_{{\text{syn}}_{i}}$ the conductance and reversal potential, respectively, of synapse i. g_res is activated by p_res(t), a unit amplitude pulse that is active for a duration t_res after a spike (declared when v_s ⩾ 10), to model the neuron's refractory period. g_K models a high-threshold potassium conductance and is only activated when a spike occurs, decaying afterward. We modeled its dynamics as

$$\begin{equation}{\tau }_{\mathrm{K}}{\dot {g}}_{\mathrm{K}}=-{g}_{\mathrm{K}}+{p}_{\text{res}}(t){g}_{\mathrm{K}\infty },\end{equation} \tag{ 6 }$$

where τ_K is the decay time constant.

The dendrite compartment's dimensionless model is

$$\begin{equation}{\tau }_{\mathrm{d}}\dot {{v}_{\mathrm{d}}}=-{v}_{\mathrm{d}}+{i}_{\text{din}}+{i}_{\text{bp}}{p}_{\text{res}}(t)+\sum {\zeta }_{i}{g}_{{\text{syn}}_{i}}\left({e}_{{\text{syn}}_{i}}-{v}_{\mathrm{d}}\right)+\sum {g}_{{\text{ch}}_{i}}\left({e}_{{\text{ch}}_{i}}-{v}_{\mathrm{d}}\right),\end{equation} \tag{ 7 }$$

where v_d is the membrane potential, τ_d the membrane time constant, i_din the input current, ${g}_{{\text{syn}}_{i}}$, ${e}_{{\text{syn}}_{i}}$ and ζ_i the conductance, reversal potential and spatial decay factor, respectively, of synapse i, and ${g}_{{\text{ch}}_{i}}$ and ${e}_{{\text{ch}}_{i}}$ the conductance and reversal potential, respectively, of ion-channel population i. A current i_bp is injected for a duration t_res (same as in equation (5)) to model a backpropagating spike.

The decay factor, which models the spatial decay in dendritic trees, is given by [12]

$$\begin{equation}\zeta (n)=\frac{1}{4\sqrt{\pi }}{\left(1+{\left(\frac{1}{1-{\gamma }^{2}}\right)}^{\frac{1}{4}}\right)}^{2}\frac{{\gamma }^{n}}{\sqrt{n}},\end{equation} \tag{ 8 }$$

where γ is the silicon dendritic tree's decay constant and n the distance traveled in number of neurons.

The synaptic population's dimensionless model [13] is

$$\begin{equation}{\tau }_{\text{syn}}{\dot {g}}_{\text{syn}}=-{g}_{\text{syn}}+{p}_{\text{rise}}(t){g}_{\text{sat}},\end{equation} \tag{ 9 }$$

where g_syn is the instantaneous synaptic conductance, τ_syn the synaptic time constant and g_sat the maximum conductance. The unit amplitude pulse p_rise(t)'s width t_rise models the duration for which neurotransmitter is available in the cleft following a pre-synaptic spike.

An ion-channel population's gating variable is modeled as:

$$\begin{equation}{\tau }_{\mathrm{c}\mathrm{h}}\dot {c}=-c+{c}_{\mathrm{s}\mathrm{s}},\end{equation} \tag{ 10 }$$

where c_ss is its steady-state activation (or inactivation) and τ_ch its time constant. c_ss is given by

$$\begin{equation}{c}_{\text{ss}}=\frac{\alpha }{\alpha +\beta }\quad \mathrm{o}\mathrm{r}\quad \frac{\beta }{\alpha +\beta },\end{equation} \tag{ 11 }$$

where α and β model a channel's voltage-dependent opening and closing rates and are described by

$$\begin{equation}\begin{aligned}\alpha & =\frac{{v}_{\mathrm{d}}-{v}_{\text{th}}}{2}+\frac{\sqrt{{({v}_{\mathrm{d}}-{v}_{\text{th}})}^{2}+\frac{1}{4{s}^{2}}}}{2}\\ \beta & =-\frac{{v}_{\mathrm{d}}-{v}_{\text{th}}}{2}+\frac{\sqrt{{({v}_{\mathrm{d}}-{v}_{\text{th}})}^{2}+\frac{1}{4{s}^{2}}}}{2}.\end{aligned}\end{equation} \tag{ 12 }$$

Here, v_th is the membrane potential at which c_s = 1/2 and s is the slope at this point. α and β satisfy a difference relation, α − β = v_d − v_th, and a reciprocal relation, $\alpha \beta =1/\left(16{s}^{2}\right)$, resulting in a sigmoidal dependence of c_ss on v_d. The gating variable's time constant is given by

$$\begin{equation}\begin{aligned}{\tau }_{\text{ch}}& =\frac{{\tau }_{\mathrm{max}}-{\tau }_{\mathrm{min}}}{{\tau }_{\mathrm{max}}}{\tilde {\tau }}_{\text{ch}}+{\tau }_{\mathrm{min}}\enspace \\ {\tilde {\tau }}_{\text{ch}}& =\frac{1}{2s\left(\alpha +\beta \right)}{\tau }_{\mathrm{max}}\enspace \end{aligned}\end{equation} \tag{ 13 }$$

τ_ch is bell-shaped with a maximum value of τ_max when v_d = v_th and a minimum value of τ_min when |v_d − v_th| ≫ 1/(2s), to avoid unphysiologically short time constants.

An ion-channel population may have one or two gating variables. When using one gating variable, its associated maximum conductance g_max may be set by a synapse population's g_syn to model an ion-channel population that is voltage as well as ligand-gated (e.g., NMDA receptors). When using a pair of gating variables, the effective conductance is given by:

$$\begin{equation}{g}_{\text{ch}}=\frac{{g}_{\text{max0}}{c}_{0}{g}_{\text{max}1}{c}_{1}}{{g}_{\text{max0}}{c}_{0}+{g}_{\text{max}1}{c}_{1}},\end{equation} \tag{ 14 }$$

where g_max0 and g_max1 are the maximum conductances associated with gating variables c₀ and c₁, respectively. In this case, g_max0 and g_max1 may be set by two synapse populations' g_syn to model neuromodulation, and the gating variables' thresholds set low to eliminate voltage-dependence. Alternatively, g_max0 and g_max1 may be programmed to the same value, g_max, to model a channel that activates and inactivates. In this case, the above expression simplifies to ${g}_{\text{ch}}={g}_{\mathrm{max}}\frac{{c}_{0}{c}_{1}}{{c}_{0}+{c}_{1}}$. The pair of gating variables always share a programmable reversal potential e_ch.

The parameter values used to obtain the simulation results presented in figures 2–4 are listed in tables 1–4. Reversal potentials and gating-variable thresholds were converted to dimensionless form by assuming that the leak's reversal potential is −80 mV and the spike threshold is −50 mV (i.e., all voltages are normalized by V_th = 30 mV).

Table 1. Parameter values for models of cortical neuron types. Current injection for FS, RS and CH neuron types was modeled with i_sin and for IB neuron type with i_din. Parameter values were chosen to produce the best fit to the in vitro data. All time values are in ms.

									Channel
Soma					Dendrite			Conductance		Activation			Inactivation
τs	τK	gK∞	tres	isin	τd	ibp	idin	gmax	ech	τmax	vth	s	τmax	vth	s
FS	3	200	0.005	0.8	3.7–9.8	—	—	_	—	—	—	—	—	—	—	—
RS	15	200	50	0.1	0.08–1.42	—	—	—	—	—	—	—	—	—	—	—
IB	18	200	50	1	0	54	0	1.09–2.5	1	7.5	1	0.5	1.25	50	0.2	−0.5
CH	13	50	250	2	30–39	12	100	0	—	—	—	—	—	—	—	—

Table 2. Parameter values used to model NMDA spikes in basal dendrites. NMDA and AMPA reversal potentials as well as NMDA threshold were chosen to match physiological values [14]. Other parameter values were chosen to produce the best fit to the in vitro data. All time values are in ms.

Soma					Dendrite			NMDA						AMPA
τs	τK	gK∞	tres	isin	τd	idin	ibp	τsyn	trise	gsat	esyn	vth	s	τsyn	trise	gsat	esyn
20	—	—	1	0	30	0	0	150	4	500	2.7	2.3	1	7.25	0.6	25–240	2.7

Table 3. Parameter values used to model coincidence detection in apical dendrites. Soma compartment was modeled with RS neuron parameter values (see table 1). Sodium channel parameters are similar to Mainen et al [15], adjusted together with other parameter values to produce the best fit to the in vitro data. All time values are in ms.

							Sodium channel
Dendrite			Synapse				Conductance		Activation			Inactivation
τd	idin	ibp	τsyn	trise	gsat	esyn	gmax	ech	τmax	vth	s	τmax	vth	s
3	0	100	8	4	80	2.7	10	5	0.4	0.9	1.2	35	0.5	−2

Table 4. Parameter values for attentional modulation simulation. Excitatory and inhibitory neurons' soma compartments were modeled with RS and FS neuron parameter values, respectively (see table 1). NMDA, AMPA and dendrite were modeled with the parameter values for NMDA spikes in basal dendrites (see table 2). All time values are in ms.

FEF excitatory
Soma					AMPA					GABA
τs	τK	gK∞	tres	isin	τsyn	trise	gsat	esyn	γ	τsyn	trise	gsat	esyn	γ
15	200	50.0	0.1	0	7.25	0.6	8	2.7	0.65	15	1	0.01	0.1	0.97

FEF inhibitory
Soma					AMPA
τs	τK	gK∞	tres	isin	τsyn	trise	gsat	esyn	γ
3	200	0.005	0.8	0.15	7.25	0.6	10	2.7	0.6

V4 excitatory
Soma					Dendrite		AMPA					NMDA
τs	τK	gK∞	tres	isin	τd	idin	τsyn	trise	gsat	esyn	γ	τsyn	trise	gsat	esyn	vth	s
15	200	50	0.1	0	30	0.001	7.25	0.6	2.0	2.7	0.65	150	4	50	2.7	2.3	1

2.2. Programming Neurogrid

Neuronal and synaptic model parameters and connectivity are described in Python, similar to PyNN [16], and a C++ back-end uses this description to program Neurogrid's electronic circuit parameters and lookup tables. Examples of Python descriptions of soma, dendrite, synapse and gating variables are listed in table 5. A calibration procedure is used to obtain a set of conversion factors to map the neuronal and synaptic model parameters to the circuit parameters [17]. In a population of silicon neurons, these parameters are lognormally distributed with a coefficient of variation (CV) of about 10% (e.g., 1/τ_s has a CV of 7.2%). Lookup tables for intracolumn connections are stored in the Neuorocores' 256 × 16 bit RAM and for intercolumn connections in the daughterboard's 8 M × 32 bit RAM.

Table 5. Python constructs to program Neurogrid. The table lists example Python function calls for modeling a soma, dendrite, synapse and channel. All time values are in seconds.

Soma
Dendrite
Synapse
Channel

A neuron's axonal arbor distributes signals by hierarchically replicating them at branch points. Placing an axon's branch points as close as possible to its terminals replaces multiple axon segments with a single segment [18] (figure 1(a), top). This principle minimizes the axonal arbor's wiring.

Zoom In Zoom Out Reset image size

Moreover, extending several dendritic branches to meet a terminal branch of an axon allows that branch to make multiple synapses (bouton cluster) [19] (figure 1(b), top). Synaptic signals from many axons sum in a dendritic branch and these branches' signals aggregate hierarchically. This principle minimizes the dendritic tree's wiring [18].

Hierarchical distribution is emulated by interconnecting silicon-neuron arrays with a tree-like (rather than a mesh-like) network [20, 21]. That allows address-events to be replicated close to their destination and reduces traffic on the array-to-array links (figure 1(a), bottom). Emulating an axon terminal's bouton cluster allows a single delivered address-event to evoke postsynaptic potentials in several silicon neurons (figure 1(b), bottom). This reduces traffic further.

Hierarchical aggregation is emulated efficiently by modeling multiple overlapping dendritic trees with a single two-dimensional resistive grid. This shared analog circuit replicates the exponential spatial decay that transforms postsynaptic potentials into current delivered to a dendritic tree's trunk [22].

Emulating axonal arbors and dendritic trees enables Neurogrid to simulate cortical models scalably and efficiently. A cortical area is modeled by a group of Neurocores (figure 1(c)); Neurogrid has sixteen of these chips. Each cell layer (or cell type) is mapped onto a different Neurocore's two-dimensional silicon neuron array. Circular pools of neurons centered at the same (x, y) location on these Neurocores model a cortical column [24].

Intercolumn axonal projections are routed by using the presynaptic neuron's address to retrieve the target columns' centers (x, y) from an off-chip random-access memory (first distribution level). This memory is programmed to replicate the neocortex's function-specific intercolumn connectivity [25, 26].

Intracolumn axon collaterals are routed by copying a retrieved (x, y) address to all of a cortical area's Neurocores using the interchip tree network; unneeded copies are filtered using an on-chip memory (second distribution level). This memory is programmed to replicate the neocortex's stereotyped intracolumn connectivity [27].

Intralayer dendritic trees—arborizing over a circular disc centered on the cell body—are realized using the resistive grid mentioned earlier. Its space constant is adjusted electronically to match the arbor's radius; a transistor-based implementation makes this possible [28].

With intercolumn projections, intracolumn collaterals, and intralayer dendrites, five thousand synaptic connections may be made by: (1) routing an address-event to ten columns; (2) copying it to six layers in each column; and (3) evoking postsynaptic potentials in a hundred neighboring neurons in all but one layer (a 5.6 neuron arbor radius).

To demonstrate these routing mechanisms' functionality and scalability, we used Neurogrid to simulate a recurrent network of a million interneurons with billions of inhibitory synaptic connections (figure 1(d)). The interneurons were organized in sixteen 256 × 256-neuron arrays that formed a torus; the first and last layers were neighbors. Intracolumn collaterals routed an interneuron's spike to the three nearest layers on each side as well as back to its own layer. The inhibition evoked decreased exponentially with distance, due to spatial decay in the intralayer dendritic trees. Half of the inhibition a soma received came from 8000 neurons in a cylinder centered on it, 7 layers in height and 19 neurons in radius.

These local inhibitory interactions synchronized the interneurons globally, reproducing previous findings [29]. Spiking activity waxed and waned periodically (3.7 Hz), with individual interneurons skipping several cycles (0.42 spikes s⁻¹ mean). Nevertheless, their spikes were entrained to the global rhythm across all 16 layers, demonstrating the functionality and scalability of Neurogrid's routing mechanisms.

Software infrastructure facilitates performing large-scale neural simulations on Neurogrid. Firmware and drivers initialize and configure the hardware. Calibration procedures establish correspondence between parameters of neuronal models and their silicon analogs [17, 30]. A mapping tool translates multiscale model descriptions (written in Python) into hardware configurations (programmed in memory). A user interface interacts with simulations in real-time. And visualization widgets render activity or plot spikes of up to a million neurons in real-time. For the details of this software infrastructure, see section 2.

Neurogrid's neuron model has separate soma and dendrite compartments (figure 1(e)). The soma compartment expresses conductances that generate spikes and adapt their rate. The dendrite compartment expresses ligand and voltage-gated conductances (up to four types each). The programmed maximum synaptic or membrane conductance is scaled by a 'gating variable' [31]. A pair of gating variables may be used to model a population of ion-channels that are activated and inactivated by voltage (e.g., T-type Ca channels), activated by ligand and voltage (e.g., NMDA receptors), or modulated by a second ligand (e.g., dopamine). For these components' equations and parameters, see section 2.

To demonstrate the expressiveness of this neuron model, we simulated four prominent neocortical cell-types classified by Nowak et al [32], namely fast spiking (FS), regular spiking (RS), intrinsic bursting (IB), and chattering (CH).

FS neurons, normally associated with inhibitory interneurons, respond to a depolarizing current input with a high-frequency spike-train, showing little or no adaptation. RS neurons, generally associated with excitatory stellate cells, respond to a depolarizing current input with a low-frequency spike-train, showing adaptation. For both types of neurons, increasing the injected current increases the spike rate.

We modeled FS and RS neurons with the soma compartment alone, including adaptation for the latter (figures 2(a) and (b)). These models matched the FS neuron's high spike-rate and relatively constant inter-spike intervals; the RS neuron's low spike-rate and increasing inter-spike intervals; and both cell types' overall spike-rate increase with injected current.

Zoom In Zoom Out Reset image size

IB and CH neurons display more complex responses. IB neurons, generally associated with pyramidal cells, respond to a sufficiently large increase in injected current with a burst followed by single spikes, as observed in guinea-pig cingulate neocortex [33]. Bursting does not occur if the current is injected when the membrane is depolarized. CH neurons, also associated with pyramidal cells, exhibit fast rhythmic bursting, as observed in cat striate cortex [34]. Their burst frequency does not increase much with increasing injected current.

We modeled IB and CH neurons with the soma and dendrite compartments. Slow calcium channels that are deinactivated near the resting potential contribute to bursting [35]. Hence, we equipped the dendrite compartment with an activating as well as an inactivating gating variable (figure 2(c)). This model matched the onset of bursting in an IB neuron and the increase in tonic spike-rate with increasing injected current.

Rhythmic bursting could arise from bidirectional interactions between the soma and the dendrite [36]. Neurogrid's neuron model supports this by allowing spikes to back-propagate from the soma compartment to the dendrite compartment (purely passive in this case). In turn, the dendrite depolarizes the soma. These reciprocal interactions produce repetitive spiking [37]. A build-up of K⁺ conductance with each spike terminates the burst; this conductance's decay with time determines the inter-burst interval. This model matched the high spike-rate within a burst, the number of spikes per burst, and the moderate decrease in inter-burst interval with increasing injected current (figure 2(d)).

Several lines of experimental evidence point to a critical role of active dendritic properties in sensorimotor processing. NMDA spikes in basal dendrites [38] have been shown to sharpen sensory tuning of neocortical neurons [40]. Such spikes double an excitatory post-synaptic potential's (EPSP) amplitude when they occur in basal dendrites of neocortical layer 5 pyramidal neurons. Coincidence of EPSPs (pre-spikes) and backpropagating action potentials (post-spikes) in apical dendrites [39] has been shown to amplify correlated sensory and motor signals during active sensation [41]. Thus, active dendrites enhance sensorimotor processing.

To demonstrate that Neurogrid's neuron model could express active dendritic behavior, we replicated the NMDA spike's all-or-none behavior and the backpropagating spike's amplifying action. An NMDA spike occurs when a small increment of synaptic strength, above some threshold, more than doubles the amplitude of the EPSP, which then remains unchanged for further increments [38].

To replicate the NMDA spike, we expressed an NMDA synaptic conductance in a dendrite compartment as well as an AMPA conductance. When the AMPA conductance drove the dendrite's potential above the NMDA conductance's voltage-activation threshold, the postsynaptic potential doubled in amplitude, replicating the all-or-none behavior of an NMDA spike (figure 3(a)).

Zoom In Zoom Out Reset image size

To replicate the amplifying action of a backpropagating spike, we expressed sodium channels that activate and inactivate in a dendrite compartment as well as a synaptic conductance with time-constant between NMDA and AMPA conductances' to mimic a mix. When a backpropagating spike coincided with a postsynaptic potential, the sodium conductance activated and increased the dendrite's voltage further. The pre-spike's timing relative to the post-spikes' when voltage peaked as well as the voltage's asymmetric increase and decrease matched experimental observations in rat layer 5 pyramidal neurons [39] (figure 3(b)). Replicating these in vitro results did not require NMDA synapses' voltage-dependence, consistent with previous models [39].

To demonstrate that Neurogrid can help relate cognitive phenomena to biophysical mechanisms, an important goal of neuroscience, we replicated the effects of spatially selective top-down attention on the responses of visual cortical neurons recorded in awake, behaving macaques [42]. Amplification of postsynaptic potentials by NMDA or sodium conductances in apical dendrites of pyramidal cells could contribute to this multiplicative interaction. The former may underlie spike-rate changes associated with feedback-driven figure-ground modulations in macaque primary visual cortex [43]. This evidence led us to hypothesize that NMDA receptors could also account for the multiplicative changes in spike rate observed in visual cortical neurons during spatially selective attention [42].

Unlike most computational models, Neurogrid models are subjected to heterogeneity. Once calibrated, the median parameter value across a Neurocore's population of silicon neurons matches the specified model parameter value. Whereas the variance of parameter values is determined by the fabrication process, resulting in a spike-rate distribution that spans a decade. For comparison, firing rates vary by over three decades across a population of cortical neurons [44]. Hence, showing that a proposed mechanism is robust to heterogeneity increases its biological plausibility.

We tested robustness of NMDA-mediated multiplicative gain to heterogeneity by simulating top-down and bottom-up interactions between a visual cortical area (V4) and a frontal cortical area (the Frontal eye field, FEF) on Neurogrid (figure 4(a)). V4 was modeled with an excitatory neuronal population (128 × 128) while FEF was modeled with an excitatory and an inhibitory population (both 128 × 128). Besides V4's excitatory population, which had dendrite compartments, all others had just soma compartments.

Zoom In Zoom Out Reset image size

FEF activity, temporally sustained by local recurrent excitation and spatially constrained by local recurrent inhibition, represented the locus of attention in a 2D map of visual space, inspired by Ardid et al's (1D) attention model [45]. To explore NMDA receptors' role in attentional modulation, unlike in the previous model, columnar feedback projections from FEF to V4 modulated activity of V4 neurons exclusively through NMDA synapses.

We discovered that heterogeneity in AMPA conductances realizes gradual gain modulation from abrupt NMDA threshold crossings. This heterogeneity (CV of 26%) caused dendrite compartments to cross NMDA's voltage-activation threshold at different visual drives (figure 4(b)). This distributed threshold-crossing increased the V4 population's spike-rate gradually with visual drive; it would otherwise have increased abruptly. This increase was steeper with FEF feedback, matching multiplicative changes with attention observed in macaque V4 [42] (figure 4(c)).

Conductance-based synapses, active membrane conductances, multiple dendritic compartments, spike back-propagation, and cortical cell types have been emulated in neuromorphic chips [37, 46–53]. Our focus here was on the next step: deploying these in silico neuronal components in multiscale modeling. We simulated cortical models with up to 1 M neurons and 8G synaptic connections using 1800-fold less energy per synaptic activation than a GPU [4, 8] (120 pJ versus 210 nJ) ⁸ .

A hybrid analog–digital architecture made this possible. It implements hierarchical distribution and aggregation and it maps cortical circuitry columnarly. Following these principles the neocortex uses to minimize its wiring reduced traffic between cores greatly. This neuromorphic approach enabled Neurogrid to simulate multiscale neural models that integrate findings across five levels of experimental investigation energy- and time-efficiently. Importantly, Neurogrid achieves scale and efficiency without sacrificing biophysical detail, and thus truly supports multiscale modeling.

Neurogrid's neuronal model is sufficiently detailed to simulate various cortical cell types as well as active dendritic behavior and neuromodulatory effects. This combination of scale and detail enabled us to discover a biologically plausible gain mechanism for attentional modulation. This advance in simulation capacity could not only lead to a better understanding of neurological disorders such as ADHD, it could also help reverse-engineer the brain's hybrid analog–digital computational paradigm, which could lead to artificial intelligence systems with billions of units and trillions of parameters that consume watts instead of kilowatts.

BVB developed the synchrony, NMDA-spike, coincidence detection, and cortical neuron models. NAS developed the attention model. NNO and JJA developed calibration, mapping, interaction and visualization software. KB directed the research effort. BVB, NAS and KB wrote the manuscript.

KB and NNO are co-founders and equity owners of Femtosense Inc. The remaining authors declare that they have no competing financial interests.

We thank E Kauderer-Abrams for calibration algorithm development; B Softky and H S Seung for assistance in editing the manuscript; and K Chin for administrative support. This work was supported by an NIH Director's Pioneer Award (DPI-OD000965) and an NIH/NINDS Transformative Research Award (R01NS076460).

The data that support the findings of this study are available from the corresponding author upon reasonable request. Correspondence and requests should be addressed to KB (email: boahen@stanford.edu).