Excitable neuronal assemblies with adaptation as a building block of brain circuits for velocity-controlled signal propagation

The time scale of neuronal network dynamics is determined by synaptic interactions and neuronal signal integration, both of which occur on the time scale of milliseconds. Yet many behaviors like the generation of movements or vocalizations of sounds occur on the much slower time scale of seconds. Here we ask the question of how neuronal networks of the brain can support reliable behavior on this time scale. We argue that excitable neuronal assemblies with spike-frequency adaptation may serve as building blocks that can flexibly adjust the speed of execution of neural circuit function. We show in simulations that a chain of neuronal assemblies can propagate signals reliably, similar to the well-known synfire chain, but with the crucial difference that the propagation speed is slower and tunable to the behaviorally relevant range. Moreover we study a grid of excitable neuronal assemblies as a simplified model of the somatosensory barrel cortex of the mouse and demonstrate that various patterns of experimentally observed spatial activity propagation can be explained.


Author summary
Models of activity propagation in cortical networks have often been based on feedforward structures. Here we propose a model of activity propagation, called excitation chain, which does not need such a feedforward structure. The model is composed of excitable neural assemblies with spike-frequency adaptation, connected bidirectionally in a row or a grid. This prototypical neural circuit can propagate activity forwards, backwards or in both directions. Furthermore, the propagation speed is slow enough to trigger the generation of behaviors on the time scale of hundreds of milliseconds. A two-dimensional variant of the model is able to generate different activity propagation patterns, similar to spontaneous activity and stimulus-evoked responses in anesthetized mouse barrel cortex. We propose the excitation chain model as a basic component that can be employed in various ways to create spiking neural circuit models that generate signals on behavioral time scales. In contrast to abstract models of excitable media, our model makes an explicit link to the time scale of neuronal spikes.

Introduction
Reliable propagation of activity is necessary for processing and transmitting sensory signals in the brain. During the last two decades, two prominent types of computational models have been studied to address this issue. First, the synfire chain consists of groups of spiking neurons connected in a feedforward architecture [1][2][3][4] potentially embedded in recurrent networks [5,6]. Second, rate propagation models [6][7][8][9] use a similar feedforward architecture, but instead of spikes they propagate fluctuations of the firing rate. In a synfire chain the refractory behavior of neurons after firing a spike is the crucial element of stable activity propagation [10]. The derivative of the membrane potential shapes spike density and sharpens the activity pulse [10,11]. Statistical methods [12,13] have been proposed for detecting the firing pattern of synfire chains in spike train recordings. Reliable neuronal firing patterns compatible with synfire chains have been observed in area HVC of the song-bird [14], while the statistical significance of synfire chains in cortical neurons is questionable [15]. Synfire chain models have been used for reproducing behavioral functions such as bird song generation [16] and monkey scribbling [17]. Systems of interacting synfire chains were also used for building a large-scale model of cortex [18,19]. While both synfire chain and rate propagation succeeded to model fast behaviors (behaviors on the order of milliseconds), they are expensive in terms of neuron numbers and therefore not suitable for reproducing behavioral phenomena that need a slower, and sometimes tunable, speed of activity propagation. Essentially synfire chains implement a clock, set by the delay of spike propagation and the rise time of excitatory postsynaptic potentials on the millisecond time scale [20,21]. In order to address this issue, we propose an excitation chain model, for activity propagation in a bidirectional chain of neuronal assemblies. Our model can be considered as a spiking version of excitable media [22][23][24][25][26], with an explicit link to the neuronal time scale of spikes.
In contrast to previous synfire chain and rate propagation models, our model does not require an explicit feedforward architecture. Specific feedforward structures have not been observed in experiments in the neocortex so far. Here we propose an excitation chain model that is consistent with the following experimental connectivity data: first, inter-and intraassembly connection probability and synaptic weights with values in the experimentally observed ranges [27,28]; and second, clustered connectivity of neurons [29]. Hence, although the concept of an excitation chain as such is a rather generic and abstract model, its basic connectivity features are consistent with experimentally observed properties.
In the next section, we describe our model and its dynamics in detail. We also illustrate how activity propagation can be made faster or slower by changing synaptic weights. Then we analyze the behavior of the model and explain the role of its elements in forming the dynamics. Finally, we extend the excitation chain to a two-dimensional model, which we may call an excitation grid. We argue that this grid of assemblies can be considered as the skeleton of barrel cortex, which can generate different spatio-temporal modes of activity propagation observed experimentally in barrel cortex [30,31].

The speed of activity propagation in a chain of excitable bistable assemblies
In order to propagate a signal of activation through several excitable assemblies, we first connect several groups of neurons in a bidirectional chain (Fig 1A). Each group contains an excitatory assembly and a population of inhibitory neurons. Assemblies are defined as small populations of neurons with high connection probability. More precisely, inside each excitatory assembly synapses are strong and the connectivity is high (connection probability p = 50%) whereas each inhibitory population has smaller synaptic weights and lower connectivity. Within each group, the excitatory assembly and inhibitory population are mutually connected. Moreover, the excitatory assembly is connected to the inhibitory population and excitatory assembly of neighboring groups. In contrast, the inhibitory populations do not have inter-group connections (see Materials and methods for the details and parameters).
Basic mechanism and functionality. If we stimulate the excitatory assembly of the first group of the chain with a transient stimulus of 25ms duration (see Materials and methods), all neurons of this group fire several spikes. The activation of this first assembly is then propagated step by step through the chain of assemblies until the last group. Fig 1B shows the raster plot of all excitatory neurons as well as the population-averaged activity of excitatory and inhibitory populations of the chain. One can see that despite the reciprocal connections between excitatory assemblies, the activity is propagated in a feedforward manner. We repeated the simulation several times with different transient stimuli (S1 Fig). Whenever the transient input stimulus is able to activate the first excitatory assembly, the activity is reliably propagated through the chain to the last group. For a vast range of parameters (Fig 1C), we have observed neither a termination of the activity wave nor an instability in the propagation dynamics (such as convergence to synchronous firing of all excitatory neurons).
If we stimulate the excitatory assembly of the last group, we see that the activity propagates backwards (Fig 1.inset). The excitation wave can also spread in both direction simultaneously. Stimulating an excitatory assembly in the middle of the chain produces two traveling waves, one towards the beginning of the chain and another one towards its end (Fig 1.inset). The property of activity propagation in different directions has been observed in multi-electrode extracellular recordings of the neocortex. For example, based on the place of local application of glutamate, neural firing is initiated in a forward, backward or bidirectional manner [32].
Speed depends on synaptic weights. The speed of activity propagation in our excitation chain is much lower than in a synfire chain [3]. In synfire chains the time needed for the activation to jump to next group is on the order of the synaptic transmission delay (1-5ms), while in our model this time is roughly 13-55ms, although we have used a short synaptic transmission delay of 1ms. This slow propagation allows us to describe phenomena on the time scale of several hundreds of milliseconds or even seconds.
The speed of propagation in the excitation chain is controlled by the inter-assembly synaptic weights. Let us first define how to measure the delay (and consequently the speed) of activity propagation. For the sake of simplicity, we define the activation time of each excitatory assembly by the average time of the first spike of each neuron in the assembly. Alternatively, and without change of the measured propagation delay, we could also define the activation time of an assembly of N neurons as the time when N/2 neurons have fired a first spike after an interval of at least 100ms. Importantly, the difference of the activation time of two neighboring assemblies can be considered as the time needed for transmitting the activity signal from one group to the next group. The inverse of the activation time difference of the first and the last assembly can be used as a measure for propagation speed.
We use the symbols w exc to denote the synaptic weights between neighboring excitatory assemblies and w inh to denote the synaptic weight from excitatory assemblies to the neighboring inhibitory populations. Fig 1C, 1D and 1E show that increasing w exc increases the speed, while increasing w inh reduces it. The activity wave remains stable and propagates reliably over a broad range of parameter values.
Removing the inhibitory populations. Excitatory assemblies are the essential elements of the excitation chain, while the role of inhibitory populations in the chain consists mostly in reducing the propagation speed. Therefore, we may simplify our model by removing all inhibitory populations. A chain of excitatory assemblies only (Fig 2A) is able to propagate the activity in different directions (Fig 2B) similar to a chain containing also inhibitory neurons. The propagation speed is regulated by modifying w exc (Fig 2C). However, because of the lack of inhibition, the speed cannot be lower than a critical value below which transmission becomes unreliable. In our simulations, we found a maximum delay of 34ms instead of 56ms with inhibition.

Analysis of excitation chain dynamics
Excitation chains rely on bistable assembly dynamics. In order to understand the dynamics of the chain and identify the components that determine the propagation speed, we first focus on one assembly of excitatory neurons. The dynamics of each assembly can be described by self-consistent equations relating the firing rate of the assembly to the average synaptic input of the neurons (see Materials and methods). If the assembly has a high network feedback coefficient C fb (Eq 9), the dynamics of the system has three fixed points (Fig 3A-top): the low point which is a stable fixed point with zero firing rate, the switch point which is the unstable middle fixed point and a high-rate fixed point which is called the high point. If the assembly is driven by synaptic input greater than the switch point current (I s ), it approaches the high point and produces a high firing rate. Since the intra-assembly synaptic weights (w exc ) and connection probability (p) are high, the network feedback coefficient (C fb / pw exc ) is also high. Therefore, the dynamics of the assembly can be explained by this three-fixed-point configuration, which we call the excitable mode of the assembly. However, because of spike-frequency adaptation of our excitatory neuron model, the frequency of spike emission progressively decreases during the high-rate state. Consequently, the neurons' gain function changes gradually (Fig 3A-bottom) and the system goes to a new configuration which has only one fixed point, the low point. (Note that, for the same reason, the assembly does not fully converge to the high point as it is shown in Fig 3A-top. While the assembly is approaching the high point, changes of the gain function move the position of the fixed point.) Therefore the assembly eventually becomes quiet and stops firing. We refer to this configuration as the dormant mode of the assembly. In this mode, receiving synaptic input does not activate the assembly. It takes a while for the assembly to recover from the dormant mode and return to the three-fixed-point configuration, which is the excitable mode. Previous work [33] addressed the dynamics for a simpler adaptive integrate-and-fire neuron model with similar analytical approaches.
Synaptic weights determine the propagation speed. For the second step of the analysis, we take into account the interaction of neighboring assemblies in the chain. Consider the case of a chain of excitatory assemblies only. When assembly 1 goes to the high point and each of its neurons fire several spikes within a short time, it sends strong synaptic input to assembly 2. However, since the synaptic weights between assemblies are relatively low, early volleys of spikes of assembly 1 do not yet suffice for assembly 2 to cross the switch point. Therefore it takes some time for assembly 2 to accumulate enough input from assembly 1. This is the reason why the propagation of activity is slow in the excitation chain and the difference of activation time is much higher than the synaptic delay. If we increase the inter-assembly weight, then a smaller number of spikes of assembly 1 is needed to produce the switch current in assembly 2, which means it will be activated sooner (Fig 3C). This explains the increase of the propagation speed along the chain when we strengthen the inter-assembly synapses (Fig 2C).
Inhibitory neurons delay the activation of each group. Let us now discuss the effect of adding the inhibitory populations to the chain. Each inhibitory population sends inhibitory  (Fig 1B). The propagation speed can be tuned (C) by modifying the synaptic weights between excitatory assemblies (w exc ), albeit in a smaller range, than in  The network feedback (C fb , Eq 9) affects the dynamics of the system. The red curve is the noisy gain function of the GIF neuron model (mean spike count in a group of 50 independent neurons over 10ms, divided by 50 × 10ms, shaded area marks ±3 SEM) measured during the initial 10ms after switch-on of a synaptic current of mean hI syn i. The green lines (solid, dashed and dash-dotted) show the relation of firing rate and synaptic current caused by network feedback (see Materials and methods, Eq 8) for increasing C fb . The slope of the green lines is inversely proportional to the effective feedback coefficient C fb of the population. Intersections of the red curve with one of the green lines indicate potential stationary states (fixed points) of a network of non-adapting neurons. A-Bottom) The noisy gain function of adapting neurons during the first 10ms after stimulus onset (solid red curve) is different from that later (dashed red curves). B) Synaptic current received by the second assembly (averaged over the assembly's neurons) in two chains with different inter-assembly synaptic weights. The thin lines show the averaged synaptic current each neuron receives from the previous assembly, while the thick lines show the total synaptic current received from both the assembly itself and the previous assembly. When the thick line separates from the thin line, the assembly starts to fire spikes on its own. This is the moment when the assembly crosses the switch point and approaches the high-activity fixed point. This occurs earlier if inter-assembly synaptic weights are increased (blue). Consequently the propagation speed along the excitation chain increases. C) The total synaptic current received by the second assembly (black) is the sum of excitation (blue) from the first assembly and inhibition (red traces show the absolute value of inhibitory currents) from the second inhibitory population. When the total input (black) crosses the threshold, the assembly switches to the high-point and becomes activated. Increasing synaptic weights from the first assembly to the second inhibitory population increases the inhibition and decreases the total input received by the second assembly. Hence it delays the activation time. D) The analytical approach (green line, see Materials and methods) approximates the difference between the activation time of two consecutive assemblies with an error of less than 4ms. The black line shows the value of the difference obtained by averaging over 10 simulations (error bars indicate standard deviation).
input to the excitatory assembly in the same group and thus reduces the effect of the synaptic input provided by the neighboring excitatory assemblies. Therefore, the inhibition delays the time at which the switch current is reached. If we increase the inter-group synaptic weight onto inhibitory neurons w inh , the inhibitory neurons fire more often and produce more inhibition. This, in turn, increases the amount of synaptic input needed from the neighboring excitatory assembly and further delays the activation time. Fig 3C illustrates the effect of w inh on the synaptic inputs.
Adaptation gives rises to directed activity propagation. Suppose that assembly 1 has become active and sends synaptic current to assembly 2. When excitatory assembly 2 receives enough synaptic current to cross the switch point, it becomes active at a high firing rate. It then sends synaptic current to its neighbors, excitatory assembly 1 and 3. However, by that time, excitatory assembly 1 is about to fall into the dormant mode and cannot generate many spikes. Therefore, only excitatory assembly 3 switches to the high point. This procedure repeats until the end of the chain. Therefore activity propagates in one direction only although connections between assemblies are bidirectional. An analogous situation happens for the backwards propagation. In case of stimulating an assembly in the middle of a previously quiet chain, since both of its neighbor assemblies are in the excitable mode, they both switch to the high point. Consequently, the activity propagates in both directions from then on.
Spike-frequency adaptation is responsible for the transition to the dormant mode by progressively changing the gain function ( Fig 3A-bottom) during the active phase of an assembly, and eventually for its termination. By modification of the adaptation parameters of excitatory neurons, we are able to adjust the duration of the activate phase of each assembly [34].
The effect of other synaptic weights on the dynamics. After having analyzed the effect of inter-group synaptic weights (w exc and w inh ) on the propagation speed of the chain, we now focus on the weights of inhibitory to excitatory neurons inside the same group. Since intragroup inhibition contributes to the total inhibition of the assembly, it can have effects similar to w inh on the propagation speed. In order to avoid redundant parameter search, we kept the intra-group inhibition constant and explored w inh . Likewise, the intra-group excitatory to excitatory connections (connections inside each assembly) are also important. As we mentioned earlier, assemblies should have a high network feedback coefficient (C fb / pw). Otherwise, they would not be able to switch to the high point and produce high firing rates in case of receiving relatively low synaptic input.
Both connection probability and synaptic weight of the intra-group excitatory to excitatory connections affect the network feedback coefficient and therefore shape the core of the excitation chain. Other intragroup connections (excitatory to inhibitory and inhibitory to inhibitory connections) are less important. However too much inhibition may shut down the assembly by finishing its activation rapidly so that not enough synaptic input arrives at the next assembly. Therefore, it may lead to a loss of signal propagation.
Obtaining the propagation speed using an analytical approach. The self-consistent approach relating the firing rate to the average synaptic input which we mentioned earlier is useful for a qualitative explanation of the dynamics of the excitation chain. However, it is not suitable for a quantitative calculation of the propagation speed, because we cannot calculate the exact gain function for this neuron model (see Materials and methods) if the effects of spike-frequency adaptation become strong. Therefore, we developed another analytical approach (see Materials and methods) in order to obtain the difference between activation times of two consecutive assemblies in the chain of only excitatory assemblies (Fig 2). Fig 3D  compares the results of simulation and the analytical approach for different values of w exc . Our theory estimates the activation time difference with an error of less than 4ms.
Embedding short-term plasticity in the model. We can also add short-term facilitation and depression [35,36] (see Materials and methods) to the model. In these cases, we are able to adjust the propagation speed by manipulating the parameters of facilitation and depression, while keeping w exc and w inh fixed. In the first variation, we added facilitation to the excitatory synapses between assemblies. The case of no facilitation that we considered earlier corresponds to a choice of U = 1 for the usage parameter U of short-term plasticity. We observed that decreasing the value of U to values below one increases the difference between activation times of two consecutive assemblies ( Fig 1F). In the presence of facilitation, the amplitude of postsynaptic currents (PSC) is initially lower compared to case of not having facilitation. When a presynaptic neuron fires several spikes within a short interval, the amplitude of each PSC in the post-synaptic neuron increases. Only after a suitable number of spikes, the PSC amplitude reaches a stationary value. Therefore, it takes time for an assembly to provide sufficient amount of synaptic input to activate the next assembly. Modification of the recovery time constant τ rec , however, does not affect the propagation speed (S2 Fig).
In the second variation, we neglected facilitation and added depression in inter-group excitatory to inhibitory connections. We also increased the intragroup inhibitory to excitatory synaptic weight by * 7 times (1.07mV instead of 0.16mV). The propagation delay decreased as we increased the value of U ( Fig 1F). The reason is that a large amount of inhibition from the inhibitory subpopulation does not allow the excitatory assembly to become active. Depressing synapses from one assembly to the neighboring inhibitory subpopulation reduce the amount of PSC onto inhibitory subpopulation over time so that the activity of an inhibitory subpopulation and its projecting inhibition drop off. Similar to the previous case, the time constant (τ facil ) does not affect the propagation speed (S2 Fig).
Note that we still need spike-frequency adaptation in the above cases because it is the adaptation that terminates the activation of the assemblies. Without adaptation, an assembly switches to active and remains active for the rest of the simulation. It is important to notice that our model is not based on a competition between assemblies via inhibition (see Discussion) and every assembly switches to the low point on its own. If we want to remove the adaptation and preserve the functionality of the chain, we can add short-term depression to intraassembly synapses. We will come back to this point in the Discussion section.
In addition to connectivity, the duration of the active phase of assemblies plays an important role in the propagation speed. The duration should be long enough, such that each assembly can provide the synaptic input needed to activate the next assembly in the chain. Therefore, if we want to achieve a slower propagation speed, we need to increase the duration of the active phase.
The duration of the active phase can be adjusted by modification of adaptation parameters of excitatory neurons [34]. Each time a neuron fires a spike, several adaptation processes on several time scales are triggered and generate both a spike-triggered current and an increase in firing threshold [37]. To describe these adaptation effects we use two exponentially decaying kernels, η k (t) (k = 1, 2) for spike-triggered currents and γ k (t) (k = 1, 2) for dynamic threshold. In Fig 1G we decreased the amplitude of kernel γ 1 (t) (the kernel with the shorter time constant). Then we used short-term facilitation for changing the propagation speed (similar to Fig 1F). For each value of the kernel amplitude, we found the value of the usage parameter U which causes slowest propagation. Fig 1G summarizes the activation times of assemblies for different combinations of parameters. We achieved a maximum delay of 95ms between one assembly and the next using this approach. Hence, we conclude that a linear chain of 11 assemblies is sufficient to cover a typical behavioral time scale of above 1 second.
The dynamic threshold (γ(t)) has two decaying kernels, one with shorter time constant (tens of milliseconds) and the other one with longer time constant (several hundreds of milliseconds). The shorter time constant determines the duration of the active phase of assembly (see [34]). The longer one affects the time that an assembly needs to recover from the dormant mode after termination of active phase. The kernels of the spike-triggered current (η(t), see Materials and methods) have similar effects. In fact, we could remove one of these two mechanisms (dynamic threshold or spike-triggered current) and preserve the functionality of the chain. However, since the original model used both mechanisms [38,39] to fit experimental data, we preferred to keep both of them.
Reduction of firing rate by using a second type of inhibitory neurons. The firing rate of assembly neurons during the active state is high (* 200Hz). In order to regulate the firing rate, we use a second type of inhibitory neurons for each assembly (Fig 4A). These additional inhibitory populations do not receive connections from the neighboring assemblies, but form only intra-group connections. Hence, they generate inhibition only after the assemblies become active, in contrast to the inhibitory groups in Fig 1. The idea of having the second type of inhibitory neurons is consistent with the observation that cortical networks contain different types of inhibitory neurons with different electrophysiological properties, different connectivity schema and probably different roles [40]. Fig 4B shows that in the presence of second inhibitory populations peak firing rates of excitatory neurons are reduced to * 100Hz.

A grid of assemblies as a skeleton for barrel cortex
While the previous section focused on a one-dimensional structure, in this section we consider an excitation wave in a two-dimensional grid structure inspired by the layout of barrel cortex. Barrel cortex processes sensory information from the whiskers and is a part of the rodent somatosensory system. It consists of vertical modules called barrel columns, each of which relates to one principle whisker [41]. Here we make a multicolumn model of barrel cortex which contains 25 columns, organized in the shape of 5 arcs while each arc contains 5 rows. The actual mouse barrel cortex includes 33 columns with an arc of 4 rows, 4 arcs of 5 rows and 3 arcs of 3 rows [41]. For simplicity, we only consider one cortical layer of the barrel cortex. In our model, every column consists of excitatory and inhibitory neurons. Excitatory neurons are divided into two groups, a minority of assembly neurons and a majority of non-assembly neurons (see Table 2). Fig 5 shows the schematic of the model. While assembly neurons have high internal synaptic weights and connection probability, the connections between assembly and non-assembly neurons as well as connections inside non-assembly neurons are sparse and weak. Inside a column, all three groups have connections to each other, but inter-column connections are different from intra-column connections. (see Table 2). Columns are identical in terms of number of neurons, neural parameters and connections between neurons. Inhibitory neurons in our model form only intra-column connections and are not connected to the neurons of other columns, consistent with the common assumption that inhibitory neurons send short axons and contact only local targets [42]. In contrast, excitatory neurons of each column of our model connect to the excitatory neurons of the four nearest-neighbor columns. However, due to the relatively long distance between two neighboring columns, the connection probability between columns is taken as low (p = 10%), consistent with experimental observations [43].
Just as in the excitation chain of the previous section, the dynamics of each assembly can be explained by the three-fixed-point configuration. However, the subpopulations of nonassembly neurons have low network feedback coefficients (Eq 9) so that their dynamics exhibit only the low-activity fixed point. Hence the assemblies govern the dynamics of the grid model while non-assembly neurons follow the dynamics at a lower firing rate. Neglecting nonassembly and inhibitory neurons for a moment, we can consider this barrel cortex model  (Fig 1B), while the chain is still able to propagate the activity forward, backward and in both directions. merely as a grid of assembly neurons. We may think of this grid as the skeleton of our model for barrel cortex. Just as in the chain model, the synaptic weight between assemblies determines the propagation speed. Here, we adjusted the value of the inter-column synaptic weight (w exc = 0.32mV) such that the speed of activity propagation in the model is similar to that observed in the experimental data [30,31,44]. We did not need to use short-term facilitation or depression to achieve the desired speed. The remaining parameters of the model are reported in section Materials and Methods.
Qualitative comparison with experiments. The grid model is able to reproduce several aspects of the dynamics of anesthetized barrel cortex in the stimulus-evoked and spontaneous regime. In the stimulus-evoked experiments [30,[45][46][47], a sensory signal was triggered by a brief deflection of a whisker. Voltage-sensitive dye imaging showed that the neural activity started in the barrel column corresponding to the stimulated whisker and propagated to the neighboring columns. After spreading over the whole field of view of barrel cortex, the activation vanished. Fig 6A shows the simulated evolution of neuronal firing rate in the grid model after stimulating the neurons of the central column. The dynamics of the model are similar to the experimental recording. For better visibility, the voltage traces of the model were temporally filtered with a Gaussian function (σ = 30ms). We show only 64 neurons of each column: each panel in Fig 6A shows 25 squares (= columns) and each square contains 64 pixels (= neurons). These neurons are randomly selected from all neurons of the column (excitatory assembly and nonassembly neurons and inhibitory neurons). While the assembly neurons receive a high amount of synaptic input (due to their strong synaptic weights) and show a high firing rate, nonassembly and inhibitory neurons receive weaker weights and show lower firing rates (0-20Hz). Since assembly neurons form only a minority of all neurons in a column (see Materials and methods), the mean firing rate averaged across all neurons in a column is close to the firing rate of non-assembly and inhibitory neurons.
For the case of spontaneous activity, experiments showed that neural activity started on one side of the barrel cortex and only in a few columns. Then it moved from column to column in a circular fashion [31,44]. With an appropriate initial state, our model is able to reproduce the circulation of activity, on the same time scale as in experiments. If we stimulate a model column on the upper-left side, the activity starts circulating to the down-and the leftward direction (Fig 6C). The activity orbits around the central column several times before it terminates.
Petersen et al. [30] found that in some experiments, after the stimulation, the activity spread faster along the row (horizontal spread in our model) than the arc (vertical spread in our model). The authors suggested that this could be caused by axons of excitatory neurons that extended farther along a row than along an arc. For an alternative interpretation of this phenomenon we assume that the connection probability between columns connected in the row direction (p = 15%) is greater than the connection probability along the arc direction (p = 10%). Simulations of this modified model exhibit a higher propagation speed along rows comparing to arcs (Fig 6B).
Qualitative analysis of activity wave. To understand the model dynamics we focus again on a regular grid of assemblies. As an aside we note that non-assembly neurons are subsidiary and passively follow assemblies. Since they receive weak synaptic input, they fire at a lower rate than neurons inside an assembly. We added non-assembly neurons in order to have the same number of neurons as observed in experiments [27] and to lower the mean firing rate averaged across a column to realistic values; but the dynamics of the network is solely driven by the assemblies. Similar to the assemblies of the excitation chain, each grid assembly is either in the dormant or the excitable mode depending on the adaptation level of its neurons. The dormant mode corresponds to high values of neuronal adaptation variables (spike-triggered current and firing threshold). Neurons of an assembly enter the dormant mode once they have emitted a burst of spikes (or if high initial values have been assigned to these variables). After recovering from the dormant mode, neurons of the assembly have a weak spike-triggered current and relatively low firing threshold. Hence the assembly has returned to the excitable mode and can now generate a burst of spikes in case of sufficient excitation. After activation, it switches again to the dormant mode.
In our model, we manipulate the initial value of the firing threshold kernel with longer time constant (γ 2 (t)) of neurons in the assembly in order to set the initial mode of each assembly. High initial values cause the dormant mode, while low initial values correspond to the . The difference in connectivity causes different propagation speed along the row compared to the speed along the arc. Note that in (A) both connectivities were 10% and the propagation speeds were identical. C) Circulation of activity in the multi-column model. After stimulation of a corner column, the activity propagates between columns and circulates across the model for several rounds. Eventually activity vanishes. This type of dynamic has been observed in spontaneous activity in mouse barrel cortex in vivo [31]. D) Initial values of the second adaptation kernel (γ 2 ), which describes spike-triggered movement of the firing threshold (see Eq 3). For the stimulus-evoked response (shown in A) all initial values are drawn from a normal distribution with zero mean and standard deviation of 3mV. Any negative values are clamped to zero. For activity circulation (shown in C) we select the initial values to create a preferred direction of motion of the activity wave as follows. For each column, we first choose a mean value for the normal distribution, while the standard deviation of all distributions is the same (3mV). Negative values are clamped to zero. A high value of the mean adaptation variable for a column makes it dormant and does not allow it to activate until it recovers from the dormant mode. A low value, however, makes it excitable and ready to propagate the activity.
https://doi.org/10.1371/journal.pcbi.1006216.g006 excitable mode. Note that the kernel with shorter time constant (γ 1 (t)) is not suitable for this kind of manipulation, because the value of the kernel converges rapidly to zero.
For simulating the stimulus-evoked response we initialize all assemblies in the excitable mode. The initial value of γ 2 (t) of each neuron is randomly selected from a Gaussian distribution with zero-mean and σ = 3mV. All negative values are clipped to zero (Fig 6D (left)). After stimulation of the central assembly, it converges to the high point and produces a burst of spikes. As a result, neighboring assemblies receive some synaptic input. Although the value of this input current is low due to low inter-column connectivity, it suffices for the neighboring assemblies to pass the switch point and rapidly converge to the high point. This scenario repeats, and so the activation spreads over all assemblies. The assemblies transit to the dormant mode after a burst of activity, consistent with experimental data [30].
The simulation of activity circulation is more complicated and needs careful tuning of the initial values of the adaptation variables of the neurons. Figuratively speaking we carve a path for the activation by choosing suitable initial values for the firing thresholds (Fig 6D (right)). For neurons inside each column, we choose again initial values randomly from a Gaussian distribution with σ = 3mV. The distribution's mean is different in each column. The means are selected such that there is path of excitable assemblies from the top-left of the grid to the column below the center. The activity propagates only along the path of excitable assemblies, whereas dormant assemblies do not switch to the high point even when they receive synaptic current. Once the activity has passed through the initially excitable assemblies, these become dormant. At the same time, the assemblies that were initially dormant recover so that the activity continues its path. This phenomenon repeats and causes activity circulating of the whole grid. After several rounds the circulation terminates because of a shortcut problem (see below). Non-assembly and inhibitory neurons do not contribute to shape the circular activation pattern. Instead, when they receive the activation from the assembly of their own column, they show some depolarization and a few spikes.
One may think that pre-shaping an activation pattern by initialization is artificial. However, any distribution of initial values which make several neighboring assemblies excitable and leave others in the dormant mode has potentially such a "shaping" effect on the potential trajectory of the activity. Different patterns of initial values lead to different patterns of activity propagation observed in the cortex. In other words, we suggest that the great trial-to-trial variability of activity previously observed via voltage-sensitive dye imaging in the cortex [48,49] might be due to different initial values of neurons in each trial. Sofar, we have tuned our pattern to construct a specific path; but in order to investigate the role of initial values more systematically, we have simulated 1000 trials. For the initial values of each trial, we randomly chose a shuffled version of initial values shown in Fig 6D (right). Specifically, we shuffled the mean of initial values of each column, but used a random distribution with σ = 3mV around the desired mean for neurons inside each column. Fig 7 indicates the types of activity propagation and the durations that activity survives in the grid. In * 23% of trials the activity does not spread over all columns. In * 35% of trials the activity spreads over all columns without repetition. In the rest of trials, cortical columns form a path of activation and the activity circulates for one or several rounds.
Termination of activity. Even with a highly tuned path, the activity terminates after several rounds. In order to understand the reason for this, we compare the activity path in early rounds and late rounds (S3 Fig). In the early rounds, the activity orbits around the central assembly and passes through all marginal assemblies to eventually reach the starting point at the corner again. Since, in our simulations, the central assembly has been assigned a very high initial value for the neuronal firing thresholds, it does not synchronize with its neighbors and does not initially follow them. However, after a while the central assembly becomes excitable and eventually activates when its, say, right-sided neighbor becomes active. Thereby it creates a shortcut for the activity which reaches the left-sided assemblies before the recovery from the dormant mode is complete. Since these assemblies are not recovered yet, they cannot become active. Therefore, the activity terminates and does not circulate any further. If we remove the central assembly or prevent it from being excitable (e.g. by decreasing its network feedback or increasing its intra-column inhibition), the activity can circulate for a very long time (green dots in S4 Fig). This is also the reason why we used absorbing boundaries for the grid. In case of reflecting boundaries (i.e., the activation wave reflects when it reaches the boundary) assemblies in the dormant mode do not have enough time for recovery and cause termination of the wave. We consider the alternative of periodic boundary conditions (i.e., when the activation wave passes a boundary, it jumps to the other side of grid) as not realistic, because once activity passes the boundary of a cortical area in the cortex, it generates (weak) propagation into a neighboring area rather than returning to the same area. Hence, we consider absorbing boundary as a good compromise and a step toward a very large scale model of neocortex, but simulating a large scale model is beyond this study.
As we described above, in our model the circulation of activity terminates spontaneously. This is consistent with spontaneous barrel cortex dynamics observed in vivo [31]. More  Fig 6D (right) and shuffled the columns. In *23% of trials the activity ceases before reaching all columns. In *35% of trials we observed that the activity spreads instead of circulating. In the remaining trials, we observed circulation with different patterns and durations depending on the path carved by the choice of initial values. The durations longer than 1000ms are clamped to 1000ms for better visibility.
https://doi.org/10.1371/journal.pcbi.1006216.g007 specifically, every time the activity vanishes in the model, there is need for stimulation of a corner assembly which is in the excitable mode. If such stimulation is provided, the activity can circulate again for several rounds. Such a stimulus could potentially be provided by other cortical areas adjacent to the barrel cortex by reverberations of activity in thalamo-cortical loops or by novel whisker stimulation.

Discussion
We have shown that a network of spiking neurons consisting of a chain of 11 assemblies can generate reliable temporally structured activity, that extends around one second, similar to excitable media on discrete structures [50]. While two well-known models of activity propagation, the synfire chain and rate propagation, require a feedforward network structure, our excitation chain works with a bidirectional activity pattern, similar to excitable media [22][23][24]50]. This is advantageous because there is no direct experimental confirmation for the existence of systematic feedforward connectivity patterns in the brain [1][2][3]7]. Moreover, the biophysical prerequisites of the excitation chain model, namely neuronal clustering, spike-frequency adaptation, recurrent connectivity and short-range inhibitory connections, are in principle consistent with experimental findings [29,32,38,42]. Recent experimental evidence suggests a temporal evolution of activity over 500-1500ms [51,52].
The excitation chain carries the signal by a wave of synchronous activity. Seen from this perspective, it is more similar to the synfire chain than to the fluctuating rate propagation which relies on asynchronous activity and rate coding [6]. However, the excitation chain is fundamentally different from the synfire chain and rate propagation in its ability to adjust the propagation speed without modifying the synaptic delay or synaptic time constant. While the propagation speed is nearly constant in a synfire chain [53], and is very close to the synaptic time scale in fluctuating rate propagation [7] and two-dimensional systems of neurons [54], we can slow down and adjust the speed in the excitation chain by changing the strength of inter-assembly synapses. Hence, we suggest that the excitation chain can be used for population coding across slow time scales [55,56]. It is also possible that changes in synaptic weights are caused by neuro-modulators or various forms of activity-dependent synaptic plasticity. Suppose that a part of a neuronal structure which produces a complex cognitive behavior is an excitation chain, such that each assembly is responsible for performing one primitive of the behavior. Repeating the behavior, the assemblies become active one after each other in a systematic and reliable manner. In this case, a Hebbian learning rule will strengthen the synapses between excitatory assemblies and thus will work towards an increase of the execution speed of the behavior as we observed in Figs 1 and 2. This may explain why practicing sequential or rhythmic movements, such as in playing a musical instrument, can increase the speed at which the movement can be performed.
The bistability of neuronal assemblies is a key component of the excitation chain model. While the mechanism of assembly activation (jumping to the high point) is based on network feedback and is independent from the neuron model, the return of the assembly to the low point (in the dormant mode) needs an element of fatigue. In our model, we use spikefrequency adaptation for this purpose. However, replacing adaption by short-term depression [35] of intra-assembly synapses yields similar results. In a previous abstract [57], we removed spike-frequency adaptation and built a chain using leaky integrate-and-fire neurons while the synapses within each assembly expressed short-term depression. With strong connections inside assemblies, and relatively weak connections between them, the chain could propagate the activation both forward and backward. Similar to what is presented here, the propagation speed was tunable by changing the inter-assembly maximum synaptic weights. In a related model, a feedforward chain of neural assemblies with strong recurrent synapses which express short-term depression will propagate activity with propagation latency over 10 assemblies between 90ms and 200ms [9]. In contrast to this earlier model, we use symmetric connectivity so that propagation has no preferred direction. Another related work [58] studied the influence of short-term synaptic dynamics on switching activity patterns in networks that contain attractors.
One may wonder about the plausibility of neural assemblies in cortical networks. Applying synaptic plasticity rules in a randomly connected network in a balanced state [59,60] leads to subgroups of neuron with strong recurrent connections [61]. A recent study used spike-timing dependent plasticity rules acting on excitatory to excitatory [62] and inhibitory to excitatory connections [63] as well as a synaptic scaling rule [64]. It was found that neurons which initially fire at higher firing rates compared to other neurons form subgroups with strong inward synapses and larger mean of outward synaptic weights after applying plasticity rules [61].
One way of extending our work is embedding the excitation chain in a bigger background network [5,16,17]. Then, excitatory assemblies and inhibitory populations receive synaptic input from the neurons of the larger background network. In case of a balanced background network, we would expect that the chain will exhibit results similar to those shown here, except that due to variance of synaptic input, the activation times will show more jitter compared to the result we obtained (Figs 1D, 1E and 2C). Constructing a balanced network for adapting neurons will be possible using recent mean field methods [65] and is left for future work. It would also be interesting to further explore bidirectional synaptic connections between excitatory and additional inhibitory neurons within the same group. Our expectation is that, if parameters are tuned to get a balanced state, the peak firing rate of neurons inside the assembly would be significantly reduced without a substantial change in the speed of signal propagation across the chain; cf. In the grid model, assemblies are responsible for transmitting the activity between barrel columns. Removing the inhibitory and non-assembly neurons from the model, the grid of assemblies is able to propagate and circulate the activity on its own. Therefore we consider the grid of assemblies as the skeleton of the model, while the other neurons just follow the activity wave. This, however, is not meant to imply that inhibitory and non-assembly neurons do not play an important role in the cortex. Whereas in our model, we only consider the anesthetized state, information processing in the awake cortex very likely involves more than just the assembly neurons and requires the contribution of the other neurons.
In the literature, there exist several competition-based models for reproducing cortical trial-to-trial variability [66][67][68][69] and for modeling working memory using continuous attractors [70,71]. In these models, individual neurons or neural assemblies (also called neural clusters) try to become active and to suppress others by direct inhibitory connections or global inhibition. However, in the two models considered here (the chain and grid of assemblies) the situation is completely different. Instead of competition, the assemblies cooperate with each other to propagate activity signals through the tissue. Hence, in this type of model there is no need for global inhibition.
Similar to bump attractor models proposed for implementing working memory [70,71] and attractor maps used for encoding spatial location in hippocampus [72][73][74], our model relies on strong recurrent excitatory connections. However, while these encoding-models are designed to stabilize neuronal activity in time in order to provide a stationary memory of a continuous variable, our model's aim is to propagate the activity with a specified speed and path to relay sensory information or to perform slow behavioral tasks.
Activity propagation was studied previously using the Wilson-Cowan model of neuronal population dynamics [22,23] and the dynamics of neural fields [24,25,50,75,76]. However, it is often hard to link network parameters used in these abstract models to neuronal and synaptic parameters. In our approach, we start directly from neural parameters extracted from experiments [38,39] and explore the influence of synaptic weight or short-term plasticity [35]. Despite the richness of the biological parameter space, the essential mathematical features of excitable media [22][23][24][25] are still apparent in the models examined here. In the literature, there exist several studies [77,78] that observed traveling activity waves in a balanced network with distance dependent connectivity structure.
The multicolumn model presented in this work is able to produce dynamics similar to the spontaneous regime of cortical activity [31]. However, it is different from resting-state networks [79][80][81][82][83] designed for explaining the synchrony of different brain regions. In these models, each brain region is typically considered as a nonlinear oscillator. The coupling strength between oscillators [79] as well as synaptic transmission delay and noise [80,81] are tuned such that different brain regions show synchrony similar to observed data. While such models are very successful in representing resting-state activity, their focus on macroscopic network structure may limit the range of dynamics produced by the models. In contrast, in our model we focus on a particular element of neuronal circuit-level connectivity, namely local assemblies of neurons, which may be contained in a small part of each cortical column. While such assemblies may make only a small contribution to the average activity, macroscopic propagation of activity across barrel columns or brain regions may depend on such assemblies. In systems with excitable elements such as ours, the type of dynamics strongly depends on initial conditions of the neurons and nonlinearly on the inputs [84]. Therefore embedded excitable assemblies provide nolinear mesoscopic processing characteristics to the circuits that may easily be overlooked in resting-state or other macroscopic models relying on averaged measures of connectivity between areas.
To conclude we would like to highlight the predictive aspect of our work. So far, neural assemblies have been used for explaining working memory [70,71], cortical trial-to-trial variability [66][67][68][69] and slow oscillations in the cortex [34]. We suggest that beside these roles, neural assemblies are also responsible for activity with adjustable speed in the cortex [52]. Consider a complex behavioral task which includes a sequence of several subtasks. If activation of each assembly of the chain initiates a subtask, the whole chain is able to perform the complex tasks with desirable speed. In case of repeating the tasks, Hebbian learning strengthens interassembly synapses and faster propagation. Therefore, the cortical network is able to perform the task faster after practicing the task. Because of symmetric connectivity, we could explain bidirectional activity propagation observed in the experiments [32]. Furthermore, our work predicts that inhibitory neurons reduce the propagation speed. In a multi-array recording [32] with GABA receptor antagonists near the electrodes of array, we therefore predict an increase in the propagation speed compared to the control condition.

Neuron model and population parameters
Neuronal parameters used in the simulations are reported in Table 1. Table 2 summarizes the network parameters.
As our neuron model we use a current-based generalized integrate-and-fire (GIF) model which implements spike-frequency adaptation using a spike-triggered current and a moving firing-threshold mechanisms [38]. The dynamics of the neuron's sub-threshold membrane potential (V(t)) is described by: where parameters C, g L and E L are the passive parameters of the neuron. I(t) is the synaptic input and η(t) is the shape of the spike-triggered current caused by spikes of the neuron itself at timest j . After each spike emission, the membrane potential is reset to V reset , integration of Eq 1 restarts and the neuron goes through an absolute refractory period of duration τ ref .
Spikes are produced stochastically (similar to an inhomogeneous Poisson process) with the firing intensity: where λ 0 is the stochastic intensity at the firing threshold V T , and ΔV is a constant which defines the level of stochasticity. The threshold V T follows the dynamic: where V Ã T is constant and γ(t) describes the time course of the threshold after a spike emission. In Eq (1), the synaptic input I i (t) received by neuron i is determined by the spikes of synaptically connected neurons: where w ij is the weight of the synapse connecting neuron j to neuron i, and aðtÞ ¼ e À ðtÀ DÞ=t syn for t ! Δ is the post-synaptic current (PSC) shape. The synaptic transmission delay (Δ) in all simulation is 1ms. S j ðtÞ ¼ P k dðt À t k j Þ is the spike train of neuron j, δ denotes the Dirac   ( ÃÃ ) This value is 13.05e −t/45.0ms for the simulation of the multicolumn model (Fig 6). δ-function and t k j is the k th spike of neuron j. The synaptic weight w ij indicates the PSC amplitude. Given the neuronal parameters, one can relate the PSC to the post-synaptic potential (PSP). We report values both of PSC and PSP amplitudes used in our simulations.
We ran simulations using the Brian simulator [85] for simulating the chain and NEST [86] for simulating the grid.

Transient stimulus
In order to initiate the activity in the chain, we stimulate one assembly of the chain (for example the first assembly, the last assembly or an assembly in the middle of the chain). For stimulating an assembly, we connect each of its neurons to 25 Poisson neurons which fire with the rate of 5Hz each for 25ms. The synaptic weight between the Poisson neurons and the assembly neurons is 0.18nA.  (Fig 2) ÃÃ Only applicable for the chain containing two types of inhibitory neurons (Fig 4) Simulation of multicolumn barrel cortex model (Fig 6) Size of exc. assem. 70 Size of exc. non-assem.

Rate-current relations
Here we present some theory necessary to explain the dynamics of excitable assemblies. The dynamics of a neuronal population can be described by two equations relating the firing rate averaged over all neurons and the mean of the synaptic input received by them. The first relation is called the neuron's gain function. Injecting a weakly fluctuating current I syn into a neuron produces an average firing rate of where hI syn i and σ I are the average and the standard deviation of the synaptic current over time, respectively, and g is the gain function. Although there are ways to compute the firing rate of adaptive integrate-and-fire neuron models in closed-form [87,88] or by using a selfconsistent numerical approach [89][90][91], there is no straightforward analytical solution for computing the gain function of the GIF model that we use here. We obtain the gain function (5) by numerical simulation [34]. For the simulations to determine the gain function numerically the injected current is given by where ξ(t) is white noise with mean hξ(t)i = 0 and covariance hξ(t)ξ(t 0 )i = δ(t − t 0 ), α(t) is the shape of the PSC defined above and q 2 ¼ R 1 0 a 2 ðtÞdt. Depending on the duration of injection, the neuron goes into different adaptation states. By injecting the current (6) for a short episode of 10ms, we can estimate the firing rate in the non-adapted state. In case of a longer stimulation period, we can divide the time into intervals of 10ms and extract the rate-current relation in the different, progressively more adapted states. This method has been used to obtain the gain functions displayed in Fig 3A and 3B.
The network activity gives rise to the second relation between the average firing rate and the average synaptic current. The synaptic input of an arbitrary neuron i is described by: where w ij is the weight of the synapse connecting neuron j to neuron i and S j (t) is the spike train of neuron j. The sum runs over all other neurons j in the assembly. Averaging both sides over time and input neurons gives the average input current: hI syn i = Npqwr, where N is the number of neurons inside the population, p is the connection probability between neurons, w is the synaptic weight and q is the total charge of one PSC pulse: q ¼ R 1 0 aðtÞdt. Rearranging this equation yields: We refer to the denominator of Eq 8 as the network feedback coefficient (C fb ) of the population [34]: We use these two relations for analysis of the behavior of excitatory assemblies in the "Results" section.

Analytical approach for obtaining the propagation speed of the excitation chain
Our aim is to estimate the difference between activation times of two consecutive assemblies (which is a measure of propagation speed) using an analytical approach. We assume that the values of the parameters are given and the time course r(t) of the population rate signal of an arbitrary assembly in the chain is known. Furthermore, we assume that the populations are silent before activation (r(t) = 0 for t < t a ). Suppose that we have two excitatory assemblies, assembly 1 and assembly 2. We refer to their activation time as " t 1 and " t 2 respectively. The aim of the calculation is to find the difference in activation times x ¼ " t 2 À " t 1 . Recall that the activation time has been defined above as the expected time at which all assembly neurons have spiked once. Assuming independent Poisson firing of the assembly with rate r(t), we can find " t 1 as the moment when the expected spike count reaches 1, nð0; " t 1 Þ=N ¼ 1, where n(a; b) is the number of spikes in time interval [a, b) and N is the number of neuron in the assembly. Inserting the shape of assembly population rate we obtain which we can solve for " t 1 . Note, however, that the Poisson assumption made here does not ensure that no neuron fires more than one spike before " t 1 . Therefore Eq 10 yields an approximate value of " t 1 .
Finding the value of " t 2 is more complicated. After activation of assembly 1, its neurons send synaptic input to assembly2. The average input received by each assembly2 neuron can be computed: where p exc and w exc are inter-assemblies connection probability and synaptic weight respectively, Ã denotes the convolution operator and Δ is synaptic transmission delay. However, this is not the only synaptic input received by neurons in assembly2. Even before " t 2 , several neurons of assembly2 may already fire spikes (due to random fluctuations) and send some feedback current to other neurons. This averaged current can be computed similarly: where p self and w self are intra-assemblies connection probability and synaptic weight respectively, and x is the difference of activation times. Note that the shape of r(t) is assumed to be the same for all assemblies of the chain. However, we must be careful about the timing of each current. Suppose that at the time t = 0 assembly 1 starts to fire, therefore assembly2 receives a first input at the time t = Δ. After a while, assembly2 starts to fire. The difference between these two starting times is denoted by x. Therefore the feedback current is received by assembly2 neurons at the time t = Δ + x.
The total synaptic input received by neurons of assembly2 is the summation of I syn_fwd and I syn_self . Consequently, we can write the total input received by assembly2 neurons as I syn ðtÞ ¼ I syn fwd ðtÞ þ I syn self ðtÞ ð13Þ Now we can calculate the subthreshold membrane potential of neurons in assembly2 by solving Eq 1: Note that since we want to calculate the time of first spikes of neurons, we can neglect the spike-triggered current and the moving firing threshold. The firing intensity is then given as a function of V(t) by Eq (2), except that V T ðtÞ ¼ V Ã T because we assumed that neurons are not adapted and all γ(t) equal zero.
Then using the distribution of first spikes P(t), we are able to calculate the average time of first spikes of assembly2: Finally, using Eqs 10 and 16, we can calculate the time difference between activation times: Note, however, that " t 2 in Eq (17) depends on the value of x through Eq (12). Therefore, we have formed a self-consistent equation for x. We feed this value in Eq 12 and get it back in Eq 17. If the output value of x equals its input value, we found the proper value. Using a simple search, we are able to find this value numerically.
We apply this approach for our chain and calculate the value of x for different values of w exc ; these results are presented in Fig 3D. Note that we need to obtain the shape of the assembly population rate r(t) by neural simulation beforehand. However, the same shape r(t) can be used for all values of w exc , because for p exc w exc ( p self w self the time course of the initial rise is dominated by the self-feedback (see Fig 3B).

Short-term plasticity
We use short-term plasticity [35,36] in one series of our simulations (Fig 1F). This synaptic model supposes that each synapse has a certain amount of resources (such as neurotransmitter packed in vesicles ready for release) denoted by x, with dynamics where U (jump of release fraction), τ rec (recovery time constant) and τ facil (facilitation time constant) are three parameters of the model. u (release fraction) and x are the two variables of the system. Whenever a neuron fires a spike (t f denotes the firing time), it produces a PSC with amplitude of uxw (w is the synaptic weight). Then, the amount of resource is decreased by ux and the release fraction is increased by U(1 − u). In our simulations, we either used facilitation or depression. We chose the values τ rec = 0.001ms and τ facil = 500ms for the facilitation case and τ rec = 800ms and τ facil = 0.001ms for the depression case. The value of U is different in each simulation. Note that for the depression case, we have to fix the amplitude of the first PSC regardless of value of U. We did that by adjusting the value of the synaptic weight.

Shuffling initial values on grid
We simulate different trials of the grid with shuffled initial values (Fig 7). The aim of shuffling is to show that many initial patterns will lead to activity circulation and not just the tuned pattern of Fig 6D (right). For the shuffling, we take mean values of γ 2 (0) (initial value of the moving threshold kernel with longer time constant) of each column in Fig 6D (right) and randomly assign these means to columns. Then, using these means, we randomly choose the value of γ 2 (0) for each neuron in the columns. The activity circulation in the grid terminates after several rounds because of a short circuit in central assemblies (C4!C3!C2). In early rounds the difference of activation time between C4 and C3 are shorter compared to the late rounds (red ellipses). In other words, C3 becomes active sooner than it is expected. Therefore, C3 is able to activate C2 (red arrow), while C2 is supposed to be activated by B2. This short circuit generates activity before assemblies recover from the dormant mode. Hence, other assemblies are not able to become active and the circulation ceases. The initial values are the same as shown in Fig 6D (right). Column A1 is stimulated at t = 100ms in order to start the circulation.