1 Introduction

We examine and build on a cortical microcircuit model, previously used in a working memory model by Fiebig and Lansner (2017) that implemented a BCPNN (Bayesian Confidence Propagation Neural Network) learning rule. We then expand on this functional columnar architecture by integrating GABAergic double bouquet cells (DBCs), which may play a key modulatory role in the cortical microcircuit (Krimer et al. 2005; Kelsom and Lu 2013). Generally speaking, the BCPNN learning rule processes spike trains of pre- and postsynaptic neurons, and computes synaptic traces of activation and coactivation, which are then used to calculate the updated weights (see Section 2.3). In other words, the BCPNN is based on spike train correlations (Tully et al. 2014).

The foremost points of this work are first, the electrophysiological modeling of the inhibitory DBCs, and secondly their integration with the previous cortical memory model and its learning rule, which yields a novel model with improved biological plausibility and maintained functionality. The previous implementation suffers from the problem that learned weights among excitatory pyramidal cells in competing functional columns become negative (inhibitory), thus violating Dale’s Principle which states that neurons release the same neurotransmitters at all of their synapses (Strata and Harvey 1999). The biological plausibility of our model is here improved by the introduction of DBCs that provide disynaptic inhibition among pyramidal cells (Silberberg and Markram 2007).

2 Materials and methods

2.1 Neuron models

We use an AdEx IAF (Adaptive Exponential integrate-and-fire) neuron model with spike-frequency adaptation (Brette and Gerstner 2005). The neuron model has been modified for compatibility with a BCPNN synapse model (Tully et al. 2014) and reparameterized for simulation of several different neuron types. The model describes the temporal development of the membrane potential Vm and the adaptation current Iw, given by the following equations:

$$ C_{\mathrm{m}}\!\frac{dV_{\mathrm{m}}}{dt} = -g_{\mathrm{L}}\!(V_{\mathrm{m}} - E_{\mathrm{L}}) + g_{\mathrm{L}} {\varDelta}\tau e^{\frac{V_{\mathrm{m}}-V_{\mathrm{t}}}{{\varDelta}\tau}} - I_{\mathrm{w}} + I\!_{\text{ext}} + I\!_{\text{syn}} $$
(1)
$$ \frac{dI_{\mathrm{w}}}{dt}=\frac{-I_{\mathrm{w}}}{\tau_{\mathrm{w}}}+b\delta(t-t_{\text{sp}}) $$
(2)

Vm represents cell membrane potential, Iw stands for the adaptation current, Cm is the membrane capacitance, gL is the leak conductance, EL is the leak reversal potential, Vt is the spiking threshold, Δτ is the spike slope factor, b is the spike-triggered adaptation, τw is the adaptation recovery time constant and Iext is the stimulation current.

Adaptation enriches neural dynamics particularly in pyramidal cells (Brette and Gerstner 2005), and we take advantage of this in modeling DBCs as well. As in the model we extend, the neuron model was somewhat simplified by excluding the subthreshold adaptation. AdEx models provide a phenomenological description of the neural behaviour, yet feature limitations in predicting the precise time course of the voltage during and after a spike or the underlying biophysical causes of electrical activity (Gerstner and Naud 2009).

2.2 Synapse model

Besides the stimulation current Iext, neurons receive synaptic currents \(I_{\text {syn}_{j}}\) from AMPA and GABA synapses summed at the membrane. The model features conductance based AMPA (reversal potential EAMPA) and GABA (reversal potential EGABA) synapses:

$$ I_{\text{syn}_{j}} = \sum\limits_{syn}\sum\limits_{i}\!g_{ij}^{syn}\!(t)\!(V_{\mathrm{m}} - E_{ij}^{syn}) = I_{j}^{AMPA}(t) + I_{j}^{GABA}(t) $$
(3)

Plastic AMPA synapses under the spike-based BCPNN learning rule (see Section 2.3), are also subject to synaptic depression following the Tsodyks-Markram formalism (Tsodyks and Markram 1997) :

$$ \frac{dx_{ij}^{dep}}{dt} = \frac{1-x_{ij}^{dep}}{\tau_{rec}} - Ux_{ij}^{dep} \sum\limits_{sp}\delta(t-t_{sp}^{i}-t_{ij}) $$
(4)

2.3 BCPNN learning rule

Excitatory AMPA weights develop according to the BCPNN learning rule (Lansner and Ekeberg 1989; Sandberg et al. 2002; Tully et al. 2014). This is a Hebbian type of learning rule used in many previous works, most recently in Fiebig and Lansner (2017). It was derived from Bayes rule, assuming that a postsynaptic neuron employs some form of probabilistic inference to decide whether to emit a spike or not. It is a more complex learning rule than the more standard spike-timing-dependent plasticity (STDP) learning rule (Caporale and Dan 2008), but can replicate the main features of such plasticity. As other spiking synaptic learning rules, it is so far insufficiently validated against quantitative experimental data on biological synaptic plasticity.

A BCPNN synapse calculates three synaptic memory traces, Pi, Pj and Pij, implemented as exponentially weighted moving averages of pre-, post- and co-activation. As old memories deteriorate they are gradually replaced by newly learned patterns, so exponentially moving averages prioritize recent patterns. Specifically, BCPNN implements a three-stage procedure of exponential filters which defines Z, E and P traces. The method then estimates Pi (normalized presynaptic firing rate), Pj (normalized postsynaptic firing rate) and also Pij (coactivation) from these traces. In the final stage, Pi, Pj and Pij update the Bayesian weights wij and biases βj. It is worth adding that E traces that enable delayed reward learning, are not used here because such conditions are not applicable. Some of the key equations are highlighted in this chapter; yet for further information and deeper understanding of the BCPNN learning rule, see Tully et al. (2014).

To begin with, BCPNN receives pre- and postsynaptic spike trains (Si, Sj) so as to calculate the traces Zi and Zj:

$$ \left\{\begin{array}{l} \begin{array}{l} \tau_{\mathrm{z_{\mathrm{i}}}} \frac{dZ_{\mathrm{i}}}{dt}=\frac{S_{\mathrm{{i}}}}{f_{\max}t_{\text{spike}}}-Z_{\mathrm{i}}+\epsilon \\ \tau_{\mathrm{z_{\mathrm{j}}}} \frac{dZ_{\mathrm{j}}}{dt}=\frac{S_{\mathrm{{j}}}}{f_{\max}t_{\text{spike}}}-Z_{\mathrm{j}}+\epsilon \end{array} \end{array}\right. $$
(5)

\(f_{\max \nolimits }\) denotes the maximal neuronal spike rate, 𝜖 is the lowest attainable probability estimate, tspike denotes the spike duration while \(\tau _{\mathrm {z_{\mathrm {i}}}}=\tau _{\mathrm {z_{\mathrm {j}}}}\) are the pre- and postsynaptic time constants respectively (here 5 ms).

P traces then are estimated from the Z traces as follows:

$$ \begin{array}{@{}rcl@{}} \left\{\begin{array}{l} \begin{array}{l} \tau_{\mathrm{p}}\frac{dP_{\mathrm{i}}}{dt}=\kappa(Z_{\mathrm{i}}-P_{\mathrm{i}}) \\ \tau_{\mathrm{p}}\frac{dP_{\mathrm{j}}}{dt}=\kappa(Z_{\mathrm{j}}-P_{\mathrm{j}}) \\ \tau_{\mathrm{p}}\frac{dP_{\text{ij}}}{dt}=\kappa(Z_{\text{ij}}-P_{\text{ij}}) \end{array} \end{array}\right. \end{array} $$
(6)

The parameter κ adjusts learning speed, and by setting κ = 0 there are no weight changes. To give prominence to the stability of memory networks with BCPNN learning rule, we set κ = 1 during the whole simulation.

Finally, Pi, Pj and Pij are used to calculate intrinsic excitability βj and synaptic weights wij:

$$ \begin{array}{@{}rcl@{}} \left\{\begin{array}{l} \begin{array}{l} w_{\text{ij}}=\log\frac{P_{\text{ij}}}{P_{\mathrm{i}}P_{\mathrm{j}}}\\ \beta_{\mathrm{j}} = \log(P_{\mathrm{j}}) \end{array} \end{array}\right. \end{array} $$
(7)

2.4 Columnar network architecture

The proposed architecture principally follows several previous spiking neural network implementations (Fiebig and Lansner 2017; Lansner 2009) and is best understood as a subsampled cortical layer II/III model with nested hypercolumns (HCs) and minicolumns (MCs) (see Fig. 1a). This modular design attributes MCs as the core functional unit of cortex, composed of pyramidal cells with shared selectivity, forming a functional (not necessarily strictly anatomical) column. The high degree of recurrent connectivity within MCs (Thomson et al. 2002; Yoshimura and Callaway 2005) and horizontal connections between them link functional columns into larger attractors (Binzegger et al. 2009; Muir et al. 2011; Stettler et al. 2002). Densely recurrent non-specific feedback inhibition, here mediated by basket cells, implements a soft winner-take-all structure (Binzegger et al. 2009) amongst the functional columns. Recurrent excitatory gain can amplify and complete noisy inputs towards discrete embedded attractors. This approach does not address the role of infragranular layers and it does not apply directly to neural structures that do not follow the implied stereotypical columnar layout (such as hippocampus and rodent V1).

Fig. 1
figure 1

a The proposed functional columnar architecture including DBCs. The model contains 2 Hypercolumns (HC0, HC1) nesting 2 minicolumns each (MC0, MC1, MC2, MC3). Each MC contains one DBC, which delivers inhibition to all pyramidal cells in their MC. MCs are selective for different stimuli such that MC0 and MC2 function coactively and compete with another pair of coactive minicolumns MC1 and MC3. When MC0 and MC2 are active, MC1 and MC3 are silent and vice versa. Pyramidal neurons of MC0 and MC2 disynaptically inhibit MC1 and MC3 through basket cells (within HCs) and DBCs (within and between HCs). Plastic connections are drawn with 20% connection probability, while the connections between basket and pyramidal cells within the local hypercolumns are drawn with 70% connection probability. Conductance delays within and between HCs are 1.5 and 4.5 ms respectively. b Training stimuli drive pyramidal and DBCs. \(PYR_{MC0}^{HC0}\) denotes the pyramidal cells in MC0. A poisson generator (STIM0 - brown) stimulates \(PYR_{MC0}^{HC0}\), \(PYR_{MC2}^{HC1}\), \(DBC_{MC1}^{HC0}\) and \(DBC_{MC3}^{HC1}\) between 1000 ms and 2000 ms. Another poisson generator (STIM1 - purple) stimulates \(PYR_{MC1}^{HC0}\), \(PYR_{MC3}^{HC1}\), \(DBC_{MC0}^{HC0}\) and \(DBC_{MC2}^{HC1}\) between 3000 ms and 4000 ms. A zero mean noise poison generator (ZMN - green shaded area) is active throughout the simulation. c Membrane voltage of a stimulated DBC. STIM1 specifically drives this cell between 3000 ms and 4000 ms (cf. Fig. 1b). The DBC presents sustained low-rate firing throughout simulation and reaches typically reported firing rates during stimulation (see Table 1)

Full size image

2.5 Simulation tools

We use NEST (Neural Simulation Tool) version 2.4.2, and a custom-built BCPNN learning rule module (Tully et al. 2014). NEST simulates the dynamics of spiking neural models and features a convenient Python interface (PyNEST) to NEST’s simulation kernel (Gewaltig and Diesmann 2007). Simulation code and the custom module is available online in ModelDB (McDougal et al. 2017) at http://modeldb.yale.edu/257610.

3 Results

3.1 Double bouquet cells

Modeling DBCs is a key contribution of this work since the suggested cortical microcircuit model learns disynaptic inhibition through them and thereby modulates the neural activity of neurons in competing MCs. DBCs are GABAergic interneurons which may play an important role in shaping neural activity and circuitry (Krimer et al. 2005; Kelsom and Lu 2013) and are mainly located in Layer II/III featuring a bitufted dendritic conformation (Markram et al. 2004). They contact the dendrites of targeted cells (Markram et al. 2004) innervating spines (69.2% ± 4.2%) and shafts (30.8% ± 4.2%) (Tamas et al. 1997). DBCs are characterized by vertically oriented descending axons (María and DeFelipe 1995), which are generally termed “bundles” or “horse-tails” (Yáñez et al. 2005; DeFelipe et al. 1989).

The majority of double bouquet cells appear to be situated in upper layers (DeFelipe et al. 1989) and one of their unique feature is a horse-tail that fit well within the minicolumn (vertical cyclinder of tissue with a diameter of roughly 25 − 50 μm). Due to this morphology they create strong connections with pyramidal cells within their local column (DeFelipe et al. 2006). Neuroanatomical data suggests that each minicolumn contains one DBC (DeFelipe et al. 2006).

DBCs present a small degree of adaptation (Tamas et al. 1997) and unlike Basket Cells, they are characterized by sustained spiking activity (Krimer et al. 2005), which under strong stimulation is at 43 ± 13 Hz (Zaitsev et al. 2008). DBCs are intermediate spiking (IS) cells (Krimer et al. 2005) and subclassified as RSNP cells (Regular-spiking nonpyramical neurons) according to studies of their physiological properties (Kawaguchi and Kubota 1996; 1997).

Invitro experiments report DBC resting potential at − 76 ± 6 mV, spiking threshold Vth at − 44 ± 8 mV and input resistance at 626 ± 312 MΩ. They are considered inhibitory neurons with high input resistance (Zaitsev et al. 2008) and low capacitance (see Table 1). In addition, they exhibit membrane time constants (17.1 ± 7.7 ms), slope factor (− 0.64 ± 2.22) and action potential amplitude (66 ± 12.4 mV ) (Krimer et al. 2005).

Table 1 DBC model parameters
Full size table

We align the simulation model for DBCs with biological findings, yet tune some factors such as adaptation (b), leak conductance (gL), slope factor (Δt) and refractory period (τref) to achieve satisfactory electrophysiological fidelity, reproducing spike patterns under sweeps of increasing suprathreshold current steps and other typically reported activity. Figure 1c displays the membrane voltage of a stimulated DBC. The resulting model parameters are broadly consistent with experimentally reported values (see Table 1).

3.2 Adding DBCs to the columnar architecture

In the new model, connections among pyramidal cells in competing MCs are now mediated by DBCs as an additional local microcircuit component (see Fig. 1a). Their functional role is to deliver the same amount of inhibition to the respective MCs as the previous model but now entirely disynaptically, without principally changing established network learning and neural dynamics.

This extended cortical microcircuit model now contains three classes of neurons; pyramidal cells, basket cells and DBCs (see Fig. 1a). We use parameters for pyramidal and basket cells from a previous model implementation (Fiebig and Lansner 2017) and derive parameters for DBCs through electrophysiological modeling and tuning based on reported invitro characterizations (see Section 3.1). We simulate a small network of two reduced HCs from the previous larger network model.

3.3 BCPNN plasticity

Stimulation of the columnar network changes the efficacy of the plastic BCPNN synapses. To show the learned connectivity, we read out connection weights after an one second long initialization with zero mean noise (Initial Weight Distribution, IWD) and finally, after learning of the two stimuli (Learned Weight Distribution, LWD), see Fig. 1b.

Figure 2 shows histograms of plastic weights between the pyramidal cells in MC0 and their post-synaptic targets. Figure 2a and b show that the integrated DBCs in the columnar architecture do not affect the behaviour of the functional weights which are not related to them.

Fig. 2
figure 2

Weight distribution of plastic connections from pyramidal cells in MC0 before (Initial Weight Distribution, IWD) and after (Learned Weight Distribution, LWD) training in the previous and updated microcurcuit (multi-trial averages). IWD (new model) in blue, LWD (new model) in green, IWD (previous model) in orange and LWD (previous model) in purple. a Recurrent AMPA weight distribution among the pyramidal cells of MC0. Initial and converged distributions from both columnar networks are similar with means of 0 and 1.7 respectively, (see plot-box). b Distribution of AMPA weights between the pyramidal cells of the coactive MC0 and MC2 (associate connection between HCs). The distribution of the initial weights has the same behaviour in both architectures, with a mean close to 0. The distributions of the learned weights overlap, with a mean of 1.3. c The initial distribution of plastic weights in both functional columnar architectures is close to 0. The negative learned weights of the previous model (violating Dale’s principle) are functionally substituted in the new model by positive weights from the pyramidal cells of MC0 to the DBC of the competing MC1. The converged positive weight distribution has a mean of 0.5. d Plastic connections between pyramidal cells of MC0 and DBC of the competing MC3 are compared with plastic connections between pyramidal cells in MC0 and MC3 in the previous model. The behaviour remains the same as in C) because the negative weights turn positive after DBCs integration

Full size image

Two important changes can be identified in Fig. 2c and d, wherein disynaptic inhibition is introduced by DBCs. Even though the strength of the learned connections onto DBCs is weak, DBCs also feature low capacitance and dense connections with local pyramidal cells, and thus deliver comparable inhibition (see Section 3.4).

This result shows the effectiveness of DBCs involvement in the microcircuit network as they learn to mediate disynaptic inhibition between pyramidal cells in competing MCs. This outcome looks promising, but we yet have to verify its functional efficacy with regards to the total inhibition delivered (see Section 3.4).

3.4 Functionality verification

Figure 3a displays the spiking activity of neurons in a simulated HC (HC0). Although DBCs keep a low level sustained spiking activity throughout the simulation, they can reach higher firing rate during training (Zaitsev et al. 2008).

Fig. 3
figure 3

a Spike raster of neurons in HC0. Colors correspond to Fig. 1a, b. When STIM0 or STIM1 are ON (cf. Fig. 1b), pyramidal cells in MC0 excite basket cells within their local MC and DBC in competing MC which in turn inhibit pyramidal cells in MC1 and vice versa. b Population averaged total inhibitory input current received by pyramidal cells in MC0 in both architectures (100 trials, 10 ms bins). New model in purple, previous network in blue. DBC and basket cells simultaneously deliver inhibition to the neurons of MC0. The total inhibition current (IGABA) starts from zero level (2500 ms-3000 ms), then decreases (3000 ms-4000 ms) reaching a climax of 220 pA and finally stabilizes at zero (4000 ms-4500 ms) following the same pattern in the new and previous cortical model

Full size image

The cortical model learns as expected and the competing MCs inhibit each other by disynaptic inhibition mediated by DBCs and basket cells. But is this inhibition equivalent to the mono-synaptic inhibition learned by the previous model?

We tested learning in both the new and previous model using the same stimulation pattern and recorded the total inhibitory input current received by pyramidal cells in MC0 (see Fig. 3b). The proposed cortical model effectively delivers the same amount of disynaptic inhibition via basket cells and DBCs. The new model has thus improved biological credibility while maintaining the same functionality.

4 Discussion

This work aims at giving prominence to the double bouquet cells and their use as an integral part of a cortical microcircuit model. The population of DBCs is limited compared to other GABAergic neurons; however, they may play a key role in shaping neural activity. By integrating them into an established model, the new model now obeys Dale’s principle with maintained function. Indirectly, this result also verifies the biological plausibility of recent network models (Lansner 2009; Tully et al. 2016; Fiebig and Lansner 2017). The newly integrated DBCs effectively learn to mediate disynaptic inhibition between pyramidal cells thus eliminating negative learned weights between pyramidal cells which violate Dale’s principle.

In conclusion, the successful integration of an electrophysiological DBC model into an established cortical microcircuit design yields a novel functionally equivalent learning network with improved biological plausibility. This model suggests that DBCs have a quite well defined role in cortical memory networks.