Different propagation speeds of recalled sequences in plastic spiking neural networks

Neural networks can generate spatiotemporal patterns of spike activity. Sequential activity learning and retrieval have been observed in many brain areas, and e.g. is crucial for coding of episodic memory in the hippocampus or generating temporal patterns during song production in birds. In a recent study, a sequential activity pattern was directly entrained onto the neural activity of the primary visual cortex (V1) of rats and subsequently successfully recalled by a local and transient trigger. It was observed that the speed of activity propagation in coordinates of the retinotopically organized neural tissue was constant during retrieval regardless how the speed of light stimulation sweeping across the visual field during training was varied. It is well known that spike-timing dependent plasticity (STDP) is a potential mechanism for embedding temporal sequences into neural network activity. How training and retrieval speeds relate to each other and how network and learning parameters influence retrieval speeds, however, is not well described. We here theoretically analyze sequential activity learning and retrieval in a recurrent neural network with realistic synaptic short-term dynamics and STDP. Testing multiple STDP rules, we confirm that sequence learning can be achieved by STDP. However, we found that a multiplicative nearest-neighbor (NN) weight update rule generated weight distributions and recall activities that best matched the experiments in V1. Using network simulations and mean-field analysis, we further investigated the learning mechanisms and the influence of network parameters on recall speeds. Our analysis suggests that a multiplicative STDP rule with dominant NN spike interaction might be implemented in V1 since recall speed was almost constant in an NMDA-dominant regime. Interestingly, in an AMPA-dominant regime, neural circuits might exhibit recall speeds that instead follow the change in stimulus speeds. This prediction could be tested in experiments.


Introduction
Spiking neural networks (SNNs) are important microscopic substrates for realizing the nodes in more abstract physiological networks governing the global interaction between different brain regions and other physiological subsystems [1,2]. Since functional interactions within a physiological network depend on the temporal dynamics and correlated activity between the nodes of the network [1], the understanding of learning and retrieval of temporal activity in SNNs is a crucial basis for investigating the state transition and experiencedependent changes in a physiological system. Sequential activity propagating across neural tissue is common in the brain and potentially underlies important functions. Spontaneous replay of spatiotemporal sequences has been observed for instance in the 2. Materials and methods

Spiking network
The SNN model is composed of N integrate-and-fire neurons connected along a linear chain. The network structure is illustrated in the figure 1. The ratio of the number of excitatory and inhibitory neurons N N : E I is set to 4:1. All N E neurons are distributed along a line, inhibitory neurons are inserted into the chain in regular intervals. Neurons i and j are connected with a probability P ij , which decreases with the distance between the two neurons. In more detail, the connection probability from neuron j to neuron i is given by = In the network model, we consider two kinds of synaptic plasticity: long-term plasticity and short-term plasticity. Long-term plasticity is governed by a STDP rule. STDP is a kind of Hebbian learning rule that is widely believed to be the basis of experimental observed forms of long-term synaptic modification [47,48]. We assume that STDP modifies all excitatory-to-excitatory synapses (denoted by 'e2e syn'). Other synapse types remain fixed [49,[51][52][53].
Short-term plasticity is included in form of STD, which is commonly observed in the axon terminals of presynaptic neurons in the visual system [55][56][57]. It has been used widely in network modeling studies [44,45,51,58,59].
In the following, we describe the network dynamics in mathematical detail. Definition and values of all network parameters are summarized in table 1.

Single neuron dynamics
The single neuron dynamics is described byintegrate-and-fire neuron model where V i is the ith neuron's membrane potential andV L the resting potential. When ⩾ V V i th , the membrane voltage V is reset toV re . τ m and R m are the time constant and the resistance of membrane potential, respectively.   [41], the retina of mice were stimulated repeatedly on a fixed path with a moving light spot (swiped from the starting point to the end point). Activity was recorded in neurons with receptive fields on the path (indicated by round circles). After multiple repeats, spontaneous activity recall (that is a sequence of activity from starting point the end point) was observed when only stimulated briefly at the starting point. (b) A chain-like network model used to simulate the recall experiment. The network consists of excitatory (green disks) and inhibitory neurons (red disks). Green and red Gaussian bumps represent distance dependent connection probabilities for excitatory and inhibitory connections, respectively. (c) Illustration of how the moving stimulus was applied in the stimulation. Stimulus was applied simultaneously to a group of neighboring neurons with groups size S sti , which shifted by one neuron each ΔT sti along the chain. Resulting stimulus speed is ΔT 1 sti . The stimulus is conveyed to each neuron by additional input of N sti independent Poisson spike trains. magnesium block causes the maximal conductance to increase nonlinearly with the membrane voltage [43]. Following [43], the currents resulting from synaptic interactions are given by where g AMPA , g GABA and g NMDA are the maximal conductances of AMPA, GABA and NMDA, respectively, and V E andV I are the excitatory and inhibitory reversal potentials, respectively. It is = C 1 ij if there exists a connection from neuron j to i, otherwise 0. W ij is the synaptic weight from neuron j to i, which is subject to synaptic plasticity (see below).

Synaptic dynamics
The variables S ij X ( = X NMDA, AMPA, GABA) are the gating variables corresponding to the synaptic currents from neuron j to i. If no spike arrived, these gating variables decay exponentially with time constant τ X [43]. Upon spike arrival, their value increases by an amount that is regulated by the recent spike history based on the dynamics of the STD. It is where t j k denotes the time of the kth spike of neuron j.
is the transmission delay from neuron j to i. Note that the transmission delay is distance dependent.
The variable x j (t) in equations (7) and (10) is the STD variable of the excitatory neuron j, and can be interpreted as the total amount of immediately available neurotransmitter [56][57][58]. STD in inhibitory neurons is neglected here, because inhibitory neurons commonly show less depression than excitatory neurons [57]. ∈ U (0, 1) is the utilization of synaptic efficacy when a spike arrives [58]. The dynamics of x j (t) is given by [45] where τ D is the depression time constant.

STDP rules
The synaptic weight variable W ij (t) is subject to an STDP learning rule only if neuron i and j are both excitatory neurons and is initially set to 0.5 for all connections. STDP can cause the weight variable to change over time in the range from 0 to is the effective maximal synaptic conductance at time t for receptor X.
We here discuss two types of STDP rules. One is the additive STDP rule (A-STDP) [50][51][52][53] defined by  [52,54,60,61], where + A and − A are learning rates of facilitation and depression, respectively, and τ + and τ − are the corresponding time constants. Δt ij is the time difference between a spike time of the postsynaptic neuron i and the arrival time of a spike of the presynaptic neuron j.
We further divide the STDP rules into two flavors: all-to-all (AA; [52,53]) and (semi)-NN ( [50-52, , ]). In a STDP-AA rule, all existing spike time intervals Δt ij are used to compute the weight update. Alternatively, in a STDP-NN rule, only NN contributions count. More specifically, if the post-synaptic neuron spikes at time t only the interval to the last pre-synaptic spike before t is considered as a positive contribution in the weight update. Analogously, if the pre-synaptic neuron spikes at t, only the time interval Δt ij to the last post-synaptic spike before t counts as negative contribution (semi-NN STDP rule).
In generally, we assume that the STDP variables of AMPA and NMDA have the same plasticity rate [51], and are therefore governed by the same weight variable W ij (t) (see equations (4) and (6)). As shown in the main text, we found that M-STDP-NN best corresponded to the experimental observations of [41]. Hence, most of the simulations presented in the paper were done using M-STDP-NN. In the simulations, we adopted an online implementation to realize the STDP weight updates [52].

Input
In a local neural systems, neurons are always bombarded with background activity from other brain areas. One way to model this external background noise input is by assuming that each neuron in the circuit received additional excitatory and inhibitory inputs from a number of independent Poisson spike trains. To save computational time, however, we simplified this assuming a Gaussian noise process as additional current input, which can be seen as a mean-field approach for the Poisson background spike inputs. We found that either way of modeling background activity did not alter our results. Thus we assumed an additional Gaussian noise current input of the form μ δ ξ = + I t ( ) bg bg bg . In the training process, the moving stimulus was defined as follows. Independent excitatory inputs composed of N sti independent Poisson spike trains with rate ν sti (modeled as external AMPA currents analogous to equations (4) and (7)) were simultaneously loaded on each neuron in a group of nearby neurons (stimulus size S sti ) for a fixed duration ΔT sti (see figure 1(c) for an illustration). After this time duration, the stimulus slides to the next new group of neurons which is overlapping with the previous one. In the simulation the overlap was of size − S 1 sti , so that the stimulus simply shifts by one neuron along the chain. The above two steps were successively repeatedly from the starting point until the end point was reached. Hence, the stimulus speed can be calculated by ΔT 1 sti . In the retrieval process, the stimulus cue was set to the same size as the moving stimulus during training and lasted for 50 ms at the starting point.

Rate network model
Using the mean-field approach, we approximate the spiking model with a rate model. By deriving a simple relation of the predicted recall speed, we attempt to infer the recall speeds that would be expected for a given learned weight distribution in the spiking network model.
We consider two variables for each of the N rate neurons: the inner state u(t) and output firing rate r(t). For a neuron i their dynamics is given by where τ u is a time constant, τ d a constant transmission delay, I bg i the constant background current input, R the maximum rate of a neuron, and f(u) an activation function. For the numerical analysis, we take the sigmoid function κ = f u u ( ) tanh ( ), and simplify f(u) further to the Heaviside step-function ( = f u ( ) 1if > u 0, otherwise 0) in the theoretical analysis below. ω ij is the connection strength from neuron j to i, which is equal to the product of the connection probability and the synaptic weights of the trained SNN in the mean-field approximation. α is a normalization parameter to scale the total effect of all recurrent inputs on the neurons' state. Since a constant current input would have the equivalent effect of decreasing the threshold u th by a value I bg i , we set = I 0 i bg for simplicity. The function h ( · ) comprises the synaptic dynamics and will be discussed below.
This rate model greatly simplifies the original SNNs and thus allows for some theoretical results.

Theoretical analysis of the recall speed
Assume that a stable wave of spike onsets is traveling along the chain of neurons. If all connections are translation invariant ω ω = ij j , the time the activity needs from the onset time of the previous neuron to the next, τ r , is directly related to the speed of propagation We assume for simplicity that a neuron receives enough input to stay above threshold once it has reached it. Because of the delay, it is τ τ τ = + r d d y n , where τ dyn is an unknown time to be determined from the dynamics using a self-consistency method.
Assume the neuron u 0 firstly crosses the threshold at time t. Because each neuron has the same connectivity, in case of a stable propagation, all r j (t) time courses will look the same but only shifted in time by τ j r , i.e.
Since the wave here travels from negative to positive j and all neurons > j 0 do not contribute in the input at time of the threshold of neuron j = 0 (initial state is zero), we can write equation (14) for times < t t: Let us first consider that there is no synaptic dynamics, e.g. = h Id, we will add the dynamics below. Integrating this linear ODE until = t t we find The integral bounds were adjusted because the rate of each neuron is zero before reaching its threshold (that is for τ < + t t j r ). Shifting the integral bounds, we find If we assume now a Heaviside activation function for simplicity, it follows Since the weights ω i rapidly fall-off with distance (like a Gaussian, see below) and the time constant of each rate neuron (in the order of 10 ms) is an order of magnitude larger than the transmission time (order of 1ms), we can approximate the exponential function with a first order Taylor approximation, , and then solve for τ r , and find with an overall scaling factor. One notes from equation (20) that the larger the 'forward' weights, the higher the speed. All synaptic forward weights contribute to the speed but are weighted with distance.

Synaptic dynamics in the rate model
We approximate the synaptic dynamics h from the full spiking model (see equations (7)-(10)). Incorporating both NMDA and AMPA terms, we define AMPA AMPA D Note that we here have neglected the nonlinearity of the NMDA dynamics for simplicity, and set β = 1 AMPA , and β < 1 NMDA (see equation (9)). Since we integrate only from the onset of r(t) (see equation (18)) and assume that the rate function is a Heaviside function (times R, see equation (15)), we only have to integrate the equation from 0 to t with constant = r t R ( ) . In this case, one finds for ρ s t is the effective time constant of the synaptic depression. If we combine equations (25) and (18) and integrate according to equation (18) with the abbreviation Γ ρ given by is the Laplace transform of the weights in space. The equation (26) can be used to test how the speed in the rate model compares to the SNN. For that we assume that where the network weights ν j can be either fitted according to ν = + a bj tanh ( ) with parameters a and b to the weights distribution of full network after training or directly taking from the full network.

Results
We stimulated a SNN repeatedly with unidirectional moving stimuli. The network structure is shown in figure 1 and mathematical details are described in the materials and methods section. Briefly, the network consisted of a chain of interconnected spiking excitatory as well as inhibitory neurons with integrate-and-fire dynamics and conductance-based synapses. To resemble the experiments in the visual cortex [41], neurons were assumed to have equally spaced topographic receptive fields positioned on the stimulus trajectory in the visual field (see figure 1). Connection probability between neurons was distance dependent, so that nearby neurons were more likely connected than neurons further apart. Excitatory and inhibitory synapse dynamics were governed by fast and slow excitatory synaptic conductances (AMPA and NMDA) [43,51], and a fast inhibitory current (GABA), respectively. Additionally, all excitatory synaptic strengths transiently weakened on a short time scale (several hundreds of milliseconds) in an activity-dependent manner (STD) [45,51,55]. Long-term plasticity of excitatory-to-excitatory (e2e) synapses was subject to STDP [47,[49][50][51][52][53]. The STDP variable W ij controlled the amount of the long-term effects on AMPA and NMDA conductances for a synapse from neuron j to neuron i and could vary from zero toW max . The relative strength of the AMPA and NMDA conductances could be adjusted by setting the scale factors g AMPA and g NMDA (maximal conductances), respectively (see material and methods for details).
In the following sections, we first show that the network could be trained to recall the propagating activity when only a transient cue was given at its start point and describe two qualitatively different kinds of activity recall behaviors. Then we show that continued training resulted in stable weight distributions, which were stable at different levels in different parameter regimes and analyze the underlying mechanism. Finally, we show with a simplified rate model that the different recall speeds are indeed caused by the differently shaped weight distributions.

Training of temporal sequences
Before learning, all parameters in the network were set to biologically reasonable values (see table 1). STDP weights W ij were set to half of the maximum value and we used a multiplicative STDP rule with NN spike-interaction (M-STDP-NN) for long-term spike-time dependent modifications (for details see material and methods).
Initially, the network could not successfully propagate activity from one end to another when triggered by a transient cue at one end (see figure 2(a) for an illustration of the training and recall procedure). However, after repeated trials of swiping an input stimulus with constant speed and fixed size along the neurons (indicated with the red-shaded line in the figure 2(a)), synaptic weights were systematically and progressively modified (compare to the weight distribution for figure 2(b)). In particular, weights connecting neurons in the direction of stimulus movement ('feed-forward' or FF) were differentially enhanced, whereas weights in the opposite direction ('feedback' or FB) were reduced in strength (see figure 2(b)). This training-induced symmetry breaking of the weight structure allowed the activity to propagate along the motion path when a transient cue was given at the starting point after 100 repeats of training. Hence, the sequence recall was successfully trained. In the following, we defined a successful recall by whether the activation propagation reached the end of the neural sequence when only the start was transiently stimulated.

Fast and slow activity recall
Having established that spatiotemporal sequences could be entrained in our network using STDP, we next asked what parameters were important for successful recall. We observed that the size of the net current to nearby neurons, that is the total current with inhibitory currents subtracted, was very important to achieve successful recall propagation. Indeed, if the maximum conductance of AMPA g AMPA was high enough, activity could propagate from the very beginning even without training. We thus adjusted this net current by varying g AMPA (while the maximum conductance of NMDA g NMDA and other parameters were fixed). In general, mean recall speeds markedly increased with increased levels of g AMPA (see figure 3(a) inset). After training using the M-STDP-NN learning rule with a small learning rate = = + − A A 0.025, the mean recall speed reached rapidly a stable level for a sufficient number of repeats (see figure 3(a)). We defined the recall speed as the slope of a linear regression of the first spike times along the chain of neurons. The mean recall speed was calculated over all successful recalls (18 cue presentation). When the sequential activity reached the last 20% neurons of the network chain when triggered by a transient cue at the starting point the recall was deemed successful.
We found that two qualitative different recall behaviors existed when varying g AMPA (see figures 3(b) and (c)). In case of high AMPA ( = g 0.48 AMPA ) recall after training was fast and immediate, seemingly depending only on a couple of fast onset spikes per neuron and often even faster than the stimulus. Indeed, the average time difference of the first spikes of neighboring neurons during recall was only 0.25 ms and thus in the order of the decay time of the AMPA conductance (2 ms). On the other hand, in case of low AMPA ( = g 0.30 AMPA ), although sequence recall was equally successfully trained, the speed of the recalled propagation was often slow compared to the stimulus speed and seemed to rely more on an enhanced firing-rate propagation rather than single spike transmission ( figure 3(c)). As the ratio of (fast) AMPA and (slow) NMDA was different in these two cases, we call the latter slow NMDA-dominant and the former fast AMPA-dominant sequence recall.

Weight distributions converge to different stable states
The observed recall speed likely relate to the synaptic weight modifications induced by STDP. To explore this, we recorded the dynamical changes of the synaptic weight structure with increasing number of training trials (figures 3(d) and (e)). Analogous to the dynamical changes in recall speed with increasing trial number, the averaged feed-forward weights initially increased (while feedback weights symmetrically decrease) and rapidly reached a stable plateau. Moreover, the dynamics of the effective weights, that is, the maximal conductance g AMPA times the average STDP weights W ij , corresponded well to the increase in recall speeds (compare figure 3(e) with figure 3(a)): the larger g AMPA , the faster recall speed.
We found that increasing the stimulus size also affected the final value of the recall speed. Increasing the stimulus size resulted in a slower recall behavior(see figure 3(f)). Note that in the experimental observations [41], a stimulus size of about1 6 the length of the trajectory was used, which corresponded to a spatial extend of about 36 neurons in our network. In the following, we thus set = S 36 sti to allow for better comparison with the experimental study.

Different activity recall behaviors as a function of stimulus speed
Since experiments showed successful recall across a range of stimulation speeds [41], we further tested whether varying the stimulation speed would have an effect of the recall speed and weight dynamics in the network model. Indeed, both the converged mean feed-forward weight values as well as the recall speed depended on the stimulation speed (see figure 4). However, the degree of how the recall speed dependence on stimulus speed varied strongly with the level of g AMPA .
In particular, we found that in case of a NMDA-dominant slow recall, variation in the stimulus speed did only slightly affect the speed of recall. In fact, the recall speed was almost independent from the stimulus speed (see figure 4(a)): when stimulus speed was increased by a factor of 10 the recall speed only increased about 1.5 times. This is similar to the observation in the experiments [41]. On the other hand, in case of AMPA dominant sequence propagation, recall speed followed closely the stimulus speed: the ratio of both speeds was nearly 1.
Thus, our network simulation exhibited two different regimes, stimulus speed dependent or (almost) independent, regulated by the level of AMPA in relation to the level of NMDA conductances. These different phenomena are summarized in table 2.

Mechanics of learning different weight distributions under different stimuli speeds and g AMPA
While recall speeds did not change significantly for NMDA dominated propagation, the weight structure nevertheless changed profoundly with the stimulus speed in a similar manner as in the AMPA dominated case ( figure 4(c)). In particular, weights of connections between nearby neurons were increased for fast stimuli speeds. In the following, we explore the possible underlying mechanisms for the observed learning induced weight modifications.
Naturally, the higher an excitatory feed-forward weight is, the higher is the excitatory drive onto the following neurons. The higher drive increases the likelihood that the next neuron fires thus takes part in a propagating wave of activity. Since a forward connection lies in the direction of stimulus movement, the postsynaptic neuron will receive the stimulus pulse slightly later than the pre-synaptic neuron. Thus, if both neurons fire, the order of spikes is positive under the STDP rule and will result in a strengthening of the forward connections. This mechanism leads to the observed symmetry breaking between the strength of forward and backward weights (see figure 2(b)).
However, why do the forward weights do not keep growing with increasing trial number but instead remain stable at a certain level well below the maximal weightW max (see figure 3(e))? There are mainly two reasons. First, the multiplicative STDP rule automatically reduces the size of the weight update if a weight is already large (see material and methods), thus the facilitation effect decreases over time. Second, during stimulation of the network as well as during replay, neurons will fire multiple spikes, so that some pre-post spike pairs occur in the 'correct' order and cause enhancement of forward connections, whereas some will be reversed and thus cause depression. We show in the following that the overall amount of weight depression and potentiation depends on the stimulus speed as well as the level of g AMPA .
To explore the amount of facilitation versus depression in recurrent networks, we identified five different types of spike-pairs for any two neurons with fixed distance during training (see figure 5(a)). Our STDP rule only results in a weight update for the latest post and pre-synaptic spike relative to a pre-or post-synaptic spike, respectively (STDP-NN). If one selects a fixed distance of two neurons in the chain, one can categorize types of spike pairs happening in the overlapping stimulation periods and non-stimulus periods and all combinations between both (see figure 5(a) for the definition). These types of spike pairs will have different effects on the weight structure during training. For instance, because the stimulus is swiped over the chain of neurons, forward connections will be particularly strongly potentiated when the stimulus induction in the previous neuron has just ended (and thus the spike rate drops) while the next neuron is still stimulated (green color in figure 5(a)).
To understand why different weight structures are learned for different stimuli speeds and g AMPA values, it is essential to note that the contributions of these types of spike pairs changes with stimulation speeds and g AMPA Figure 4. Relationship between network recall and training speeds using the M-STDP-NN learning rule. (a) Recall speeds became stronger dependent on the stimulus speed during training if the strength of the AMPA conductances were increased. Note that small g AMPA levels result in an almost stimulus-speed independent recall similar to experimental observation in the visual cortex [41]. (b) Average effective forward weights after training (AMPA (left) NMDA (right)) increase linearly with increasing stimulus speed. Note that the slope of the weight increase does not follow the slope of the increase in recall speeds (panel (a)), and the mean forward weight becomes less changed for lower g AMPA when increasingV Sti (panel (b)). (c) Average incoming weight structure (upper panel) and connection probability (lower panel) with relative neuronal distance. Forward weights (negative relative distances) are strengthened, whereas backward weights are weakened. Note that the weights structure changes with stimulus speeds as well as with different g AMPA levels.  values. In particular, if the stimulation speed is reduced, the length of the stimulation period for each neuron is increased (because the light spot is swiped more slowly) and thus the stimulus drive is stronger (see mean rates in figure 5(b)) leading to a higher stimulus-induced weight strengthening. Although all types of contributions changes over the time of learning due to the characteristics of the multiplicative STDP rule (as shown in figure 5(d)), we can integrate the contributions from the first trial to the trial when weights have just converged ( figure 5(e)). One can see that for increasing stimulation speeds, the amount of stimulus-induced strengthening decreases (figure 5(e) right-up, blue line). However, in case of a strong and slow moving stimulus, the recurrent network activity reacts faster than the moving stimulus in that forward neurons spike before the stimulus has arrived. This results in a high amount of negative spike pairs under the STDP rule (see figure 5( The effect of increasing the levels of g AMPA is mixed as well. The main reason why higher levels of g AMPA result in weaker weights, is because the sporadic activation of neurons at times unrelated to the sequential firing pattern becomes more pronounced, especially at high stimulus speeds (see mean rates in figure 5(b) and STD variable in figure 5(c)). This higher spikes rate in the network leads to larger amount of negative spike pairs Note that stimulation is prolonged for low stimulation speeds (stimulus arrives at the selected neuron at relative time 0). Thus rebounding network activity after the stimulus is higher for higher speeds than for slower speeds because synapses are more depressed by STD. (c) STD variable corresponding to the four cases in (b). (d) Added recurrent and stimulus period spike pairs contributions (as defined in (a)) on the weight update for each successive trail. Note that changes mostly reduce over trials because of the multiplicative STDP rule, which decreases the size of the update with increased synaptic weight. (e) Integrated contributions (from the first trial to the trial when weight distribution reached a stable state) showing the effect of different types of spike pairs for different stimuli speeds and g AMPA levels. Note that potentiation and depression both occur and their relative effects depend on the recurrent activity, on the stimulus speed, as well as on the g AMPA level. Parameters as in We conclude that the complex interplay between recurrently evoked activity unrelated to the sequential activity, the overall excitability of the neural system, and the changing stimulus drive causes the network to converge to different synaptic weight settings.

Rate model predicts increase in network recall speeds
We found that the spatial weight distributions could reach different stable states under different stimuli speeds (see figure 4(c)). To establish what recall speed is to be expected from a given weight distribution, we simplified the full SNN and derived a simple rate model (see material and methods). In the rate model, we show that, in agreement with the intuition, increasing the value of even single forward weight will increase the recall speeds.
To analyze in more detail the correspondence between the learned synaptic weight structure in the SNN with the recall speed, we took the learned weights of the full spiking network model (from figure 4) and used the rate model to estimate the predicted speed. In the theoretical investigation, we found a self-consistency description of the recall speed, which could be solved numerically for a given weight structure. However, since the exact correspondence of the input drive is unknown, we fitted free parameters to best match 6 recall speeds of the full network for each AMPA level individually. We found that the rate model correctly predicted an increase in recall speed. In particular, the rate model well explained the recall speeds for low AMPA levels. However, the steep increase of the recall speed seen for high g AMPA could not be well matched (see figure 6). This mismatch indicates that the spiking network model might rely on nonlinear effects of the NMDA conductance (not included in the rate model) and coincident spike timings (as opposed to rate based activations) for the recall of the sequence.

Different STDP learning rules
We found that the use of the M-STDP-NN rule resulted recall speeds and weight distributions that agree with experimental observations, while the use of alternative STDP rules not. In particular, weight distributions learned with additive STDP (A-STDP) [52,53] where unrealistically uniform when driven by Poisson noise(see figure 7(a)). This does not correspond to experimental data where commonly a Gaussian bumb-like distribution of synaptic weights is observed [60]. Meanwhile, when training the network with A-STDP rules, all forward weights were driven to nearly the maximum weight value and recall speeds are fast regardless of different stimulation speeds (not shown). However, the ability to induce sequence recalls for very low stimulation speeds does not agree with the experimental observations, where too slow stimulation velocities did not results in any sequence recall [41].
Moreover, we noted that the NN spike interaction pairing is equally important. If one uses instead all spike pairs (M-STDP-AA), very small stimulation speeds result in very strong weight values and very fast recall speeds (see figure 7(b)). This result is also incompatible with the experimental observations. We therefore conclude that the M-STDP-NN rule best agrees with experimental observations in the visual cortex. Figure 6. Fitting of the stimulation and recall speed relationship using the rate model. Solid lines: full spiking network, dashed lines: best fit of the rate network. In the fit, the trained weights of the full network were used. For each g AMPA value, free parameters of the rate model were fitted to the original recall speed. Free parameters were: α u th (the overall input), and the rate R. The fit to the recall speed is good for weak AMPA conductances, however, poor for strong AMPA conductances, suggesting that the nonlinear NMDA activation and coincident spiking effects dominate the activity recall in the high g AMPA case. Other parameters: τ =

Discussion
In this study, we investigated the ability of spiking neuron network models to learn and recall spatiotemporal sequences. We found that is possible to induce a sequential recall similar to experimental observations using STDP as a learning rule. Two modes of recall activity were observed in the network models. A fast, AMPA-dominant dependent recall and a NMDA-dominant slow recall. Increasing the stimulus speed less affected the NMDA-dominant recall and thus does better correspond to the experimental study in the visual cortex [41]. In the anesthetized state in the experiments [41], recall speed was usually slower than that of the stimulation speed and thus similar to the NMDA-dominant behavior in the spiking network model. However, in awake animals, the recall speed increased by about three times and therefore better corresponds with our AMPA-dominant case. It might therefore be possible that recurrent dynamical properties change between anesthetized and awake states.
In the AMPA-dominant state, we found that transmission was based on a fast propagation of initial spike volleys and the dependence of recall and stimulation speed was less well explained by a rate based model. In this fast spike-based transmission (as opposed to rate-based in NMDA-dominant sequence recall), recall speed could even be faster than the stimulus speed, as was also observed experimentally in the awake state in the primary visual cortex [41] as well as in hippocampus [3]. This fast transmission seems to rely on a few nearly coincident spike arrivals in a post-synaptic neuron, which evoked synaptic AMPA current strong enough to drive the neuron above threshold. In additional control simulations, we indeed found that this fast transmission solely relied on the AMPA dynamics: If the NMDA conductances were completely removed in the network, sequence recall was still possible and recall speed was very fast, similar to the AMPA-dominant case investigated in the results.
Our model further predicts that the AMPA-dominant fast recall is more dependent on the stimulus speed than in the slow NMDA-dominant recall. Unfortunately, the experimental study [41] did not explicitly investigate the stimulus-recall speed dependence in the awake state, so that we cannot compare this prediction to the published experiments. Interestingly, this high dependence on speed of the AMPA-dominant case could be used to imprint and memorize different sequence speeds. In apparent agreement, strong synaptic terminals are observed in the Mossy-fiber pathways in the hippocampus, which was shown to be important for temporal  NN). In case of M-STDP-AA, recall speeds become very fast for slow stimulation speeds. This is in contrast to experimental results where the learning performance was instead poor in case of slow stimulation speeds [41]. sequence learning. Given our results, one might speculate that strong synapses are important for memorizing varying temporal sequence speeds.
In our network, we included STD, which can be very commonly observed in neural systems and might serve various computational functions. In our simulations, STD was beneficial in reducing the likelihood of persistent activity which occurs otherwise for strong excitatory connection settings, that is, moderate to high settings of the g NMDA and g AMPA conductances. Since strong excitatory connections are on one hand important for the ability to train and recall sequential activity in the network, but persistent activity would on the other hand wipe out any learned weight structure, including STD effectively enlarged the parameter range of g NMDA and g AMPA where successful activity recall could be found. Further simulations indicated that a network without STD hardly performed activity recall with very slow speeds, no matter how one changed g GABA and g NMDA (not shown). Furthermore, we showed that the amount of spurious recurrent network activity can have negative effects on learning the required weight structure, so that incorporating STD might be important, because it tends to reduce ongoing activity unrelated to the sequence. Other studies analyzing long-term potentiation on sequential activity propagation usually neglect STD for simplicity [25,26,38]. Our results suggest that STD might nevertheless be important for modeling realistic sequence recall.
We noticed further that the type of the STDP learning rule affected our results. In particular, only the multiplicative NN STDP rule allowed the network to learn very slow recall speeds. In [50] it was suggested that, NN spike interaction could be achieved biologically by a reset of the membrane voltage caused by the backpropagating post-synaptic spike that effectively overrides the recent spiking history. Analogously, calcium saturation or glutamate receptor desensitization caused by the first postsynaptic spike might override the effect of all following spikes on the synaptic modifications [50].
Spike dependent synaptic modifications have been observed in many areas of the brain but the details of the implemented STDP rule might differ between brain areas [52]. Our results indicate that all-to-all STDP rules are less in agreement with the experiment as even very small stimulation speeds resulted in very fast recall speeds. Thus, the learning was 'too good'. Whether other STDP rules such as STDP rules that include partial contribution of other spike pairs or rules that involve spike triplets [63][64][65][66][67] could possible reproduce the experimental results remains open and subject for future research.
We used a strong stimulus drive during learning, resulting in a high population firing rate. This high stimulus strength was chosen to be comparable to the stimulus in the experimental study [41], where similarly large light spots were used that caused a high frequency response in the recorded visual neurons in the anesthetized state. In the network models, we additionally used a smaller stimulus sizes and found that temporal sequences could also be trained. However, smaller stimulus sizes resulted in larger recall speeds, probably because the temporal stimulation profile in neighboring single neurons was less overlapping. Note that the contribution of spike pairs in the overlapping stimulation period tended to suppress the learning of the asymmetric weight structure needed for successful sequence recall (as shown in the analysis of figure 5). Whether this behavior could also be seen in experiment remains open, because the experimental study did not vary the size of the light spot in a systematic way.
In the experimental study, the sequential activity recall was not always successful [41]. The success rate was about 40% in anaesthetized state (personal communication with Shengjin Xu), and even smaller in the awake state. In our simulations, the rate of successful recall appears similarly probabilistic in that a recall could halt half-way without reaching the endpoint. Interestingly, the recall can spontaneously re-occur following a partly recall (driven by the longer-lasting NMDA effects) and only then successfully transmit to the endpoint (compare to figure 3(c)). Because we measured the recall speeds by the time when the first spike reaches the endpoint, these incomplete recalls decrease the overall recall speed of the network. Increasing the g AMPA level does mainly decrease the probability of incomplete recalls, and thus the average network recall speed becomes faster.
In our network simulation the spontaneous activity was relatively low which corresponds to the situation in the anesthetized state (see [41]). In awake animals, on the other hand, spontaneous activity is known to increase (e.g. [41]). Since spontaneous activity will randomly activate neurons, it might affect the ability to learn and recall sequences in the network. While a careful investigation of levels of spontaneous activity on the recall behaviors under different STDP rules and stimulus speeds is beyond of the scope of this study, additional control experiments indicated that, as long as the stimulus is strong enough, successful learning and retrieval was still possible with moderate levels of spontaneous activity (2-3 Hz before training). We thus conclude that increasing the amount of spontaneous activity to moderate levels would not crucially affect our results.
In the recall experiments [41], the recorded length of visual cortex areas is 2 mm. Because transmission speeds of axons are at least not less than 0.3 m s -1 in the visual cortex [62], the maximal transmission delay should be not greater than 7 ms for the whole recorded length. In our spiking network, we used a chain-like structure of = N 200 e excitatory neurons to simulate the recorded cortex area on the path of the stimuli trajectory. We adopted transmission delays that linearly depended on the neuronal distances in the first part of the simulations (0.1 ms/neuron). Since the possible longest synaptic distance is 200 (distance from the first excitatory neurons at starting point to the last excitatory neurons at the end point), the maximal delay was about 20 ms in our simulations, thus larger than in the recall experiments if length are taken as comparable. Thus the observed fast speed in AMPA-dominant simulations is not a result of a too small delay in our networks. We found that the choice of the type of transmission delay did not qualitatively influence our results that there are two different recall behaviors, one AMPA-dominant fast and reliable sequential recall and one slow NMDAdominant rate-based recall. Similarly, when increasing the length of the network (for instance by doubling the neuron number), the qualitative results presented here are not expected to change, as long as the local connections structure (i.e. the spatial width of the recurrent connections) remains identical.

Conclusion
We here investigated in simulations how sequential spike activity could be learned and recalled in SNNs. The complex interplay between long and small time scales in the network, e.g. NMDA, AMPA, STD, STDP, and stimulation speed as well as the influence of recurrent spikes on the long-term weight modifications indicates that neuronal learning in realistic networks faces many challenges. In particular, our results show that weight modification and recall patterns could strongly depend on the state of the neural network.
This state-dependence of spiking activity in local neural networks might have important implications for the global dynamics of a physiological system that functional connects several local neural networks into a complex physiological network. Since the global state of the physiological system might change dynamically depending on the state and the functional connectivity to other subsystems in the physiological system, as recently shown for e.g. the relation between several organ systems and frequency bands in the brain during switching between sleep and awake [1], our results indicate that properties of the microscopic spiking dynamics in local neural networks might be intricately linked to global physiological state transitions. However, what the functional importance of this link is, will be subject for future studies.