1 Introduction

Learning and memory are orchestrated by synaptic plasticity, the ability of synapses to change their efficacy in an activity-dependent manner. Donald O. Hebb’s postulate about how synaptic plasticity might manifest was paraphrased with the well known mantra: “cells that fire together, wire together” (Hebb, 1949; Shatz, 1992). The first proof of coincident pre- and postsynaptic population activity leading to potentiation (an increase in efficacy) came from pathway stimulation in hippocampal slices (Bliss and Lømo, 1973). It was later confirmed at the neuron pair level (Markram et al., 1997; Bi and Poo, 1998), and spike-time dependent plasticity (STDP) became a quintessential protocol to study Hebbian plasticity in vitro. In the early 2000s a plethora of cortical pathways were studied and plasticity proved to be synapse location- and therefore pathway-dependent (Sjöström and Häusser, 2006; Letzkus et al., 2006; Froemke et al., 2010). The molecular substrate of Hebbian coincidence detection is the N-methyl-D-aspartate (NMDA) receptor, which upon removal of the Mg²⁺ block by depolarization, conducts Ca²⁺ (Mayer et al., 1984). The calcium-control hypothesis, put forward by Lisman (1989) postulates that prolonged, moderate amounts of Ca²⁺ lead to depression (a decrease in efficacy) while large transients of Ca²⁺ lead to potentiation. By putting these together, it became evident that it is not necessarily the timing of the postsynaptic spike, but rather the depolarization induced entrance of calcium into the postsynaptic dendrite that is important to evoke changes in synaptic efficacy (Goldberg et al., 2002; Lisman and Spruston, 2005).

In parallel with slice electrophysiology, Hebbian plasticity was also studied through its effect on behavior via fear conditioning experiments (McKernan and Shinnick-Gallagher, 1997) and this line of research led to numerous new techniques for tagging and re-activating cells that participate in newly formed memories (Tonegawa et al., 2015). While these studies highlighted the need to study plasticity at the network level, most changes are expected to happen at the synapse level. Therefore, high-throughput methods tracking synaptic proteins like PSD95 (Ray et al., 2023) and α-amino-3-hydroxy-5-methyl-4-isoxazolepropionate (AMPA) subunit GluA1 (Graves et al., 2021; Kim et al., 2023) are currently being developed. While readily applicable to monitor synaptic efficacy in vivo, currently these techniques cannot be supplemented with recordings of neural activity, thus the reason for the changes in efficacy can only be speculated.

The bridge between in vitro pairs of neurons and in vivo behavior is often provided by complementary, simulation based-approaches. Early theoretical work explored the potential link between memories and cells that fire and therefore wire together, concentrating on the storage and retrieval of memories in strongly recurrent networks (Hopfield, 1982), which remained an active topic of research (Fusi and Abbott, 2007; Krotov and Hopfield, 2016; Widrich et al., 2020). Inspired by STDP experiments, modelers have postulated diverse plasticity rules that fit the most recent experimental findings (Gerstner et al., 1996; Kempter et al., 1999; Song et al., 2000; Pfister and Gerstner, 2006; Clopath et al., 2010). Models that take into account the crucial role of calcium in driving plasticity outcomes have also been proposed (Shouval et al., 2002; Graupner and Brunel, 2012; Rubin et al., 2005; Jȩdrzejewska-Szmek et al., 2017; Rodrigues et al., 2023). The calcium-based model of Graupner and Brunel (2012) describes the evolution of intracellular calcium concentration ([Ca²⁺]_i) given the pre- and postsynaptic spike trains and updates the efficacy of the synapse, upon [Ca²⁺]_i crossing thresholds for depression and potentiation. Bio-plausible plasticity rules have been shown to bring about Hebbian cell assemblies, i.e., groups of neurons that fire together to enable long-term memory storage (Litwin-Kumar and Doiron, 2014; Zenke et al., 2015; Fauth and Van Rossum, 2019; Kossio et al., 2021). A common theme in these models is the necessity for homeostatic plasticity to keep the networks stable. While experimental evidence exists for homeostatic plasticity (Turrigiano and Nelson, 2004), this has been shown to be too slow to induce stability (Zenke et al., 2017a) and alternative mechanism have been proposed (Vogels et al., 2011; Delattre et al., 2015). While these studies provide mechanistic explanations of learning and memory, they used point-neuron models, therefore neglecting the structural and functional importance of dendrites in plasticity (with the exception of Bono et al. (2017); Kastellakis and Poirazi (2019)). The compartmentalized nature of dendritic trees gives rise to spatial clustering of synapses (Kastellakis and Poirazi, 2019; Farinella et al., 2014; Iacaruso et al., 2017; Tazerart et al., 2020) and local, non-linear voltage events (Poirazi et al., 2003; Stuart and Spruston, 2015) both of which are thought to contribute to the removal of the Mg²⁺ block from NMDA receptors and therefore gate plasticity.

To investigate the effect of dendrites and multi-synaptic interactions for plasticity at the network level, in this study we equipped the biophysically detailed, large-scale cortical network model of Isbister et al. (2023) with our recently developed model of functional plasticity (Chindemi et al., 2022) between excitatory cells (Figure 1). At this level of detail, a calcium-based plasticity model is the natural choice since calcium dynamics are already part of the simulation (Chindemi et al., 2022). Moreover, it allowed us to model in vivo-like conditions, specifically the low extracellular calcium concentration ([Ca²⁺]_o), which has been experimentally shown to reduce plasticity (Inglebert et al., 2020) (Figure 1E). As we had access to the pre- and postsynaptic activity and efficacy of millions of synapses including their dendritic location we could characterize the rules governing plasticity at the microcircuit level. During simulated evoked activity, plasticity acted sparsely and specifically, although it was still able to reorganize the network’s dynamics, manifesting in more pattern-specificity after plasticity. Potentiation dominated in amplitude and depression counteracted it in frequency, which led to stable firing rates with-out explicitly introducing any homeostatic terms (Zenke et al., 2017a; Turrigiano and Nelson, 2004). We found plasticity was driven by the network topology beyond pairwise connectivity. Specifically, by co-firing cell assemblies, spatial clustering of synapses on dendrites, and the centrality of a connection in the network. Furthermore, we confirmed our testable predictions in the MICrONS (2021) dataset, an openly available electron microscopic reconstruction.

Figure 1:

2 Results

To study how plasticity and various network level features interact to explain structure-function relationships, we used a bio-realistic, large-scale cortical model of the rat non-barrel somatosensory cortex (nbS1). By doing so, we had a continuous readout of both the activity of all neurons and the efficacy of the synapses between them under in vivo-like, low [Ca²⁺]_o conditions (~ 1 mM instead of 2 2.5 mM in vitro). The current model improves on Markram et al. (2015) in terms of both anatomical, e.g., atlas-based cell composition and placement (described in Reimann et al. (2022)), and physiological properties, e.g., improved single cell models, multi-vesicular synaptic release, and layer-wise compensation for missing synapses (described in Isbister et al. (2023)). For this study, we used a seven column subvolume comprising 211,712 neurons in 2.4 mm³ of tissue (Figure 1A). In line with the biological variability, excitatory cells are modeled as a diverse set of morphologies (Reimann et al., 2022; Kanari et al., 2019) (Figure 1B) equipped with conductances distributed across all electrical compartments (Reva et al., 2022) (Supplementary Figure S1A). The connectivity and synaptic physiology of these cells were extensively validated (Isbister et al., 2023; Reimann et al., 2022) (Figure 1C; Supplementary Figure S1C). To deliver input with spatio-temporal precision to the circuit, we used thalamic fibers from the ventral posteriomedial nucleus of the thalamus (VPM) and the high-order posteriomedial nucleus of the thalamus (POm; Figure 1D; Meyer et al., 2010).

To simulate long-term plasticity, we integrated our recently published calcium-based plasticity model that was used to describe functional long-term potentiation and depression between pairs of pyramidal cells (PCs) (Chindemi et al., 2022). In short, the model follows the formalism of Graupner and Brunel (2012), where pre- and postsynaptic spikes led to changes in synaptic [Ca²⁺]_i (Figure 1E). Calcium entering though NMDA receptors and voltage-dependent calcium channels (VDCCs) contributes to [Ca²⁺]_i (Equation (2) in Methods). When the integrated calcium trace of a synapse crosses the threshold for depression (θ_d) or the higher one for potentiation (θ_p), synaptic efficacy (ρ) is updated (Figure 1E left; Equation (1) in Methods). Changes in ρ are then converted into changes in the utilization of synaptic efficacy (U_SE), a variable of the Tsodyks-Markram model of short-term plasticity describing the baseline probability of vesicle release (Tsodyks and Markram, 1997), and the peak AMPA receptor conductance (ĝ_{AMP A}; Equations (7) and (8) in Methods respectively). As a result of updating U_SE, short- and longterm plasticity are tightly coupled in the model (Markram and Tsodyks, 1996; Costa et al., 2015; Deperrois and Graupner, 2020). In the absence of [Ca²⁺]_i influx ρ exhibits bistable dynamics (Lisman, 1985; Graupner and Brunel, 2012) and at initialization, synapses are assumed to be at one of the two fixed points (fully depressed ones at ρ = 0 and fully potentiated ones at ρ = 1) and their assignment to these states is pathway-specific (Supplementary Figure S1C3).

We calibrated layer-wise spontaneous firing rates and evoked activity to brief VPM inputs matching in vivo data from Reyes-Puerta et al. (2015). Spontaneous activity was driven by somatic injection of a layer- and cell type-specific noisy conductance signal (Isbister et al., 2023, see Methods). Evoked activity was driven by a thalamocortical input stream already described in Ecker et al. (2024). In short, 10 VPM input patterns were repeatedly presented in random order with a 500 ms inter-stimulus interval, together with a non-specific POm input. These VPM patterns were defined with varying degrees of overlap in the sets of activated fibers (Figure 2A; see Methods). An exemplary raster plot, characterizing the evoked state of the plastic network is shown in Figure 2B.

Figure 2:

2.1 Plasticity changes promote stable network activity

In vivo, plasticity does not seem to create macroscopic network instability. Therefore, it is a critical first step to validate our model in this regard. Previous theoretical work has shown that while embedding simple STDP rules in spiking networks without fast homeostatic plasticity leads to exploding firing rates (Zenke et al., 2017a; Morrison et al., 2007), this is not necessarily the case for calcium-based rules (Higgins et al., 2014; Graupner et al., 2016; Wang and Aljadeff, 2022). To assess the stability of our network while undergoing plasticity, we simulated 10 minutes of biological time of in vivo-like, stimulus-evoked network activity and verified that excitatory firing rates remained stable throughout the entire length of the simulation (Chindemi, 2018, Figure 2C).

When comparing the amount of changes in ρ across time steps, we found that most of the plastic changes happened in the first 1-2 minutes of the simulation, after which they stabilized (Figure 2D1). While small changes were still apparent towards the end of the simulation, by visualizing individual synaptic efficacy traces we confirmed that most of them oscillated around a dynamic fixed-point and the amount of changes in the second half of the simulation were negligible (Chindemi, 2018, Supplementary Figure S5). The results were similar (Figure 2E1, black line) when we averaged ρ values within a connection, i.e., all synapses between a pair of neurons (on average 4.1 2.3 synapses per connection; Supplementary Figure S1C1). Considering changes at the connection level allowed us to compare our results against a traditional spike pair-based STDP rule (Gerstner et al., 1996; Kempter et al., 1999; Song et al., 2000, see Methods) executed on the 35 M excitatory spikes from our simulation. We observed that the STDP rule kept inducing the same magnitude of changes throughout all 10 minutes of the simulation instead of stabilizing after a transient phase (Figure 2E, pink line). To further assess whether the remaining changes in mean ρ were in the form of fluctuations around a dynamic fixed-point, we compared the changes in the second half of the simulation to a random walk with the same step size (see Methods). The magnitude of changes were well below the random walk control for our simulation, while the changes from the STDP rule were above both, but below a random walk with a step size fitted to the changes induced by the STDP (Figure 2E2, see Methods). Consequently, the lifetime of synaptic efficacy configurations in the network are longer when using a calcium-based rule than an STPD rule or a random walk. This demonstrates that not only the scale of the simulated network, its biorealistic connection probabilities and the low, in vivo-like rates contribute to the stability of changes observed in our simulation, but also the use of a calcium-based plasticity model (Chindemi et al., 2022) in line with previous modeling insights (Higgins et al., 2014).

Next, we went beyond keeping a stable network stable and investigated if an unstable network exhibiting synchronous activity (Figure 2F, left) could be brought back to a stable, asynchronous regime by the calcium-based plasticity rule. The synchronized activity state was the result of simulating a higher [Ca²⁺]_o and thus higher release probability (U_SE, see Methods). Surprisingly, the firing rate flipped from synchronous to asynchronous (Figure 2F, right) and it did so within a minute of biological time (Supplementary Figure S6B). This was likely driven by a sudden decrease in the mean ρ values over synapses in the network, which decrease could not be observed for the STDP rule (Figure 2G). As before, both rules led to smaller changes than their corresponding random walk controls (Figure 2H). Taken together, our results indicate that the calcium-based rule evolves an over-excitable network towards lower excitability, which in this case is a role that a homeostatic rule would play (Turrigiano and Nelson, 2004; Zenke et al., 2017a).

In summary, the calcium-based plasticity rule promotes stability in our network. Indeed, when applied to asynchronous activity, it generates changes in synaptic efficacy that stabilize after a transient phase without leading to exploding firing rates. Furthermore, when applied to an unstable synchronous activity regime, the mean synaptic efficacy of the whole network consistently decreases making the network less excitable and fosters stability in the macroscopic activity.

2.2 Plasticity induces sparse and specific changes driven by the stimulus and the network’s topology

We have shown that the changes in synaptic efficacy stabilize after a transient period, which raises the question: What predictions does the model make about the structure of the plastic changes in vivo? First of all, the changes were sparse, implying that the configuration of synaptic efficacy of the overall network remained largely unchanged. Indeed, we found that the number of connections changing was below 18% across the whole network and remained low when restricted to individual layers (Figure 3A1). Moreover, the propensity of changes increased as the pairwise firing rates increased (Figure 3A2), in line with previous modeling insights (Litwin-Kumar and Doiron, 2014; Graupner et al., 2016). On the level of individual synapses changes were even sparser (Figure 3B1), and in both cases depression was more common than potentiation. Layer 5 (L5) PCs contributed mostly to changes on the basal dendrites, while apical changes happened mostly on L6 PCs (Figure 3B2).

Figure 3:

So far, we characterized plastic change in terms of ρ, as this parameter is bounded to the [0, 1] interval and thus easy to interpret across pathways. Another parameter affecting the function of the synaptic connection is ĝ_{AMP A}. We found a minimal decrease (-2 pS, Figure 3C1) in its mean value due to plasticity, which was explained by frequent depression and a heavy tail of potentiation amplitudes (Figure 3C2). The distribution of ĝ_{AMP A} remained lognormal, in line with biology (Buzsáki and Mizuseki, 2014; Rößler et al., 2023, Figure 3C1, bottom). The plasticity model employed was bistable around ρ_∗ = 0.5, i.e., in the absence of further Ca²⁺. influx, values above 0.5 would converge towards 1 and values below towards 0 (see Graupner and Brunel, 2012 and Methods). The fraction of changes in ρ that crossed this threshold was slightly higher for depression (26%) than for potentiation (23%) (compare Figure 3D1 to B1). Furthermore, the vast majority of these crossings occurred early, with 60% of potentiation and 40% of depression happening within the first minute (Figure 3D2).

Additionally, we hypothesized that both the amount, i.e., the level of sparsity, and the structure of plastic changes are non-random and instead shaped by the stimuli and the underlying network’s topology. We confirmed our hypothesis first for the amount of plastic changes. We ran control simulations where we delivered random Poisson spikes on the same VPM fibers with the same mean rate, but without the spatio-temporal structure of the stimuli. We found that this reduced the number of connections undergoing plasticity by 25%, demonstrating the importance of stimulus structure over simple firing of pre- and postsynaptic neurons (Supplementary Figure S7). Additionally, without synaptic transmission, plasticity was reduced even further: In simulations of stimulus streams where intrinsic connectivity between the simulated neurons was cut we found occasional changes in ρ, but the number of changing synapses was an order of magnitude below baseline conditions (Supplementary Figure S7). Thus, our plasticity model is not strictly Hebbian, since the effect of postsynaptic firing alone could change synaptic efficacy, although presynaptic release was required for most of the observed changes (Graupner and Brunel, 2012; Graupner et al., 2016). When the external stimuli were also left out, negligible plastic changes occurred (Supplementary Figure S7).

We then confirmed that the network’s topology shaped plastic changes as well. In particular, we found layer-to-layer pathway specificity in the amount of plastic changes observed. That is, the amount of changes between layers differed from those expected in a random control with the same pre- and postsynaptic populations (Figure 4A2 vs. A3). We further quantified the difference in the subnetwork of changing connections from the random control by counting directed simplices (Fig 4B). These are motifs that have previously shown to be linked to network function (Reimann et al., 2017) as well as quantify the complexity of the network’s topology (Kahle, 2009; Bobrowski and Kahle, 2018). Succinctly, a directed simplex of dimension k, is a motif on k + 1 neurons with all-to-all feed-forward connectivity (Figure 4A1 right, see Methods). We found strong overexpression of simplices in the subgraph of changing connections compared to the control (Figure 4B). Furthermore, the maximal simplex dimension found in the subgraph was two higher than expected by chance. This result shows that the connections changing are not determined only by their pre- and postsynaptic populations but also depend on their embedding in the whole network. We will explore this further in Section 2.5.

Figure 4:

Finally, we showed that the input stimuli shaped plastic changes by running 2-minute-long plastic simulations in which we only presented a single pattern (several times, with the same 500 ms inter-stimulus interval as before) and compared the changes in mean ρ matrices. By using an input - output distance correlation analysis, we found that the distances between these mean ρ matrices were highly correlated with the input distances of the associated VPM fiber locations (r = 0.716, p < 0.0001, Figure 4C2). Furthermore, this required using Euclidean distance between the ρ matrices on the output side, implying that both the mean ρ values themselves and their position in the network are relevant. When either Hamming distance (taking only the identity of changing connection into account) or Earth mover’s distance (taking only the distribution of mean ρ values into account) were used instead, we found weaker and non-significant correlations respectively (Figure 4C3 and C4 respectively).

In summary, we observed that ~ 7% of synapses undergo long-term plasticity under realistic in vivo-like conditions in 10 minutes of biological time, and most of these synapses are on L5 PC’s basal dendrites. Moreover, the changes are not random, but depend not only on the firing rates but also on the recurrent connectivity and input stimuli. Potentiation dominated in amplitude, while depression counteracted it in frequency, keeping the network stable without needing to model homeostatic plasticity (Zenke et al., 2017a; Turrigiano and Nelson, 2004).

2.3 More frequent plastic changes within and across cell assemblies

We have shown that plastic changes are sparse and highly specific. We therefore tried to understand the rules determining which synapses changed. From the parametrization of the plasticity model, we learned that presynaptic spikes contribute orders of magnitude higher [Ca²⁺]_i than postsynaptic ones if the NMDA receptors are fully unblocked (Supplementary Figure S2B). Thus, in order to effectively depolarize the dendrites and unblock NMDA receptors, spikes at low, in vivo-like rates must be synchronized in time. Therefore, we specifically hypothesized that plasticity of connections may be structured by the membership of participating neurons in Hebbian cell assemblies (Hebb, 1949; Harris, 2005). In our previous analysis (Ecker et al., 2024), and in line with experimental results (Harris, 2005; Song et al., 2005; Perin et al., 2011), we observed that the number of afferents from an assembly, is a great predictor of a neuron’s membership in an assembly. Strong positive interactions were also found across assemblies, but only when the direction of innervation reflected the temporal order of assembly activation. These results, combined with the biophysics of the plasticity model, suggest that connections within an assembly and the ones between temporarily ordered assemblies, are expected to undergo plastic changes with a higher probability.

To test this hypothesis, we detected cell assemblies, based on their co-firing function, from the in silico spiking activity of the 10-minute-long plastic simulation using methods established by experimentalists (Carrillo-Reid et al., 2015; Herzog et al., 2021, see Methods). In short, spikes were binned and bins with significantly high firing rates were hierarchically clustered (Figure 5A). These clusters correspond to the functional assemblies of neurons (Figure 5C), with a neuron being considered a member if its spiking activity correlates with the activity of an assembly significantly stronger than chance level. Since time bins and not neurons, were clustered in the first place, this method yields one assembly per time bin and single neurons can be part of several assemblies (Figure 5B, D). Note, that in contrast to modeling studies, where assemblies are usually defined based on their strong internal connectivity (i.e., their structure Litwin-Kumar and Doiron, 2014; Zenke et al., 2015; Fauth and Van Rossum, 2019; Kossio et al., 2021) our functional assemblies are detected from the activity of the network only, eliminating a potential bias arising from links between structural properties of changing synapses and the assembly detection method.

Figure 5:

Assemblies were activated in all stimulus repetitions and a series of three to four assemblies remained active for 190 30 ms (Figure 5B). Pattern C elicited the strongest response, while pattern H the weakest, and the responses of patterns H and I were the most similar to each other, as expected, since they share 66% of the VPM fibers (Figure 2A). Assembly activations had a well-preserved temporal order - with some of them always appearing early during a stimulus, while others later - and from now on we will refer to them as early, middle, and late assemblies, and will order them in the figures accordingly.

When checking the propensity of changes within and across assemblies, we indeed found more synapses undergoing long-term plasticity (Figure 5F). While only 5% of synapses depressed in the whole dataset, we found up to 13.8% when restricting the analysis to assemblies. Similarly, compared to 2% of all synapses potentiating, we observed up to 5.2% when restricting to assemblies. Interestingly, large values were found in the off-diagonal entries (Figure 5F), i.e., synapses across assemblies underwent more plastic changes than the synapses within these assemblies. In our model, the initial ρ values are pathway-specific and highest in L4 pathways (Brémaud et al., 2007, Supplementary Figure S1C3). Therefore, early assemblies, with large numbers of L4 cells have a higher than average initial ρ (Figure 5C and E respectively), thus their synapses are more likely to depress (Figure 5F left). As early assemblies are stimulus specific, and thus not part of the same assembly sequences, synaptic depression between these cells can be seen as some kind of stimulus separation. On the other hand, late assemblies, that are predominantly composed of cells from the deep layers, have a low initial ρ (Figure 5E; Supplementary Figure S1C3) and synapses towards them are more likely to potentiate (Figure 5F right). These assemblies are mostly non-specific and participate in all assembly-sequences, thus the potentiation of their efferents means a general strengthening of the stimulus response as a whole.

Together these results indicate that, in line with 70-years-old predictions, cells that fire together wire together (Hebb, 1949). Our contribution lies in making the qualitative statement above into a quantitative one: Under in vivo-like conditions, cells that fire together more than expected have three times higher chances of changing the efficacy of their connections.

2.4 Synapse clustering contributes to the emergence of cell assemblies, and facilitates plasticity across them

In addition to co-firing, a group of presynaptic neurons is more effective in depolarizing a given dendritic branch if they all innervate the same branch, i.e., they form a spatial synapse cluster (Kastellakis and Poirazi, 2019; Farinella et al., 2014; Iacaruso et al., 2017; Tazerart et al., 2020). To quantify this, we selected the 10 most innervated L5 TTPCs (thick-tufted pyramidal cells) within a cell assembly and then detected spatial clusters of synapses, defined as at least 10 synapses within a 20 µm stretch of a single dendritic branch (Methods). Next, we grouped all synapses on these 10 selected neurons per assembly into four categories based on assembly membership of the presynaptic neuron and whether the synapse was part of a cluster or not (see exemplary clustered assembly synapses in Figure 6A). While some of the ρ values of clustered synapses within an assembly, underwent small constant changes, many of them changed at the same time (vertical stripes on Figure 6B); indicating the importance of spatial co-localization for plasticity.

Figure 6:

To systematically quantify this effect, we calculated the likelihood of plastic change in each category: (not) in an assembly and (non-)clustered, by contrasting the conditional probability of observing it in a given category with the probability of observing any change irrespective of the category (see Equation (15) in Methods; Figure 6C). Surprisingly, clustered synapses within an assembly were only likely to undergo plastic changes in late assemblies, but not in early and middle assemblies (Figure 6C black arrow). Furthermore, when comparing the amplitude of changes across conditions with a 2-way ANOVA, we found that clustered within-assembly synapses depress to a smaller degree than the other ones (Supplementary Figure S8A). To explain these trends, we repeated the same analysis on the probability of the initial value of ρ to be either 0 or 1 in each category. Initial ρ values in early and middle assembly synapses, especially the clustered ones, are very likely to be initialized as fully potentiated while synapses within the late assemblies were likely to be initialized in the fully depressed state (Figure 6D black arrow). Thus the picture emerging is as follows: early and middle assemblies are partially defined by clustered (both spatial and functional) synapses that are initialized as fully potentiated. These synapses are unlikely to change, but when they do, they depress less than the others, and would converge back to ρ = 1 in absence of activity, as they do not cross the ρ_∗ = 0.5 unstable fixed-point. These stable early assemblies can therefore function as a stable backbone amid ongoing synaptic plasticity.

In our previous investigation, we found that most changes happened across assemblies (Figure 5F), so we extended the analysis described above to cross-assembly synapses (Supplementary Figure S8B). Here, the picture was reversed: cross-assembly synapses that were part of a spatial cluster were likely to be initiated as fully depressed and then had a high chance of undergoing potentiation (Supplementary Figure S8C black arrow). Together with the previous results, this suggests that synapses between assemblies are more likely to change, which is even more pronounced if these synapses form a cluster on the postsynaptic dendrite.

2.5 Network-level metrics predict plasticity and separate potentiation from depression

So far, we have found co-firing and synapse clustering to be good predictors of plastic changes in line with previous work (Zenke et al., 2015; Kastellakis and Poirazi, 2019; Harris, 2005). Additionally, we have shown that changes are also driven by the underlying network’s connectivity (Figure 4A and B). We now elaborate on this effect using edge participation, a novel network metric of edge centrality that considers how the connections are embedded in the entire network. More precisely, the k-edge participation of a connection is the number of k-simplices that the connection is part of (Figure 7A1-Full, see Methods). Across dimensions, as edge participation increased so did the probability of undergoing plastic change (Figure 7A2). This effect was stronger for potentiation than depression (Figure 7A3). Thus, we hypothesized edge participation to be a good predictor of synaptic strength. We verified this both in our model as well as in the MICrONS (2021) mm³ dataset, an electron microscopic reconstruction of cortical tissue that combines synaptic resolution with the scale needed to calculate meaningful edge participation. Both networks had similar maximal dimension of simplices as well as ranges of edge participation values (Supplementary Figure S9A), making the results of this analysis comparable. For our model we used the sum of ĝ_{AMP A} values as a proxy for connection strength. For MICrONS, we used total synapse volume of a connection (see Methods), as a measure of connection strength (Harris and Stevens, 1989). Indeed, high edge participation led to higher connection strength in MICrONS (2021) as predicted by our modeling insights; moreover, the curves across dimensions are qualitatively the same for both datasets (Figure 7B).

Figure 7:

To gain insights into pattern-specificity, we analyzed the 2-minutes-long single pattern simulations and considered instead edge participation within the assembly-specific subgraphs (Figure 7A1-Sub). When doing so, we found an even higher predictive power, which was additionally pattern specific (Figure 7C). More precisely, we found that the probability of changes only increased with edge participation in the assembly sub-graph associated with the presented pattern or to a late assembly and that it was the strongest for the former (see Figure 7C left for pattern A for an example and right for a summary of the maxima attained for all patterns).

Finally, pattern-indegree, i.e., the number of VPM fibers belonging to a pattern that innervate a neuron is also a good predictor of the probability of a connection to change i.e., an increase of either the pre- or postsynaptic pattern-indegree of a connection leads to more frequent plastic changes (Supplementary Figure S9B1). However, given the highly non-random structure of the full network, the values of pre- and postsynaptic pattern-indegree are not independent and thus their effect in driving plastic changes is hard to disentangle (Supplementary Figure S9B2, C).

In summary, edge participation in the full network allowed us to delineate potentiation from depression given change, as well as predict connection strength, which we confirmed in the MICrONS (2021) dataset. Moreover, edge-participation within assembly subgraphs is a pattern specific predictor of plastic changes.

2.6 Increased stimulus-specificity and more reliable assembly sequences characterize the network after plasticity

As the changes are sparse, one might wonder if they cause any functional change. To study this, we compared the network’s activity before and after plasticity (see Methods). Firing rates decreased slightly as a result of plasticity (Figure 8A left), while spike correlations remained stable, in line with recent findings (Oby et al., 2019; Feulner et al., 2022, Figure 8A right). Both of these effects were more pronounced in the deeper layers (Supplementary Figure S10A). Correlations of connected pairs of neurons decreased slightly after plasticity following the decrease of firing rate (Supplementary Figure S10A). However, the highest correlations both before and after plasticity were found in the most central edges in the network i.e, those with highest edge participation in high dimensional simplices (Figure 8B1). We further verified this prediction in the MICrONS dataset (Figure 8B2, see Methods). It has been previously shown before that assemblies have many edges that are central in the network, therefore, we hypothesized that they remain functional and reliably active (Fauth and Van Rossum, 2019; Kossio et al., 2021; Pérez-Ortega et al., 2021). We showed this on the level of assembly sequences: plastic changes promote stimulus specific and reliable responses. While we found considerable overlap between assemblies detected before and after plasticity, an early assembly responding to patterns B and E before plasticity, split into two stimulus-specific assemblies (Figure 8C). Note that this is not an artifact of the thresholding applied during assembly detection as forcing the pipeline to yield the same number of assemblies instead splits a middle assembly first (Supplementary Figure S10E). Moreover, the assembly sequences responding to repetitions of a given pattern generally became more reliable after plasticity (Supplementary Figure S10B). To quantify this, we measured the Hamming similarity of assembly sequences between repetitions and found significant increases for most patterns (see Methods, Figure 8D). Interestingly, this increase on the assembly level lied in contrast with a slight decrease of the reliability of the neuronal activity on the single cell level (Figure 8E).

Figure 8:

The plastic changes we observed at the network-level are in contrast with in vitro observations at the connection level. In particular, with the redistribution of synaptic efficacy towards earlier spikes during high-frequency firing caused by changes in U_SE (Chindemi et al., 2022; Markram and Tsodyks, 1996; Costa et al., 2015; Selig et al., 1999; Sjöström et al., 2003), which under our in vivo-like, low firing rates and [Ca²⁺]_o is not in the ranges where this connection-level redistribution is relevant. Specifically, at in vitro levels of [Ca²⁺]_o the potentiated U_SE shifts the connection further into the depressing regime, causing the redistribution (Figure 8F left). In vivo the lower [Ca²⁺]_o counters this by moving the connection into the pseudo-linear regime; additionally, the low firing rates make this phenomenon less relevant in general (Figure 8F right and G respectively). Thus, while Markram and Tsodyks (1996) showed a redistribution of synaptic efficacy after plasticity at the single connection level in vitro (Figure 8F), we found a redistribution at the network level under in vivo-like conditions (Figure 8A-D, Supplementary Figure S10A).

In summary, we observed a network-level redistribution of synaptic efficacy that maintained spike correlations globally, while enhancing stimulus-specificity and the reliability of the sub-plopulations in terms of assembly activation. Furthermore, we derived and confirmed in the MICrONS (2021) dataset a prediction stating that the most central connections of the network drive the correlated activity.

3 Discussion

In this work, we aimed to understand how network structure, function, and dendritic processing interact to drive synaptic plasticity at the population level in the neocortex. Specifically, we used a plasticity model, whose parameters were determined by pairwise STDP experiments but modeled at the level of dendrites. We studied its effect at the population level under in vivo-like conditions to determine network-level rules that best predict plastic changes by simulating it in a biophysically detailed model of a cortical microcircuit that was subjected to 10 different stimuli under a variety of protocols. Our principal observations are as follows: (1) Plastic changes were sparse, affecting only 7% of the synapses. A balance between occasional, large-amplitude potentiation and more frequent depression stabilized the network without explicitly modeling homeostatic plasticity. Moreover, the changes were non-random and stimulus specific. (2) Plastic changes were largely determined by the anatomical structure of synaptic connectivity and its relation to functional units, i.e., changes were most likely between co-firing cell assemblies, at clustered synapses, between neurons highly innervated by the same input, and in central edges in the network. (3) While assemblies remained fairly stable, they become more stimulus specific and reliable which has important implications for coding and learning.

The first observation (1) is quite significant considering that we did not design the learning rule to be sparse and stable. In previous models of plastic networks, changes were not sparse and additional mechanisms were needed to keep the network’s activity stable (Litwin-Kumar and Doiron, 2014; Delattre et al., 2015; Zenke et al., 2015; Fauth and Van Rossum, 2019; Kossio et al., 2021; Zenke et al., 2017a; Turrigiano and Nelson, 2004). The machine learning community is also aware of the importance of sparse changes, as in continual learning one has to balance plasticity and stability of existing weights to avoid catastrophic forgetting (McCloskey and Cohen, 1989; Ratcliff, 1990). In recent years, they have come up with impressive techniques that mask or slow down weight updates (e.g., by a regularizer on the loss) to improve the performance of deep neural networks (Zenke et al., 2017b; Kirkpatrick et al., 2017; Mallya and Lazebnik, 2018; Frankle and Carbin, 2019), whereas in our model it emerged without any mechanisms explicitly enforcing it. Finally, demonstrating the stimulus-specificity of changes was a crucial validation that our results are not just a description of a biophysical process, i.e., plasticity, but have implications for learning in general.

For our second observation (2), we described three different kinds of plasticity rules that emerge from the underlying plasticity model based on: activity, subcellular structure, and network structure. While all of these are powerful predictors of plastic changes, none of them fully determines them. It is rather the interplay between them (and potentially additional unknown rules) that brings about non-homogeneous changes across the full microcircuit. The first two can be explained from the biophysics of the plasticity model and links our results to the classical work of Hebb (1949) as well as the recent literature on synapse clustering (Kastellakis and Poirazi, 2019; Farinella et al., 2014; Iacaruso et al., 2017; Tazerart et al., 2020). With respect to synapse clustering, we would highlight that our synapses are stochastic and the release probability between PCs is ~ 0.04 at the simulated low [Ca²⁺]_o = 1.05 mM (Markram et al., 2015; Ecker et al., 2020; Jones and Keep, 1988; Borst, 2010). Therefore, care should be taken when comparing our results with either glutamate uncaging experiments, which bypass the presynaptic machinery (Pettit et al., 1997; Losonczy and Magee, 2006), or with other modeling studies that use deterministic synapses (Farinella et al., 2014; Poirazi et al., 2003; Ujfalussy and Makara, 2020). With respect to network-based rules, previous simulation approaches have characterized the learned connectomes as having an overexpression of reciprocal connectivity and densely connected neuron motifs (Brunel, 2016; Zhang et al., 2019). We expanded on this body of previous work, by introducing cell and connection-specific metrics that directly yield the probability of observing depression or potentiation of any given connection. Most of our predictors behaved similarly with respect to depression and potentiation. Only a metric based on the embedding of individual edges in the global (or assembly-specific) network was able to separate them. This edge-based prediction was confirmed by analyzing an electron microscopic reconstruction of cortical tissue in a comparable way (MICrONS, 2021).

Our third observation (3) that assembly sequence become more reliable and stimulus specific may sound counter-intuitive considering that we observed reduced firing rates and less reliable single neurons. This shows that plasticity mechanism can contribute to the emergence of a robust population code with unreliable neurons. Crucially, this means that plasticity can make circuit responses more distinct (between stimuli by increasing specificity) but also less distinct (for a given stimulus by increasing reliability). This has important computational implications that may be explored in a future study.

To facilitate further research, we are open-sourcing our model alongside detailed instructions to launch simulations and analyze the results (see Data and code availability). Simulating the model requires a performant hardware and software infrastructure (e.g., we needed 12 M core hours to run the simulation presented in this manuscript). With respect to the second part we are continuously improving the performance and efficiency of the simulator (Kumbhar et al., 2019). The model has several assumptions and limitations, which will ideally be improved upon iteratively and in a community-driven manner. First, for our simulations we used the fitted parameters of the calcium-based model as reported in Chindemi et al. (2022). Due to various improvements of synapse parameters, the results for pairing protocols at various frequencies shift slightly while maintaining the presence of STDP (Supplementary Figure S2A). Second, experimental evidence indicates that firing bursts affects plasticity (Letzkus et al., 2006; Williams and Stuart, 1999; Larkum, 2013). However, L5 TTPC bursts were rare in our simulations, as the model is based on an early developmental stage (P14-16: juvenile rats) and firing bursts only becomes prominent as the animals mature (Zhu, 2000). Furthermore, bursts can be triggered by synaptic input onto apical dendrites (Larkum, 2013; Hay et al., 2011). Inputs to apical dendrites from high-order thalamus have been shown to gate plasticity in L23 PCs in vivo, via dis-inhibiting the distal dendrites (Gambino et al., 2014; Williams and Holtmaat, 2019). While POm input was present in our model and targeted apical dendrites, we activated them in random and non-specific ways. On the other hand, top-down inputs represent context or brain state and are thought to serve as an error or target signal for learning and thus are likely very specific (Makino, 2019). Feedback signals from high-order cortical areas also innervate apical dendrites (Harris et al., 2019), but as our model is of a single cortical area only, these inputs were absent. Missing excitatory input from other areas were compensated by somatic conductance injections (Isbister et al., 2023), which could be replaced in the future by simulating additional dendritic synapses instead. Finally, learning in the brain is orchestrated by several mechanisms in unison and some of these are beyond the scope of this paper, e.g., neuromodulation, slow homeostatic plasticity, structural plasticity, and calcium induced calcium release (Turrigiano and Nelson, 2004; Delattre et al., 2015; Magee and Grienberger, 2020).

Bottom-up modeling is often criticized for being arbitrarily complex. In our simulations, the mean value of synaptic efficacies of an over excited network reduced, comparable to Turrigiano and Nelson (2004) without introducing rate homeostasis explicitly. In the machine learning field, mechanisms are added to classical algorithms to make updates sparser, overcome catastrophic forgetting, and allow the learning of tasks sequentially (Kirkpatrick et al., 2017; Zenke et al., 2017b). The calcium-based plasticity rule employed here results in sparse updates without additional mechanisms, unlike STDP rules, which update efficacies at each spike arrival. Thus, we showed that adding biological complexity allowed us to reduce the number of mechanisms explicitly modeled. Calcium-based plasticity models existed for a long time (Shouval et al., 2002; Rubin et al., 2005), have been shown to be stable (Higgins et al., 2014; Graupner et al., 2016; Wang and Aljadeff, 2022), have been modified to mimic in vivo-like low [Ca²⁺]_o (Graupner et al., 2016), and versions more complex than ours exist (Mäki-Marttunen et al., 2020; Rodrigues et al., 2023). Our contribution lies in using it at subcellular resolution, taking dendritic processing into account (Chindemi et al., 2022) and simulating a network of this scale. Using such scales allowed us to explore not only the temporal aspects of plasticity, but also their link to spatial and network aspects, such as clusters of synapses and cell assemblies. Furthermore, the analysis methods developed on large scale functional and structural data empowered us to make and test novel predictions in the MICrONS (2021) dataset, which while pushing the boundaries of big data neuroscience, was so far only analyzed with single cells in focus instead of the network as a whole (Ding et al., 2023; Wang et al., 2023).

4 Methods

4.1 Calcium-based plasticity model

The calcium-based plasticity model is fully described in Chindemi et al. (2022), but a minimal description of it can be found below. Synaptic efficacy (ρ) is based on the Graupner and Brunel (2012) formalism, which exhibits a bistable dynamics (ρ = 0 fully depressed, ρ = 1 fully potentiated, and ρ_∗ = 0.5 unstable fixed-point) described as:

where τ = 70 s is the time constant of convergence, θ_d and θ_p are depression and potentiation thresholds, γ_d = 101.5 and γ_p = 216.2 are depression and potentiation rates and Θ is the Heaviside function. Ca^∗ is linked to the dynamics of [Ca²⁺]_i in spines (see below), which was modeled as:

where and I_{V DCC} are calcium currents through NMDA receptors and VDCCs, η = 0.04 is the fraction of unbuffered calcium, F is the Faraday constant, X is the spine volume, is the resting value of [Ca²⁺]_i, and τ_Ca = 12 ms is the time constant of free (unbuffered) calcium clearance. depends on the state of the Mg²⁺ block as follows:

where g_NMDAand E_NMDA= 3 mV are the conductance and the reversal potential of the NMDA receptor, V is the local voltage, and m describes the nonlinear voltage dependence due to the Mg²⁺ block following the Jahr and Stevens (1990) formalism:

where θ is a scaling factor of the extracellular magnesium concentration [Mg²⁺]_o, and κ is the slope of the voltage dependence. Parameters θ = 2.552 and κ = 0.072 were obtained by refitting the model to cortical recordings from Vargas-Caballero and Robinson (2003) (as opposed to the original parameters fit to hippocampal ones Jahr and Stevens, 1990).

Inspired by previous theoretical insights (Rubin et al., 2005), a leaky integrator of [Ca²⁺]_i was introduced (Ca^∗) to slow down its time course instead of modeling enzymes downstream of calcium (e.g. CamKII as others did (Mäki-Marttunen et al., 2020; Rodrigues et al., 2023)):

where τ ^∗ = 278.318 ms is the time constant of the integrator. Updates in ρ were linked to this Ca^∗ variable crossing θ_d and/or θ_p (see equation (1)). The two synapse-specific threshold were derived based on peaks in [Ca²⁺]_i caused by pre- and postsynaptic spikes, c_pre and c_post respectively. To measure these parameters for all 312,709,576 synapses, simulations of single cells were run, in which either the pre- or the postsynaptic cell was made to fire a single action potential and the local [Ca²⁺]_i was monitored in each synapse. Presynaptically evoked [Ca²⁺]_i peaks were three orders of magnitude larger, than the ones evoked by postsynaptic spikes (Supplementary Figure S2B). Postsynaptically evoked [Ca²⁺]_i peaks had a multimodal distribution in the apical dendrites (Supplementary Figure S2B right), in line with Landau et al. (2022). Since 8% of L6 PCs could not be made to fire a single action potential (only bursts), synapses on those cells (10,995,513 in total) were assumed to be non-plastic, i.e., their thresholds were set to a negative value that could not be crossed. Similarly, as the plasticity of connections between L4 spiny stellate cells was shown to be non-NMDA dependent (Chindemi et al., 2022; Egger et al., 1999) those connections were made non-plastic. For the remainder of cells, θ_d and θ_p were derived as follows:

where a_i,j and b_i,j are constants optimized during model fitting for apical and basal dendrites respectively. Changes in ρ were then converted by low-pass filtering into changes U_SE and ĝAMPA as follows:

where , and are the fully depressed (d) and fully potentiated (p) values of the given variables in-between which they evolve and τ_change= 100 s is the time constant of the conversion of changes in ρ into changes in U_SE and ĝAMPA.

In the nbS1 model U_SE is also modulated by [Ca²⁺]_o, where a reduction in [Ca²⁺]_o leads to a pathway-specific, non-linear reduction in U_SE (Markram et al., 2015; Ecker et al., 2020).

When extracting the network’s state after plasticity, not only the U_SE and ĝAMPA values, but also the peak NMDA conductances (ĝNMDA) were updated according the ρ values in the last time step of the simulation.

Model parameter optimization (Chindemi et al., 2022) and the derivation of thresholds from c_pre and c_post measurements was done with plastyfire.

4.2 In vivo-like spontaneous and evoked activity

The calibration process that leads to the in vivo-like spontaneous activity is fully described in Isbister et al. (2023), but a minimal description and a list of the parameters used in this article can be found below. As extracellular recordings are known to overestimate firing rates (Wohrer et al., 2013), a spectrum of spontaneous states at fixed percentage of the rates reported in Reyes-Puerta et al. (2015) were calibrated (Isbister et al., 2023). Matching specific firing rates in silico was achieved by iterative adjustments of layer and cell-type (excitatory/inhibitory) specific somatic conductance injection (following an Ornstein-Uhlenbeck process (Destexhe et al., 2001)). By introducing plasticity at all E to E synapses, an additional depolarizing current from VDCCs was added to the model, which made the network more active than its non-plastic counterpart (Supplementary Figure S3A). This required an algorithmic lowering of the amplitude of injected conductances from Isbister et al. (2023) to achieve the same in vivo-like layer-wise spontaneous firing rates (Supplementary Figure S3B). The spontaneous state used in the article is characterized by the parameters: [Ca²⁺]_o = 1.05 mM (Jones and Keep, 1988), percentage of reported firing rates = 40%, the coefficient of variation (CV; std/mean) of the noise process = 0.4.

The thalamic input patterns, and the spike trains delivered on them are fully described in Ecker et al. (2024), but a minimal description, highlighting the changes applied in this study, can be found below. First, the flatmap location (Bolaños-Puchet et al., 2024) of VPM fibers avoiding the boundaries of the network were clustered with k-means to form 100 bundles of fibers. Second, the four base patterns (A, B, C, and D) were formed by randomly selecting four non-overlapping groups of bundles, each containing 12% of them. The remaining six patterns were derived from these base patterns with various degrees of overlap: three patterns as combinations of two of the base ones (E, F, G), two patterns as combinations of three of the base ones (H, I), and one pattern as a combination of all four base ones (J). Third, the input stream was defined as a random presentation of these 10 patterns, in a balanced way. Last, for each pattern presentation, unique spike times were generated for its corresponding fibers following a 100 ms-long inhomogeneous adapting Markov process (Muller et al., 2007). The maximal rate of the VPM fibers was set to 17.5 Hz (compared to 30 Hz for the non-plastic circuit in Ecker et al., 2024) and half of that for POm. The overlap of the patterns is clearly visible in the firing pattern of each group of fibers corresponding to them (Supplementary Figure S4).

4.3 Network simulations

Simulations were run using the NEURON simulator as a core engine with the Blue Brain Project’s collection of hoc and NMODL templates for parallel execution on supercomputers (Kumbhar et al., 2019; Hines and Carnevale, 1997; Awile et al., 2022) (see Data and code availability). Simulating 10 minutes of biological time with reporting the state of all synapses (in every second) took 2,350,000 core hours (∼ 4x more than the corresponding non-plastic circuit without reporting), on our HPE based supercomputer, installed at CSCS, Lugano. Simulations were always repeated at least three times to assess the consistency of the results.

4.4 Control STDP rule

To compare the amount of changes induced by the calcium-based model of Chindemi et al. (2022) with classical plasticity rules, the 35,264,818 excitatory spikes from the 10-minute-long simulation were evaluated with a pair-based STDP rule (Gerstner et al., 1996; Kempter et al., 1999; Song et al., 2000). Synaptic weights evolved as follows under the STDP rule:

where t_pre and t_post are the times of pre- and postsynaptic spikes, Δt = t_post — t_pre is the difference between them; A_± = 0.05 describe the weight update, which decayed exponentially with time constants τ_± = 20 ms. The same parameters were used for the comparison with the 10-minute-long synchronous (unstable) simulation (76,025,009 excitatory spikes). The STDP rule was implemented in Brian2 (Stimberg et al., 2019).

4.5 Random walk control

Changes in ρ induced by the calcium-based model of Chindemi et al. (2022) and the STDP rule described above, were compared to a random control that takes into account global and noise-induced random changes in the network. To that end, a control model was used, where k connections of the network either increase or decrease by a noise value ɛ with equal probability in each time step. Note that, in this control, the size of the change vector at each time step is given by:

This process was modeled with a random walk in k-dimensional space of step size l, which at step N was given by:

where each is a random variable in , where k is the number of changing connections of the network and such that the step size is given by . The changes induced by the random walk at step N in the network is given by its distance from the origin at that time step, which for large enough N and K is approximately:

The step size l was determined as the average step size in the data after the 5 minute mark, when the transient activity has long passed. That is, l is given by:

4.6 Cell assembly detection

The combination of methods from Carrillo-Reid et al. (2015) and Herzog et al. (2021) yielding the assembly detection pipeline is fully described in Ecker et al. (2024), but a minimal description, highlighting the changes applied in this study, can be found below. First, spikes of excitatory cells were binned using 20 ms time bins (Harris et al., 2003). Second, time bins with significantly high firing rates were determined as crossing a threshold defined as the mean activity level plus the 95^th percentile of the standard deviation of 100 shuffled controls. These shuffled controls were less strict than in Ecker et al. (2024). Unlike in the original study, where spikes were only shifted by one time bin forward or backward (Carrillo-Reid et al., 2015), spikes were shifted by any amount. This change was introduced because the network’s response to the same patterns was more variable in the plastic simulations, and to not miss any of them, a lower threshold was more fitting. Third, based on the cosine similarity of activation vectors, i.e., vectors of spike counts of all neurons in the given significant time bins, a similarity matrix was built (Carrillo-Reid et al., 2015). Fourth, this similarity matrix was hierarchically clustered using Ward’s linkage (Pérez-Ortega et al., 2021; Montijn et al., 2016). Like for any other unsupervised clustering method, the number of optimal clusters cannot be known beforehand, thus potential number of clusters were scanned between five and twenty. The number of clusters with the lowest Davis-Bouldin index was chosen, which maximizes the similarity within elements of the cluster while minimizing the between cluster similarity (Davies and Bouldin, 1979). Fifth, neurons were associated to these clusters based on their spiking activity, and it was determined whether they formed a cell assembly or not. The correlations between the spike trains of all neurons and the activation sequences of all clusters were computed and the ones with significant correlation selected to be part of the potential assemblies. Significance was determined based on exceeding the 95^th percentile of correlations of shuffled controls (1000 controls with spikes of individual cells shifted by any amount as above; Herzog et al., 2021; Montijn et al., 2016). Finally, it was required that the mean pairwise correlation of the spikes of the neurons with significant correlations was higher than the mean pairwise correlation of neurons in the whole dataset (Herzog et al., 2021). Clusters passing this last criterion were considered to be functional assemblies and the neurons with significant correlations their constituent cells.

The reliability of assembly sequences was defined as the Hamming similarity over the repetitions of a single pattern, which attains a value of 1 if two assembly sequences are identical and 0 if they are completely different. Assemblies of neurons were compared using their Jaccard similarity.

Assemblies were detected using the assemblyfire package.

4.7 Synapse clusters and likelihood of plastic changes within them

To quantify the importance of co-localization of synapses for plastic changes, clusters of synapses were detected. To be part of a synapse cluster, a synapse was required to have at least nine other synapses on the same dendritic branch, i.e., between two branching points of the dendrite, with ≤ 10 µm (Euclidean) distance. Significance of spatial clustering was determined similarly to Druckmann et al. (2014). To that end, the distribution of synapse neighbor distances of the 10 selected synapses were compared with a Poisson model, assuming exponentially distributed inter-synapse distances, based on all (same branch) synapse neighbor distances on the given neuron. Clusters were merged in a subsequent step, thus synapse clusters with more than 10 synapses, spanning more than 20 µms were also feasible.

Plastic changes in synapse clusters were only analyzed for a small subpopulation of neurons in an assembly (10 L5 TTPCs per assembly), which were selected based on maximizing two connectivity features: assembly-indegree and synaptic clustering coefficient (SCC). The first one is the number of connections from a presynaptic assembly, while the second one quantifies the colocalization of synapses on the dendrites of a neuron from its presynaptic assembly. Maximizing these two metrics, which were introduced (and are described in detail) in Ecker et al. (2024), allows one to select subpopulations with high probability of finding synapse clusters.

Control synapse clusters, originating from non-assembly neurons were also detected on the same postsynaptic neurons.

The normalized likelihood of changes, conditioned on the four categories a synapse could fall into (assembly clustered, assembly non-clustered, non-assembly cluster, non-assembly non-clustered) were quantified using the Michaelson contrast, defined as:

where changed was split to be either potentiated (Δρ > 0) or depressed (Δρ < 0). Note that a nonzero value for one category always has to be compensated by a nonzero value with the opposite sign in another. For the normalized likelihood of initial ρ values, the same equation was used but Pr(Δρ > 0| category) was replaced with Pr(ρ₀ = 1| category) for potentiated and Pr(Δρ < 0| category) with Pr(ρ₀ = 0| category) for depressed synapses respectively.

4.8 Topological metrics

For quantifying the non-random nature of changing connections, directed simplex counts were used. A k-simplex in a network G is a set of on k + 1 nodes (neurons) of G that are all-to-all connected in a feedforward fashion. That is, there is a numbering of the nodes v₀, v₁, . . . v_k, such that for every i < j there is an edge from v_i to v_j, and k is called the dimension of the simplex. In particular, 0-simplices are the nodes of the network, 1-simplices directed edges (connections) and 2-simplices are also known as transitive triads. Random controls were defined as the same number of edges between the same nodes, resulting in the same 0- and 1-simplex counts. Given an edge in G, a notion of its edge-centrality in the network is its k-edge participation, which is the number of k-simplices the edge belongs to. This extends the classic notion of node participation from nodes to edges. This value can be computed either for the simplices in the entire network or in an specificied subnetwork e.g., the subnetwork on neurons belonging to an assembly.

Simplex counts and edge participation values were computed using the connalysis package (Egas Santander et al., 2024), based on a fast C++ implementation from flagsercount (Lütgehetmann et al., 2020).

4.9 MICrONS dataset

For the comparison of two of our findings with a rodent electron microscopy dataset with coregistered calcium-imaging traces, the minnie65_public release of the MICrONS (2021) dataset was used. Synapses from sources other than one of the 60,048 classified neurons inside the bigger (minnie65) volume were discarded. The 117 version of the structural dataset was used for these analysis and made freely available in SONATA format (Dai et al., 2020) at https://doi.org/10.5281/zenodo.8364070. For comparable analysis we restricted the volume to its central part (650, 000 x 950, 000 and 700, 000 z 1, 000, 000) and considered only E to E connections. Synapse volume was defined as the number of voxels painted by the automatic cleft segmentation (from the correspondings synapses_pni_2 table) times the 4x4x40 nm voxel size. Total synapse volume was defined as the sum across all synapses mediating a connection.

For the functional dataset, eight calcium-imaging sessions from the 661 version were used. These sessions were selected by the same criteria as in Reimann et al. (2023); Egas Santander et al. (2024), namely at least 1000 neurons were scanned in each session and at least 85% of them were co-registered in the structural data. Neurons with non-unique identifiers were dropped, and the remaining ones restricted to the central part of the structural volume described above. Deconvolved spike trains of neurons (243 neurons/session on average, 1817 unique neurons in total, 2% and 15% of all cells respectively) were used to calculate Pearson correlations. To average spike correlations from different session without any bias, they were first z-scored within each session.

4.10 Spike time reliability

Spike time reliability was defined as the mean of the cosine similarities of a given neuron’s mean centered, smoothed spike times across all pairs of repetitions of a single pattern (Schreiber et al., 2003; Cutts and Eglen, 2014). To smooth the spike times, they were first binned to 1 ms time bins, and then convolved with a Gaussian kernel with a standard deviation of 10 ms as in Ecker et al. (2024); Egas Santander et al. (2024).

4.11 Visualization

Rendering of cells from the cortical network on Figure 1A was done with Brayns, while the rendering of selected cells with their thalamo-cortical synapses on D with BioExplorer. Morphologies on Figures 1 and 6 were rendered with NeuroMorphoVis (Abdellah et al., 2018).

Acknowledgements

The authors thank Nicolas Ninin for his involvement in the early stage of this project, Elvis Boci and Cyrille Favreau for their help with visualizations, Michael Gevaert, Joni Herttuainen and Thomas Delemontex for their assistance with software engineering, and Alberto Antonietti, Christoph Pokorny, Kathryn B. Hess, Ran Levi, Wulfram Gerstner, Roberto Araya and Henry Markram for discussions.

Funding

This study was supported by funding to the Blue Brain Project, a research center of the École polytechnique fédérale de Lausanne (EPFL), from the Swiss government’s ETH Board of the Swiss Federal Institutes of Technology.

Author contributions

• Conceptualization: A.E., D.E.S., G.C., E.B.M., M.W.R.

• Formal analysis: A.E., D.E.S., M.W.R.

• Investigation: A.E., S.B.-P., G.C.

• Methodology: A.E., D.E.S., G.C., J.B.I., E.B.M., M.W.R.

• Project administration: M.W.R.

• Software: A.E., D.E.S, S.B.-P., G.C., J.K., P.K., M.W.R.

• Resources: J.B.A, I.M.

• Supervision: M.W.R.

• Validation: A.E., D.E.S., D.P.G.M, S.B.-P., J.B.I., M.W.R.

• Visualization: A.E., M.A.

• Writing - original draft: A.E., D.E.S., M.W.R.

• Writing - review & editing: A.E., D.E.S., S.B.-P., G.C., J.B.I., E.B.M, M.W.R.

Declaration of interest

The authors declare no competing interests.

Data and code availability

The 2.4 mm³ subvolume of the juvenile rat somatosensory cortex, containing 211,712 neurons and 312,709,576 plastic synapses in SONATA format (Dai et al., 2020) is freely available at: https://doi.org/10.5281/zenodo.8158471. It can be loaded and instantiated in CoreNEURON (Kumbhar et al., 2019) with neurodamus. The circuit and the simulations can be analyzed using Blue Brain SNAP and ConnectomeUtilities built on top of it. Cell assemblies were detected and can be analyzed with assemblyfire. Exemplary Jupyter notebooks using the packages above were deposited in the same repository on Zenodo.