An Indexing Theory for Working Memory Based on Fast Hebbian Plasticity

Neuron model

We use an integrate-and-fire point neuron model with spike–frequency adaptation (Brette and Gerstner, 2005), which was modified by Tully et al. (2014) for compatibility with a custom-made Bayesian Confidence Propagation Neural Network (BCPNN) synapse model in NEST (see Simulation environment) through the addition of the intrinsic excitability current . The model was simplified by excluding the subthreshold adaptation dynamics. Membrane potential (V_m ) and adaptation current are described by the following equations: (1) (2)

The membrane voltage changes through incoming currents over the membrane capacitance (C_m ). A leak reversal potential (E_L ) drives a leak current through the conductance (g_L ), and an upstroke slope factor (Δ _T ) determines the sharpness of the spike threshold (V_t ). Spikes are followed by a reset of membrane potential to V_r . Each spike increments the adaptation current by b, which decays with time constant . Simulated basket cells feature neither the intrinsic excitability current nor this spike-triggered adaptation.

In addition to external input I_ext (see Stimulation protocol), neurons receive a number of different synaptic currents from their presynaptic neurons in the network (AMPA, NMDA, and GABA), which are summed at the membrane accordingly: (3)

Synapse model

Excitatory AMPA and NMDA synapses have a reversal potential E^AMPA = E^NMDA , while inhibitory synapses drive the membrane potential toward E^GABA. Every presynaptic input spike (at with transmission delay t_ij ) evokes a transient synaptic current through a change in synaptic conductance that follows an exponential decay with time constants τ^syn depending on the synapse type (τ^AMPA ≪ τ^NMBA), as follows: (4)

H(·) is the Heaviside step function. is the peak amplitude of the conductance transient, learned by the spike-based BCPNN learning rule (next section). Plastic synapses are also subject to synaptic depression (vesicle depletion) according to the Tsodyks–Markram formalism (Tsodyks and Markram, 1997), modeling the transmission-dependent depletion of available synaptic resources by a utilization factor U, and a depression/reuptake time constant τ_rec, as follows: (5)

Spike-based BCPNN learning rule

Plastic AMPA and NMDA synapses are modeled to mimic NMDA-dependent Hebbian short-term potentiation (Erickson et al., 2010) with a spike-based version of the BCPNN learning rule (Wahlgren and Lansner, 2001; Tully et al., 2014). For a full derivation from Bayes rule, deeper biological motivation, and proof of concept, see Tully et al. (2014) and an earlier STM model implementation by Fiebig and Lansner (2017).

Briefly, the BCPNN learning rule makes use of biophysically plausible local traces to estimate normalized presynaptic and postsynaptic firing rates, as well as coactivation, which can be combined to implement Bayesian inference because connection strengths and neural unit activations have a statistical interpretation (Sandberg et al., 2002; Fiebig and Lansner, 2014; Tully et al., 2014). Crucial parameters include the synaptic activation trace Z, which is computed from spike trains via presynaptic and postsynaptic time constants , which are the same here but differ between AMPA and NMDA synapses, as follows: (6)

The larger NMDA time constant reflects the slower closing dynamics of NMDA receptor-gated channels. All excitatory connections are drawn as AMPA and NMDA pairs, such that they feature both components. Further filtering of the Z traces leads to rapidly expressing memory traces (referred to as P-traces) that estimate activation and coactivation as follows: (7)

These traces constitute memory itself and decay in a palimpsest fashion. Short-term potentiation decay is known to take place on timescales that are highly variable and activity dependent (Volianskis et al., 2015; see Discussion, The case for Hebbian plasticity).

We make use of the learning rule parameter κ (Eq. 7), which may reflect the action of endogenous neuromodulators [e.g., dopamine (DA) acting on D₁ receptors (D1Rs)] that signal relevance and thus modulate learning efficacy). It can be dynamically modulated to switch off learning to fixate the network or temporarily increase plasticity (κ_encoding, κ_normal; Table 1). In particular, we trigger a transient increase of plasticity concurrent with external stimulation.

Tully et al. (2014) showed that Bayesian inference can be recast and implemented in a network using the spike-based BCPNN learning rule. Prior activation levels are realized as an intrinsic excitability of each postsynaptic neuron, which is derived from the postsynaptic firing rate estimate p_j and implemented in the NEST neural simulator (Gewaltig and Diesmann, 2007) as an individual neural current with scaling constant β_gain: (8)

is thus an activity-dependent intrinsic membrane current to the neurons, similar to the A-type potassium channel (Hoffman et al., 1997) or TRP channel (Petersson et al., 2011). Synaptic weights are modeled as peak amplitudes of the conductance transient (Eq. 4) and determined from the logarithmic BCPNN weight, as derived from the P-traces with a synaptic scaling constant , as follows: (9)

In this model, AMPA and NMDA synapses make use of and , respectively. The logarithm in Equations 8 and 9 is motivated by the Bayesian underpinnings of the learning rule and means that synaptic weights multiplex both the learning of excitatory and disynaptic inhibitory interaction. The positive weight component is here interpreted as the conductance of a monosynaptic excitatory pyramidal to pyramidal synapse [Fig. 1, plastic connection to the coactivated minicolumn (MC)], while the negative component (Fig. 1, plastic connection to the competing MC) is interpreted as disynaptic via a dendritic targeting and vertically projecting inhibitory interneuron like a double bouquet and/or bipolar cell (Tucker and Katz, 2003; Kapfer et al., 2007; Ren et al., 2007; Silberberg and Markram, 2007). Accordingly, BCPNN connections with a negative weight use a GABAergic reversal potential instead, as in previously published models of this kind (Tully et al., 2014, 2016; Fiebig and Lansner, 2017). Model networks with negative synaptic weights have been shown to be functionally equivalent to those with both excitatory and inhibitory neurons with only positive weights (Parisien et al., 2008). In the context of this particular model microcircuit and learning rule, this was explicitly and conclusively demonstrated by the addition of double bouquet cells (Chrysanthidis et al., 2019).

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Local columnar connectivity within STM and LTM. Connection probabilities are given by the percentages; further details are in Tables 1, 2, and 3. The strength of plastic connections develops according to the synaptic learning rule described in the spike-based BCPNN learning rule. Initial weights are low and distributed by a noise-based initialization procedure (see Stimulation protocol). However, dashed connections are not plastic in LTM (besides the synaptic depression of Eq. 5), but already encode memory patterns previously learned through an LTP protocol, and loaded before the simulation using receptor-specific weights found in Table 2.

Code for the NEST implementation of the BCPNN synapse is openly available (see Code accessibility).

Axonal conduction delays

We compute axonal delays t_ij between presynaptic neuron i and postsynaptic neuron j, based on a constant conduction velocity V and the Euclidean distance between respective columns. Conduction delays were randomly drawn from a normal distribution with mean according to the connection distance divided by conduction speed and with a relative SD of 15% of the mean in order to account for individual arborization differences and varying conduction speeds as a result of axonal thickness/myelination. Further, we add a minimal conduction delay of 1.5 ms to reflect not directly modeled delays, such as diffusion of transmitter over the synaptic cleft, dendritic branching, thickness of the cortical sheet, and the spatial extent of columns, as follows: (10)

STM network architecture

The model organizes cells in the three simulated cortical areas into grids of nested hypercolumns (HCs) and MCs, sometimes referred to as macro columns, and “functional columns,” respectively. The STM network is simulated with HCs spread out on a grid with spatial extent of 17 × 17 mm. This spatially distributed network of columns has sizable conduction delays due to the distance between columns and can be interpreted as a spatially distributed subsampling of columns from the extent of dorsolateral PFC (e.g., BA 46 and 9/46, which also have a combined spatial extent of ∼289 mm² in macaque).

Each of the nonoverlapping HCs has a diameter of ∼640 μm, comparable to estimates of cortical column size (Mountcastle, 1997), contains 48 basket cells, and its pyramidal cell population has been divided into 12 MCs. This constitutes another subsampling from the ∼100 MCs per HC when mapping the model to biological cortex. We simulate 20 pyramidal neurons per MC to represent approximately the layer 2 population of an MC, 5 cells for layer 3A, 5 cells for layer 3B, and another 30 pyramidal cells for layer 4, as macaque BA 46 and 9/46 have a well developed granular layer (Petrides and Pandya, 1999). The STM model thus contains ∼18,000 simulated pyramidal cells in four layers (although layers 2, 3A, and 3B are often treated as one layer 2/3).

STM network connectivity

The most relevant connectivity parameters are found in Tables 1, 2, and 3. Pyramidal cells project laterally to basket cells within their own HC via AMPA-mediated excitatory projections with a connection probability of p_{p – B} (i.e., connections are randomly drawn without duplicates until the target fraction of all possible pre–post connections exist). In turn, they receive GABAergic feedback (FB) inhibition from basket cells (p_{B – p} ) that connect via static inhibitory synapses rather than plastic BCPNN synapses. This strong loop implements a competitive soft WTA (winner-take-all) subnetwork within each HC (Douglas and Martin, 2004). Local basket cells fire in rapid bursts, and induce alpha/beta oscillations in the absence of attractor activity and gamma, when attractors are present and active.

Pyramidal cells in layer 2/3 form connections both within and across HCs at connection probability p_{L 23 e - L 23 e} . These projections are implemented with plastic synapses and contain both AMPA and NMDA components, as explained in the subsection Spike-based BCPNN learning rule. Connections across columns and areas may feature sizable conduction delays due to the implied spatial distance between them (Table 1).

Pyramidal cells in layer 4 project to pyramidal cells of layer 2/3, targeting 25% of cells within their respective MC only. Experimental characterization of excitatory connections from layer 4 to layer 2/3 pyramidal cells have confirmed similarly high fine-scale specificity in rodent cortex (Yoshimura and Callaway, 2005) and, in turn, full-scale cortical simulation models without functional columns have found it necessary to specifically strengthen these connections to achieve defensible firing rates (Potjans and Diesmann, 2014).

In summary, the STM model thus features a total of 16.2 million plastic AMPA- and NMDA-mediated connections between its 18,000 simulated pyramidal cells, as well as 67,500 static connections from 9000 layer four pyramidals to layer 2/3 targets within their respective MC, and 1.2 million static connections to and from 1200 simulated basket cells.

LTM network

We simulate two structurally identical LTM networks, referred to as LTMa and LTMb. LTM networks may be interpreted as a spatially distributed subsampling of columns from areas of the parietotemporal cortex commonly associated with modal LTM stores. For example, inferior temporal cortex (ITC) is often referred to as the storehouse of visual LTM (Miyashita, 1993). Two such LTM areas are indicated in Figure 2.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Schematic of modeled connectivity within and across representative STM and LTM areas in macaque. STM features 25 HCs, whereas LTMa and LTMb both contain 16 simulated HCs. Each network spans several hundred square millimeters, and the simulated columns constitute a spatially distributed subsample of biological cortex, defined by conduction delays. Pyramidal cells in the simulated supragranular layers form connections both within and across columns. STM features an input layer 4 that shapes the input response of cortical columns, whereas LTM is instead stimulated directly to cue the activation of previously learned long-term memories. Additional corticocortical connections (feedforward in brown, feedback in dashed blue) are sparse (<1% connection probability) and implemented with terminal clusters (rightmost panels) and specific laminar connection profiles (bottom left). The connection schematic illustrates laminar connections realizing a direct supragranular forward-projection, as well as a common supragranular backprojection. Layer 2/3 recurrent connections in STM (dashed green) and corticocortical backprojections (dashed blue) feature fast Hebbian plasticity. For an in-depth model description, including the columnar microcircuits, please refer to Materials and Methods and Figure 1.

We simulate HCs in each area and 9 MCs per HC (Tables 1, 2, 3, for further details). Both LTM networks are structurally very similar to the previously described STM, yet they do not feature plasticity among their own cells, beyond short-term dynamics in the form of synaptic depression. Unlike STM, LTM areas also do not feature an input layer 4, but are instead stimulated directly to cue the activation of previously learned long-term memories (see Stimulation protocol). Various previous models with identical architecture have demonstrated how attractors can be learned via plastic BCPNN synapses (Lansner et al., 2013; Tully et al., 2014, 2016; Fiebig and Lansner, 2017). We load each LTM network with nine orthogonal attractors [see Fig. 4B, 10 in the example (which features two sets of five memories each)]. Each memory pattern consists of 16 active MCs, distributed across the 16 HCs of the network. We load in BCPNN weights from a previously trained network (Table 2), but thereafter set κ = 0 to deactivate plasticity of recurrent connections in LTM stores.

In summary, the two LTM models thus feature a total of 7.46 million connections between 8640 pyramidal cells, as well as 435,456 static connections to and from 1152 basket cells.

Interarea connectivity

In this model, we focus on layers 2/3, as its high degree of recurrent connectivity (Thomson et al., 2002; Yoshimura and Callaway, 2005) supports attractor function. The high fine-scale specificity of dense stellate cell (Yoshimura et al., 2005) and double-bouquet cell inputs (DeFelipe et al., 2006; Chrysanthidis et al., 2019) enable strongly coding subpopulations in the superior layers of functional columns. This fits with the general observation that layers 2/3 are more input selective than the lower layers (Crochet and Petersen, 2009; Sakata and Harris, 2009) and thus of more immediate concern to our computational model.

The recent characterization of supragranular feedforward (FF) and FB projections (from large cells in layer 3B and 3A, respectively), between association cortices and at short and medium cortical distances (Markov et al., 2014), allows for the construction of a basic cortical hierarchy without explicit representation of infragranular layers (and its long-range FB projections from large cells in layer 5 and 6). This is not to say that nothing would be gained by explicitly modeling infragranular layers, but it would go beyond the scope of this model.

Accordingly, our model implements supragranular FF and FB pathways between cortical areas that are at a medium distance in the cortical hierarchy. The approximate cortical distance between ITC and dlPFC in macaque is ∼40 mm and with an axonal conductance speed of 2 m/s, distributed conduction delays in our model (Eq. 10) average just >20 ms between these areas (Girard et al., 2001; Thorpe and Fabre-Thorpe, 2001; Caminiti et al., 2013).

In the forward path, layer 3B cells in LTM project toward STM (Fig. 2). We do not draw these connections one by one, but as branching axons targeting 25% of the pyramidal cells in a randomly chosen MC (the chance of any layer 3B cell to target any MC in STM is only 0.15%). The resulting split between targets in layer 2/3 and 4 is typical for FF connections at medium distances in the cortical hierarchy (Markov et al., 2014) and has important functional implications for the model (LTM-to-STM forward dynamics). We also branch off some inhibitory corticocortical connections as follows: for every excitatory connection within the selected targeted MC, an inhibitory connection is created from the same pyramidal layer 3B source cell onto a randomly selected cell outside the targeted MC, but inside the local HC. This way of drawing random forward-projections retains a degree of functional specificity due to its spatial clustering and yields patchy sparse forward-projections as observed in the cortex (Houzel et al., 1994; Voges et al., 2010), with a resulting interarea connection probability of only 0.0125% (648 axonal projections from L3B cells to STM layers 2/3 and 4 results in ∼20,000 total connections after branching, as described above.

In the FB path, we draw sparse plastic connections from layer 3A cells in STM to layer 2/3 cells in LTM: branching axons target 25% of the pyramidal cells in a randomly chosen HC in LTM, simulating a degree of axonal branching found in the literature (Zufferey et al., 1999). Using this method, we obtain biologically plausible sparse and structured FB projections with an interarea connection probability of 0.66%, which, unlike the forward pathway, do not have any built-in MC specificity but may develop such through activity-dependent plasticity. More parameters on corticocortical projections can be found in Table 3. On average, each LTM pyramidal cell receives ∼120 corticocortical connections from STM. Because ∼5% of STM cells fire together during memory reactivation (see Results), this means that a mere 6 active synapses per target cell are sufficient for driving (and thus maintaining) LTM activity from STM (there are 96 active synapses from coactive pyramidal cells in LTM).

Notably LTMa and LTMb have no direct pathways connecting them in our model since we assume that the use of previously not associated stimuli in our simulated multimodal tasks and further, that plasticity of biological connections between them are likely too slow (LTP timescale) to make a difference in WM dynamics. This arrangement also guarantees that any binding of long-term memories across LTM areas must be the result of interaction via STM instead. Overall in our model, corticocortical connectivity is very sparse, <1% on a cell-to-cell basis.

Stimulation protocol

The term I_ext in Equation 1 subsumes specific and unspecific external inputs. To simulate unspecific input from nonsimulated columns, and other areas, pyramidal cells are continually stimulated with a zero mean noise background throughout the simulation. In each layer, two independent Poisson sources generate spikes at rate and connect onto all pyramidal neurons in that layer, via nondepressing conductances (Table 2). Before each simulation, we distribute the initial values of all plastic weights by a process of learning from 1.5 s low, unstructured background activity (Table 2; ). To cue the activation of a specific memory pattern (i.e., attractor), we excite LTM pyramidal cells belonging to a memory patterns component MC with an additional excitatory Poisson spike train (rate, r_cue, length, t_cue; conductance, g_cue). As LTM patterns are strongly encoded in each LTM, a brief 50 ms stimulus is usually sufficient to activate any given memory.

Synthetic field potentials and spectral analysis

We estimate local field potentials (LFPs) by calculating a temporal derivative of the average low-pass filtered (cutoff frequency at 250 Hz) potential for all pyramidal cells in local populations at every time step, similarly to the approach adopted by Ursino and La Cara (2006). Although LFP is more directly linked to the synaptic activity (Logothetis, 2003), the averaged membrane potentials have been reported to be correlated with LFPs (Okun et al., 2010). In particular, low pass-filtered components of synaptic currents reflected in differentiated membrane potentials appear to carry the portion of the power spectral content of extracellular potentials that is relevant to our key findings (Lindén et al., 2010). As regards the phase response of estimated extracellular potentials, the delays of different frequency components are spatially dependent (Lindén et al., 2010). However, irrespective of the LFP synthesis, the phase-related phenomena reported in this study remain qualitatively unaffected since they hinge on relative rather than absolute phase values.

Most spectral analyses have been conducted on the synthesized field potentials with the exception of population firing rates, shown in Figure 3, A and B. Spectral information is extracted with a multitaper approach using a family of orthogonal tapers produced by Slepian functions (Slepian, 1978; Thomson, 1982), with frequency-dependent window lengths corresponding to five to eight oscillatory cycles and frequency smoothing corresponding to 0.3–0.4 of the central frequency, which was sampled with the resolution of 1 Hz (this configuration implies that two to three tapers are usually used). To obtain the spectral density, spectrotemporal content is averaged within a specific time interval.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

Basic network behavior in spike rasters and population firing rates. Activity in the untrained network under strong background input. A, Subsampled spike raster of STM (top) and LTM (bottom) layer 2/3 activity. HCs are separated by gray horizontal lines. Global oscillations in the alpha range (10–13 Hz) characterize this activity state in both STM (top) and LTM (bottom) in the absence of attractors. Inset, Power spectral density of LFP of each network. B, Cued LTM memory activation express as fast oscillation bursts of selective cells (50–80 Hz), organized into a theta-like envelope (4–8 Hz), see also power spectrum inset. The gamma band is broad due to the varying lengths of the underlying cycles (i.e., noticeably increasing over the short memory activation period). The underlying spike raster shows layer 2/3 activity of the activated MC in each HC, revealing spatial synchronization. The brief stimulus is a memory-specific cue. C, LTM-to-STM forward dynamics as shown in population firing rates of STM and LTM activity following LTM activation induced by a 50 ms targeted stimulus at time 0. LTM-driven activations of STM are characterized by an FF delay. Shadows indicate the SD of 100 peristimulus activations in LTM (blue) and STM (orange) with and without plasticity enabled (dashed, dark orange). Horizontal bars indicate the activation half-width (Materials and Methods). Onset is denoted by vertical dashed lines. The stimulation of LTM and the activation of plasticity is denoted underneath. D, Subsampled spike raster of STM (top) and LTM (middle) during forward activation of the untrained STM by five different LTM memory patterns, triggered via specific memory cues in LTM at times marked by the vertical dashed lines. Bottom spike raster shows LTM layer 2/3 activity of one selective MC per activated pattern (colors indicate different patterns). Top spike raster shows layer 2/3 activity of one HC in STM. STM spikes are colored according to each cells dominant pattern selectivity (based on the memory pattern correlation of individual STM cell spiking during initial pattern activation, see Materials and Methods, Spike train analysis and memory activity tracking). Bottom, The five stimuli to LTM (colored boxes) and modulation of STM plasticity (black line). Extended Data Figure 3-1 shows basic network behavior in spike rasters and population firing rates under low-input feature fluctuations in membrane voltages and low-rate, asynchronous spiking activity, while Extended Data Figure 3-2 shows network activity during plasticity-modulated stimulation with 20% spatial extent, illustrating the impact of conductance delays on cortical dynamics (see Model robustness).

The coherence for a pair of synthesized field potentials at the spatial resolution corresponding to a hypercolumn was calculated using the multitaper auto-spectral and cross-spectral estimates. The complex value of coherence (Carter, 1987) was evaluated first based on the spectral components averaged within 0.5 s windows. Next, its magnitude was extracted to produce the time-windowed estimate of the coherence amplitude. In addition, phase-locking statistics were estimated to examine synchrony without the interference of amplitude correlations (Lachaux et al., 1999; Palva et al., 2005). In particular, the phase-locking value (PLV) between two signals with instantaneous phases Φ₁(t) and Φ₂(t) was evaluated within a time window of size N = 0.5 s as follows: The instantaneous phase of the signals was estimated from their analytic signal representation obtained using a Hilbert transform. Before the transform was applied, the signals were narrow band filtered with low-time domain spread, finite-impulse response filters (in the forward and reverse directions to avoid any phase distortions). The analysis was performed mainly for gamma-range oscillations. A continuous PLV estimate was obtained with a sliding window approach, and the average along with SE were calculated typically over 25 trials.

Spike train analysis and memory activity tracking

We track memory activity in time by analyzing the population firing rate of pattern-specific and network-wide spiking activity usually using an exponential moving average filter time constant of 20 ms. We do not use an otherwise common low-pass filter with symmetrical window, because we are particularly interested in characterizing activation onsets and onset delays. As activations are characterized by sizable gamma-like bursts, a simple threshold detector can extract candidate activation events and decode the activated memory. This is trivial in LTM due to the known nature of its patterns. In STM, we decode the stimulus specificity of each cell individually by finding the maximum correlation between input pattern and the untrained STM spiking response in the 320 ms following cue onset (which is the stimulation interval during the plasticity-modulated stimulation period; Fig. 3D) following the pattern cue to LTM. Thereafter, we can filter the population response of cells in STM with the same selectivity on that basis to obtain a more robust readout. We validate the specificity by means of cross-correlations, which reveal that the pattern-specific populations are rather orthogonal according to the covariance matrix (off-diagonal magnitude, <0.1). In all three networks, we measure the onset and offset of pattern activity by thresholding each individual activation at half of its population peak firing rate. In LTM, we further check pattern completion by analyzing component MC activation. Whenever targeted stimuli are used, we analyze peristimulus activation traces. When activation onsets are less predictable, such as during free STM-paced maintenance, we extract activation candidates via a threshold detector trained at the 50th percentile of the cumulative distribution of the population firing rate signal.

Code accessibility

We used the NEST simulator (Gewaltig and Diesmann, 2007) version 2.2 for our simulations (RRID:SCR_002963), running on a Cray XC-40 Supercomputer of the PDC Centre for High Performance Computing. The custom-built spiking neural network implementation of the BCPNN learning rule for MPI (message passing interface) parallelized NEST is freely available on github (https://github.com/Florian-Fiebig/BCPNN-for-NEST222-MPI) and is included in the Extended Data 1. Further, the model is also available on ModelDB (https://modeldb.yale.edu/257610).

Model robustness

Our model incorporates a plethora of biological constraints, such as estimates of the extent and distance of areas (e.g., STM patch size approximates macaque dlPFC and is 40 mm from either LTM patch), laminar cell distributions ( , ,…), and hypercolumnar size. The model also abides by various electrophysiological constraints, such as plausible EPSP, IPSP sizes, estimates on laminar connection densities, laminar characterization of cortical FF/FB pathways with remote patchy connectivity, estimates on axonal conductance speeds, dendritic arbor sizes (branching factors), commonly accepted synaptic time constants for various receptor types, depression, adaptation, and builds on top of established models, such as the neuron model or the synaptic resource model. References to many of these constraints can be found throughout the Materials and Methods.

Because our model is quite complex and synthesizes many different components and processes, it is beyond the scope of this work to perform a detailed parameter sensitivity analysis. However, from our extensive simulations we conclude that it is robust and degrades gracefully. Almost all uncertain parameters can be varied ±30% without breaking WM function. The model is dramatically subsampled, and scaling up would be possible. This could be expected to further improve overall robustness. Highly related modular cortical network models have been studied extensively previously (Lundqvist et al., 2010, 2011; Tully et al., 2013, 2014; Fiebig and Lansner, 2017). For example, the model sensitivity to important short-term plasticity parameters affecting active maintenance mechanisms and intermittent gamma bursts (e.g., neural adaptation and synaptic depression time constants) were specifically explored in a single-network model (Fiebig and Lansner, 2017; see Fig. 8).

In the following, we briefly address new aspects of model sensitivity, previously unexplored, such as the parameterization of corticocortical connectivity, spatial scale (and associated conduction delays), as well as the transient modulation of Hebbian plasticity during rapid WM encoding.

In the FB pathway, a mere 0.6% connectivity is sufficient to support LTM activation in maintenance and recall. As rigorous testing (data not shown here) revealed, lower connectivity degrades WM capacity, unless we increase the total number of coactive STM cells by other means. FF connectivity can be even lower (0.015% in this model) because terminal clusters in STM are smaller and provide more information contrast (corticocortical connectivity). In both cases, our model uses very sparse connectivity, yet it could be increased or decreased if single synaptic currents were reduced/increased, respectively. Somewhat peculiarly, we also found that we needed to increase the corticocortical conductance of the backprojections ( ) by the same factor of 1.8 (over the local conductance gain ) as another highly detailed multiarea model of macaque visual cortex (Schmidt et al., 2018) to achieve functional WM at the stated long-distance connection probabilities.

There are upper and lower limits on conduction delays in our model. When corticocortical conduction delays exceed 65 ms (corresponding to 130 mm in distance), STM FB can no longer activate the LTM network because bursts desynchronize before they arrive. STM and LTM could be adjacent, as we briefly mention at the end of the Results section, but there is a minimum spatial scale for each component network. The length of gamma bursts decreases if we reduce the spatial extent (and thus the connection delays between HCs) by 45%. At 20%, when the largest inter-HC delays fall to <5 ms (Extended Data Fig. 3–2), the spiking activity of activated memories collapses into a single brief burst, which degrades learning and effective information transmission both within and across networks. Networks may be much smaller, however, if this is compensated by slower axonal conductance velocities (<2 mm/ms). Furthermore, we verified that the relative temporal delay dither in Equation 10 can be varied considerably (0–30%) without noticeable effects on memory performance.

The Hebbian plasticity of the model can be modulated via the parameter κ (Eq. 7). While κ is normally 1 (κ_normal, a transient increase of κ = κ_encoding; Table 1), it enables rapid, one-shot encoding in STM (Fig. 3D). Halving or doubling κ_encoding affects the overall working memory performance of the model only slightly, as measured by the number of items maintained during the delay period, or the overall rate of gamma bursts (Extended Data Fig. 4–4). It is, however, not possible to maintain normal WM operation without upregulating plasticity during encoding (leaving κ_encoding = κ_normal), unless additional compensatory changes are made to increase STM excitability, background excitation, or excitatory long-range connectivity. The strong correlation between working memory load and gamma-burst rate was previously discussed by Fiebig and Lansner (2017) in the context of evidence from multi-item WM recordings in macaque by Lundqvist et al. (2016).