## Abstract

Spike-timing-dependent plasticity (STDP) is an important synaptic dynamics that is capable of shaping the complex spatiotemporal activity of neural circuits. In this study, we examine the effects of STDP on the spatiotemporal patterns of a spatially extended, two-dimensional spiking neural circuit. We show that STDP can promote the formation of multiple, localized spiking wave patterns or multiple spike timing sequences in a broad parameter space of the neural circuit. Furthermore, we illustrate that the formation of these dynamic patterns is due to the interaction between the dynamics of ongoing patterns in the neural circuit and STDP. This interaction is analyzed by developing a simple model able to capture its essential dynamics, which give rise to symmetry breaking. This occurs in a fundamentally self-organizing manner, without fine-tuning of the system parameters. Moreover, we find that STDP provides a synaptic mechanism to learn the paths taken by spiking waves and modulate the dynamics of their interactions, enabling them to be regulated. This regulation mechanism has error-correcting properties. Our results therefore highlight the important roles played by STDP in facilitating the formation and regulation of spiking wave patterns that may have crucial functional roles in brain information processing.

## 1. Introduction

Spike-timing-dependent plasticity (STDP) is a temporally asymmetric form of Hebbian learning induced by tight temporal correlations between the spikes of pre- and postsynaptic neurons. In this type of plasticity, which involves both long-term potentiation (LTP) and long-term depression (LTD) of synapses, firing of a presynaptic neuron immediately before a postsynaptic neuron results in LTP of the synaptic transmission and a reversal of firing order results in LTD (Abbott & Nelson, 2000; Dan & Poo, 2006; Han, Caporale, & Dan, 2008). It has been shown that STDP can have significant effects on the collective dynamics of neural circuits; STDP tends to enhance population synchrony (Arieli, Shoham, Hildesheim, & Grinvald, 1995; Levy, Horn, Meilijson, & Ruppin, 2001; Nowotny, Zhigulin, Selverston, Abarbanel, & Rabinovich, 2003), and other interesting effects can arise depending on what model features STDP is combined with. For instance, in randomly coupled networks, STDP, together with heterosynaptic competition, facilitates the formation of multiple spike timing sequences (Fiete, Senn, Wang, & Hahnloser, 2010), and STDP with an update rule that is multiplicatively dependent on synaptic weights facilitates the formation of balanced states with irregular spiking activity (Morrison, Aertsen, & Diesmann, 2007). In randomly coupled networks with axonal conduction delays, it has been found that STDP mediates the development of strongly connected neuronal groups (Izhikevich, Gally, & Edelman, 2004) and enables neural circuits to self-organize to a state at the border between randomness and synchrony when there are neural population bursts in such networks (Lubenov & Siapas, 2008). It has been found that STDP, when applied to feedforward networks, can enhance the formation of synfire chains (Suri & Sejnowski, 2002; Jun & Jin, 2007). Although these studies have shown the importance of STDP in shaping the collective dynamics of neural circuits, they have primarily focused on randomly coupled or feedforward neural circuits that do not have complex spontaneous spatiotemporal dynamics. Accordingly, there is relatively little theoretical understanding of the effects of STDP on the pattern dynamics of spatially extended neural systems in which the local dynamics are coupled between neighboring neurons in two dimensions, a property that is conspicuously absent in systems such as feedforward or randomly coupled networks but widespread in real neural systems (Kandel, Schwartz, & Jessell, 2000).

Spatiotemporal activity emerging from spatially extended neural systems exhibits intriguingly organized patterns in both space and time, which are commonly observed in multiunit electrophysiological recording, EEG local field potential recording, MEG, and optical imaging, in both spontaneous activity (Arieli et al., 1995; Han et al., 2008) and evoked responses (Rubino, Robbins, & Hatsopoulos, 2006; Ferezou et al., 2007; Wu, Huan, & Zhang, 2008). In space, these patterns often take the form of localized patches of activity. Recordings over large populations of neurons have shown that several of these localized patterns can occur simultaneously across cortical regions (Kleinfeld & Delaney, 1996; Freeman & Barrie, 2000) and that these patterns often propagate in space. Despite many efforts, a detailed understanding of the neurobiological factors contributing to the formation of such dynamic wave patterns in populations of neurons, particularly in spiking neural circuits, is still developing (Pinto & Ermentrout, 2001; Folias & Bressloff, 2004; Coombes, 2005; Coombes & Owen, 2005; Gong & van Leeuwen, 2009; Bressloff & Kilpatrick, 2011; Lu, Sato, & Amari, 2011; Gong & Robinson, 2012). It has been suggested that such wave patterns might play important functional roles in transferring information (Rubino et al., 2006; Gong & van Leeuwen, 2009) and in carrying out distributed dynamical computation based on their collisions (Gong & van Leeuwen, 2009). A question that is at least as important as the formation of such dynamical patterns is the learning and regulation of their dynamics; how this can be achieved by neurophysiological mechanisms remains unclear.

The propagation of these dynamic spiking patterns through a substrate of neurons can naturally generate spike timing sequences. Such sequences have been observed in various parts of the brain, including cortex (Moran & Schwartz, 1999; Ikegaya et al., 2004; Tang et al., 2008), hippocampus (Nádasdy, Hirase, Czurkó, Csicsvari, & Buzsáki, 1999; Louie & Wilson, 2001), basal ganglia (Barnes, Kubota, Hu, Jin, & Graybiel, 2005), and the songbird high vocal center (Hahnloser, Kozhevnikov, & Fee, 2002), and they are supposed to be important for brain information processing (Abeles, 1991). STDP is an attractive mechanism for producing these sequences as it depends on the precise firing order of neurons and has been used to create spiking sequences in networks with synfire chains (Suri & Sejnowski, 2002; Jun & Jin, 2007) and in random networks (Fiete et al., 2010). However, these studies did not concern themselves with how multiple sequences can emerge from spatially extended neural circuits with complex ongoing dynamics.

Here we use a spatially extended, two-dimensional spiking neural circuit model representing the gross architecture of the cerebral cortex to study the effects of STDP on the dynamics of spatiotemporal patterns in neural systems. We show that STDP can promote the formation of localized, propagating wave patterns or multiple spike timing sequences in a broad parameter space of the neural circuit. We illustrate that the formation of such dynamical patterns is due to the interactions between ongoing pattern dynamics of the neural circuit and STDP, which gives rise to symmetry breaking of initial pattern shapes and the coupling topology, happening spontaneously over time. In addition to facilitating this self-organized symmetry breaking, STDP provides a synaptic mechanism to learn the paths taken by spiking wave patterns and modulate the dynamics of their interactions. This therefore enables these dynamic patterns to be regulated.

## 2. Model

*n*-by-

*n*two dimensional lattice of integrate-and-fire (IF) neurons. We denote the voltage of a neuron at integer coordinates (

*i*,

*j*) at time

*t*as

*V*(

_{ij}*t*), which has dynamics governed by where ms is the neural time constant, is the time step,

*I*is a constant external input, and

_{ex}*I*(

_{ij}*t*) is the input received from other neurons at time

*t*. Each neuron receives both excitatory and inhibitory input from other neurons. A spike is generated whenever the neuron reaches a threshold voltage

*V*=1. After emitting a spike, the neuron remains quiet for a refractory period

_{th}*t*. For simplicity, we choose ms. The coupling strength between any two neurons located at (

_{ref}*i*,

*j*) and , respectively, is denoted as , constructed from a Mexican hat function: where is the Euclidean distance between neurons on a lattice with periodic boundary conditions (i.e., a torus structure), although we have found that the symmetry breaking discussed later occurs in networks with open boundary conditions as well. are constants that determine the shape of the coupling function. As interactions between neurons in real neural systems have finite ranges, connections between neurons in the model are constrained to . For these constants, coupling is excitatory if and inhibitory if

*D*

_{1}>

*d*>

_{ij}*D*

_{0}; this therefore results in a center-surround coupling structure, as is commonly found in real neural systems (Schuett, Bonhoeffer, & Hübener, 2002). Note that we have found that the collective dynamics of the neural circuit is not sensitive to the values of any of the above parameters. To enable the total excitatory and inhibitory coupling strengths to be changed without varying their spatial extent, they are normalized separately using the following normalization constants: The coupling strength between neurons (

*i*,

*j*) and is then given by where

*W*and

_{E}*W*are the total excitatory and inhibitory coupling strengths, respectively. With this coupling function, the input the neuron (

_{I}*i*,

*j*) receives from other neurons to which it is coupled is where is the Heaviside step function: if and if

*x*<0.

*t*is the

^{k}_{ij}*k*th firing time of neuron (

*i*,

*j*) and is the

*l*th firing time of neuron , the sum is over all spike pairs before time

*t*, and is the STDP window function where

*t*and

_{pre}*t*are the firing times of the pre- and postsynaptic neurons respectively, and . As typically used in other studies (Song, Miller, & Abbott, 2000; Izhikevich et al., 2004), the parameters of STDP are chosen such that

_{post}*H*(

*t*, 0) is odd, that is, and ms. We have varied the values of

*A*

_{+}and

*A*

_{−}with values such as and have found that the main results discussed below do not change. We obtain the time-dependent coupling strength used in this study by adding the synaptic dynamics caused by STDP (see equation 2.7) to the time-independent coupling (see equation 2.5): STDP may cause some coupling strengths to reach values that are not within biological limits. To prevent this, STDP requires a hard limit of synaptic modification (Song et al., 2000; Izhikevich et al., 2004), without any further constraints. We limit changes to coupling strengths to be within 50% of their initial value, although the results are not sensitive to this value, which bounds the time-dependent coupling strength as follows: where is equal to the initial connection strengths as described in equation 2.5.

*I*(

_{ij}*t*) in equation 2.1 is then given by . We have found that none of our results are strongly affected by adding this or other types of noise to the model.

## 3. Network Dynamics Without STDP

Before studying the 2D spiking neural circuit with STDP, we briefly review its collective dynamics without STDP as described by equations 2.1 and 2.6. Starting from random initial conditions, the system self-organizes into localized spiking patterns with complex collective dynamics. For the network with a size of , we observe two particular types of spatiotemporal patterns. First, multiple; localized spiking patterns with a near-isotropic shape with slight variability (type 1) wander around a fixed spot in a seemingly random way without direction preferences; they do not move together, as shown in Figure 1A. Second, multiple, localized wave patterns with crescent-shaped fronts (type 2) have long-range propagation. As patterns are traveling in different directions, they collide when they approach each other, as shown in Figure 1B. Figure 1C maps out where these types of patterns occur in the parameter space of *W _{E}* and

*W*for a network of size 100 × 100 neurons. For smaller networks with a size of 80 × 80 neurons, a third type of pattern, multiple localized wave patterns with crescent-shaped fronts propagating in the same direction, can exist.

_{I}*X*(

_{i}*t*),

*Y*(

_{i}*t*)) of the

*i*th pattern at time

*t*: where

*x*and

^{i}_{j}*y*are the

^{i}_{j}*x*and

*y*positions of the

*j*th neuron of the

*i*th pattern and

*N*is its total number of firing neurons. We use mean squared displacement (MSD) to characterize the motion of these patterns: where is the time increment and represents averaging over all time values

_{i}*t*and all patterns. For type 2 patterns, is linear on a log-log scale, as shown in Figure 2A. This can be verified using linear regression, indicating that is a power function of time increment (Gong & van Leeuwen, 2009): For type 1 patterns, appears to have two different scaling behaviors, as shown in Figure 2B. When , we find that with , indicating normal diffusion. When , we find that with , indicating anomalous subdiffusion. Generally the exponent of subdiffusion () increases as

*W*decreases; it increases from 0.1 to 0.7 as

_{E}*W*is decreased from 2 to 1.3. Such subdiffusive dynamics indicate the existence of coherence in pattern dynamics. To confirm that these patterns are subdiffusive over long timescales, we use an additional method for characterizing subdiffusive dynamics: the power spectrum of the pattern's motion where is the Fourier transform of the time series of the

_{E}*x*-component of the center-of-mass position of the

*i*th pattern. We obtain the averaged power spectrum

*P*(

*f*), by averaging over all patterns. For patterns with subdiffusive dynamics, the power spectrum should have a part that approximates a power law— where is the same as in equation 3.3 (Golding & Cox, 2006). As shown in Figure 2C, this is indeed the case, showing that the pattern dynamics over long time scales is subdiffusive.

In the parameter regime of type 1 patterns in our model, if only a single pattern is present, we find that it would move in a normal diffusive way; similar dynamics for a single pattern have been found in other deterministic spiking neural circuits (Chow & Coombes, 2006) as well. Multiple, localized patterns with normal diffusive dynamics have been found in a stochastic spiking neural circuit (Burak & Fiete, 2009). In our deterministic integrate-and-fire network, multiple, localized patterns, however, have normal diffusion when but subdiffusion when . The normal diffusion over short timescales arises mainly from the intrinsic fluctuations of spikes within individual patterns. To gain a further understanding of the subdiffusive motion, we calculate the total input that individual patterns get from their surrounding neighbors, *I _{total}*(

*t*). As shown in Figure 2D, there is long-range autocorrelation in the input. Such correlated inputs or interactions from other patterns may result in the subdiffusive motion, as analogous to that occurring in other physical systems (Viñales & Despósito, 2007).

## 4. Effect of STDP on the 2D Spiking Neural Circuit Model

After reviewing the collective dynamics of the 2D neural circuit without STDP, we now show how they are affected by STDP. We first focus on the parameter regime where the original system without STDP has multiple, localized patterns with subdiffusive dynamics. After applying STDP to the model, equation 2.9, the pattern dynamics changes over time, and they are characterized by . As shown in Figure 3A, initially the dynamics are identical to that generated in the model without STDP; that is, localized patterns with wandering behavior are present. After 2000 ms, the system dynamics rapidly shifts into a new regime; the transition process happens over a short time period of around 200 ms. In the new regime, all patterns have long-range propagation, with dynamics characterized as superdiffusive over long timescales—, as shown in Figure 3B. The shape of these patterns, however, is still approximately rotationally symmetric but with slight variability. This can be characterized by using an order parameter (see section 4.1) that has a small but nonzero value. The next change to the dynamics of patterns occurs gradually over the next 2000 ms. During this period, some patterns will begin to propagate with crescent-shaped fronts, while others maintain their rough circular shape.

In the final stage, all dynamic patterns have changed from being approximately circularly symmetric localized patterns to crescent-shaped propagating wave fronts. Their propagation speed approaches a constant value of approximately 1 grid unit per time step, although the precise speed is dependent on *W _{E}* and

*W*. This speed is not influenced by the time discretization; changing the size of does not significantly affect the propagation speed of these wave fronts. In this stage, , indicating that these patterns now have a more regular, ballistic propagation. A collective preferred propagation direction has emerged: all patterns propagate in approximately the same direction. This process persists until synaptic strengths have reached the limit given in equation 2.10. Figure 3G shows that during this process, the collective dynamics of these patterns has changed from being subdiffusive to superdiffusive and then to approximately regular motion—initially for , before increasing until

_{I}*t*=10

^{4}ms, at which point . As wave patterns travel with constant speed, is an upper bound indicating ballistic propagation. This confirms that all patterns are propagating in the same direction for large

*t*. Note that we have observed that these changes to the synaptic dynamics are not sensitive to changes in the system parameters: for all (

*W*,

_{E}*W*) values that produce type 1 patterns parameterized in Figure 1C, a similar three-stage transition occurs, resulting in traveling crescent-shaped wave fronts for sufficiently large

_{I}*t*. To study the robustness of these dynamics against random perturbations, we also study the stochastic model, with synaptic coupling strengths described by equation 2.11. We find that a similar three-stage transition in the collective dynamics of the neural circuits happens, as shown in Figures 2D to 2F. This therefore indicates that this transition is not sensitive to perturbations added to the model.

### 4.1. Self-Organized Symmetry Breaking.

*k*ranges over the firing neurons in an individual localized pattern, (e.g., for a typical pattern shown in Figure 3,

*N*=7),

_{j}*j*is the index of the wave pattern, and is the azimuthal angle of the

*k*th firing neuron of the wave pattern with respect to a fixed reference axis. For patterns with a high degree of rotational symmetry, the magnitude of the order parameter should be small and very close to 0; for nonrotationally symmetric patterns, it should be significantly higher. Similar-order parameters are often used to study the degree of rotational symmetry in a wide variety of physical systems (Sethna, 2006). To construct an order parameter for the coexisting, multiple patterns, is averaged over all concurrent patterns, and then magnitude is taken:

Figure 4 shows how changes as time proceeds. For *t*<*t*_{1}=2000 ms, and changes little, this small but nonzero value indicates some heterogeneity or variability in the pattern at first, mainly due to the random initial conditions of the system; we refer to this as a roughly symmetrical pattern. For *t*>*t*_{1}, starts to increase, indicating that the spiking neurons within individual localized patterns tend to condense to some preferred angle and break the initial rough rotational symmetry. As we will show, this increase is due to the interaction of pattern dynamics with STDP. Finally, reaches a maximum value of 0.18. We note that without STDP, no such increase occurs, and for all time. To show the robustness of this behavior, we calculate the order parameter in the network with stochastic synapses. As shown in Figure 4, the magnitude of the order parameter undergoes similar dynamics as in the deterministic model. The robustness of this symmetry breaking is also shown by the fact that it persists across a wide range of parameter values; specifically, it occurs for all parameter values producing type 1 patterns in Figure 1C, in both the deterministic and stochastic models. It is interesting to note that as quantified by the order parameter shown in Figure 4, the symmetry breaking happens spontaneously as time proceeds in a self-organizing manner. We therefore refer to this as self-organized symmetry breaking, in contrast to previous studies where fine tuning of system parameters is required to break the rotational symmetry of two-dimensional localized patterns (Bressloff & Kilpatrick, 2011; Gong & Robinson, 2012).

We next examine how STDP drives self-organized symmetry breaking. Without loss of generality, we choose a single neuron located at (50, 50) and examine the changes of the time-dependent coupling strength, . We note that at *t*=0, ; coupling strengths are perfectly symmetrical (see Figure 5A). As shown in Figure 5, for *t*<*t*_{1}=2000 ms, coupling strength changes appear to be homogeneous, with no clear direction preference. For *t*_{1}<*t*<*t*_{2}=4000 ms, however, coupling strength changes begin to show a slight direction preference. For *t*>*t*_{2}, a clear directionality in coupling strengths exists: the connections in the northwest direction have been increased and those to the southeast have been decreased. This is consistent with the time frame of symmetry breaking characterized by the order parameter as shown in Figure 4. This therefore indicates that STDP has led to asymmetrical coupling strengths and the resultant self-organized symmetry breaking of the pattern shapes.

### 4.2. Interaction Between STDP and Ongoing Pattern Dynamics.

We now propose a mechanism by which STDP interacts with the pattern dynamics to cause self-organized symmetry breaking. As we have noted, each localized pattern wanders around; a typical trajectory is shown in Figure 6. Without loss of generality, we label the starting center-of-mass position of a localized pattern as *A*. As the pattern wanders, it leaves and then returns to this location many times, as shown in Figure 6. The center-of-mass location of the pattern immediately after it leaves *A* for the first time is then denoted as *B*_{1}. The pattern wanders around before eventually returning to *A* at some time *t*. We denote *C*_{1} to be the center-of-mass location of the pattern at time *t*−1, that is, the location of the pattern immediately prior to returning to *A* for the first time. The pattern then leaves *A* for the second time, wandering to location *B*_{2} at time *t*+*d* for some *d*>0. As shown in Figure 6, this process repeats many times. We denote the sequence *B _{i}* to be the location of the pattern immediately after it leaves

*A*for the

*i*th time. Similarly, we denote the sequence

*C*to be the location immediately prior to returning to

_{i}*A*for the

*i*th time.

*A*. Let the neurons firing when the wave pattern is at

*A*,

*B*and

_{i}*C*be , , and , respectively. We first consider the case when the pattern leaves position

_{i}*A*for the first time—when it moves from

*A*to

*B*

_{1}. By definition, the firing times of the neurons in will be before those in . As a result, neurons in are the presynaptic neurons, and those in are postsynaptic. According to the STDP rule, equation 2.8, the coupling strengths from neurons in to neurons in will be increased. Because connection strengths along the path from to have increased, the pattern is now more likely to follow this path to

*B*

_{1}when it next leaves

*A*. Similarly, as neurons in fire before those in and our STDP window function, equation 2.8, is odd, connection strengths from to will be decreased; this process creates a feedback mechanism, which we detail below. As a result, the probability that

*B*

_{2}=

*C*

_{1}has decreased. Based on the understanding of changes in the coupling strength and the resulting changes in the probability of propagation in a particular direction, we develop a simple stochastic model capturing the key dynamics. We assume that initially, all propagation directions for randomly wandering patterns are equally likely. We first consider the probability distribution for

*B*

_{2}: the next location after the pattern has moved away from

*A*and then returned for the first time. As a result of changes to coupling strengths discussed above, the probability for moving along the path from

*A*to

*B*

_{2}is given by where

*d*

_{0}=0.5 is a small distance that judges how far apart two location are considered to be identical (i.e., when

*d*(

*B*

_{2},

*b*

_{1})<

*d*

_{0}, we treat

*B*

_{2}is identical to

*B*

_{1}),

*d*(

*B*

_{2},

*B*

_{1}) is the Euclidean distance between

*B*

_{2}and

*B*

_{1},

*p*

_{0}is an initial uniform probability for all directions, is the Heaviside step function, where

*M*is the total number of possible locations for

*B*

_{2}based on a maximum velocity

*v*of 1 grid unit per ms, and is a small constant that reflects the change in probability caused by moving along the path from

*A*to

*B*

_{1}or

*C*

_{1}to

*A*. Assuming that changes in probability are proportional to coupling strength changes, equation 4.3 can be applied twice to obtain We generalize equation 4.4 to obtain the probability of a particular path being chosen as

*B*

_{n+1},

*P*(

*B*

_{n+1}), by applying equation 4.3

*n*times:

We now investigate the evolution of *P*(*B _{n}*) as

*n*increases. As shown in Figure 7A, for

*n*<10

^{4}, the probabilities are roughly uniformly distributed near their initial values. As

*n*increases, some directions become dominant until eventually only a single direction has nonzero probability (see Figure 7B):

*P*(

*B*) approaches a delta function as

_{n}*n*increases, indicating that localized patterns now propagate in only a single direction. This is repeated for different sets of random initial conditions, and we find that eventually, a single propagation will become dominant, although which direction is chosen is dependent on initial conditions, as shown in Figures 7A and 7B. This is consistent with what we observe in the full model. To further show convergence to a delta function over many different initial conditions in the simple model, we calculate the variability in

*P*(

*x*) for a fixed

*x*as time progresses, which we then average over many trials, each with a different initial condition. As show in Figure 7C, this variability is initially nonzero, before converging to 0, indicating a single dominant propagation direction. The dynamics of this stochastic model suggests that the interactions between STDP and the underlying dynamics of the localized patterns generate a feedback mechanism; that is, stronger (weaker) connections result in more (less) propagation probability along that path, which is further strengthened (weakened). This feedback mechanism eventually gives rise to the symmetry breaking happening in a self-organizing way in the full network model, self-organized symmetry breaking.

In the model, as lateral inhibitory coupling is present and the range of coupling is larger than the size of spiking patterns, the interactions between patterns are repulsive. This repulsive interaction causes type 1 patterns to form a rough hexagonal grid. When STDP is included, if one particular pattern begins to move away from its location in the hexagonal grid of patterns, it will influence other patterns. The patterns that it moves toward are repelled, and those that it moves away from are effectively pushed toward the newly created gap in the lattice, as the force repelling them from the original pattern is weakened. As a result, patterns will tend to move in the same direction as the original pattern. This movement is then amplified by the feedback mechanism created by STDP, which amplifies coupling strengths in the direction of motion. Because these patterns always interact through this repulsive interaction, we see all patterns moving in the same direction after symmetry breaking, as shown in Figure 3C.

## 5. Learning and Regulation of Dynamical Wave Patterns

We have shown that STDP facilitates the formation of localized spiking wave patterns without fine tuning of parameters. As wave patterns might play important functional roles in carrying out distributed dynamical computation in neural systems (Gong & van Leeuwen, 2009), it is interesting to study how the dynamics of these patterns such as their propagating paths can be learned. On the other hand, the mechanism of learning of spiking sequences is still unclear (Wörgötter & Porr, 2005), so learning the propagating paths of these patterns might provide some insight into understanding the learning of spiking sequences as well.

First, we show that the propagation paths of these dynamic patterns can be learned using STDP. We first set *I _{ex}*=0.045, so that no spontaneous firing of neurons occurs in the system. We then apply a brief, localized stimulus that is able to evoke a localized crescent-shaped pattern propagating away from the stimulus location, as shown in Figure 8A; the evoked spiking wave pattern behaves similarly to those shown in Figure 1B; it is capable of long-range propagation. When STDP is included in the circuit, this procedure is repeated to enable learning to happen. After a single trial, the time-dependent coupling strengths (see equation 2.9) have been modified by the passage of the pattern. As shown in Figure 8B, coupling strengths parallel to the propagation direction of the wave pattern are increased, and those in the opposite direction are decreased. After around 15 repetitions, the network learns the propagation path of the pattern. As a result, after applying external stimuli around the initial location, a wave pattern will return to the learned path and then propagate along it. If the location of the start of the wave is shifted by up to 6 grid units, we find that despite this relatively large perturbation, the wave pattern rapidly (within 10 ms) returns to the original path, as shown in Figure 8A, and then propagates along it. Such perturbations in stimulus locations can be regarded as errors in the initial stimulus. This dynamic behavior of returning to the learned path therefore indicates that the dynamics of the patterns now has an error-correcting property. This error-correcting property is analogous to that shown in conventional attractor neural networks (Hopfield, 1982), in which small perturbations of state away from a learned attractor are corrected as the state quickly returns to the attractor. When waves return to learned paths, they have also been regulated, as propagation will now occur preferentially along the learned paths.

Next, we illustrate that STDP can be used to modulate the collision dynamics of two wave patterns. We apply two stimuli simultaneously, which evoke two wave patterns that collide with each other, as illustrated in Figure 8C. We repeat this learning process 15 times and notice changes to the behavior of the wave patterns after the collision. Initially waves travel through the collision with a small change in propagation direction (). However, after repeating several trials, as shown in Figure 8C, this perturbation angle has increased so that the waves travel perpendicular to their original propagation after the collision. After learning, when the positions of the two wave patterns are shifted, they quickly converge to the learned paths and collide at the original collision location; that is, the whole collision process has the error-correcting property as well.

Even when there are multiple spontaneous dynamic patterns of type 2, as shown in Figure 1B, it is possible to use STDP to regulate their propagating paths. We illustrate this with several artificial examples. First, we try to regulate all waves to propagate in a specific direction . To illustrate this, we add a constraint to our implementation of STDP. Let be the instantaneous propagation direction of wave pattern *i* at time *t*. We now apply STDP if the wave pattern is traveling close to the desired direction: . As a result, the wave paths that propagate in the desired direction are learned by the network, and any other paths are not learned. Figure 9A shows that after learning, all waves propagate in the same direction. This occurs if the initial parameters produce either type 1 or type 2 patterns. STDP can also be used to regulate collisions so that they occur at a particular location when there are multiple dynamic patterns. To demonstrate this, we apply STDP only if a wave is traveling along one of two perpendicular paths meeting in the center of the network (or any other location). Applying STDP in this way encourages waves to travel along these paths until they collide where the paths meet, as shown in Figure 9B. Waves that are not traveling along the paths may propagate freely. As a result, waves traveling along the desired path are frequently involved in collisions.

## 6. Discussion and Conclusion

By constructing a two-dimensional spiking neural circuit with Mexican hat connectivity, we have studied theoretically the effect of STDP on the dynamics of complex spatiotemporal patterns in neural systems. In our network model, the dynamics of individual neurons are represented by commonly used integrate-and-fire spiking neurons (Burkitt, 2006), and the network topology is based on the center-surrounded intracortical interactions that have been widely observed in real neural systems (Buzás, Eysel, Adorján, & Kisvárday, 2001). We find that the interactions between the dynamics of STDP and that of the ongoing patterns in the system can result in the formation of crescent-shaped, localized wave patterns by the self-organized symmetry breaking of synaptic weights. This occurs naturally as time proceeds in a broad parameter space for both deterministic and stochastic networks; in other words, it occurs in a fundamentally self-organizing manner. Aside from facilitating the formation of the dynamic patterns, we find that STDP also provides a necessary neurophysiological mechanism to learn the propagating paths of these patterns and modulate their colliding dynamics. These mechanisms also have error-correcting properties, as has been shown in this study for a few patterns and interactions between a pair of patterns.

The effect of STDP on the dynamics of neural systems has been mainly studied in feedforward or randomly coupled neural circuits that do not exhibit significant, complex spatiotemporal dynamics (Suri & Sejnowski, 2002; Kitano, Câteau, & Fukai, 2002; Izhikevich et al., 2004; Morrison et al., 2007; Jun & Jin, 2007; Lubenov & Siapas, 2008; Fiete et al., 2010). This letter therefore extends the study regarding the effect of STDP on neural dynamics by considering spatially extended, two-dimensional neural systems, which certainly come closer to describing real neural systems than previously studied randomly coupled or feedforward systems. As illustrated in this study, the interactions between the dynamics of STDP and that of ongoing patterns result in symmetry breaking that happens in a self-organizing manner. These interactions, as outlined above, act as a positive feedback process shaping the pattern dynamics: larger (smaller) weights along some direction increase (decrease) the chance of a localized spiking pattern moving in this particular direction and are further strengthened (weakened) by STDP. The repulsive interactions of many patterns caused by lateral inhibitory coupling then lead to their approximately parallel paths.

Since the propagation of these localized spiking patterns through a substrate of neurons can naturally generate spiking sequences, this study also proposes a novel mechanism to generate multiple spiking sequences, clearly demonstrating that the formation of such spiking sequences is feasible in spatially extended neural systems with STDP (see Itskov, Vinnik, & Diamond, 2011, in which a two-dimensional network with spike threshold adaptation instead of STDP is used to form spiking sequences). This is in contrast to previous studies of multiple spiking sequences in randomly coupled neural circuits, where STDP has to be combined with other synaptic dynamics such as synaptic competition (Fiete et al., 2010), or fine tuning of system parameters is necessary to prevent runaway excitation when forming multiple spiking sequences (Mehring, Hehl, Kubo, Diesmann, & Aertsen, 2003; Kumar, Rotter, & Aertsen, 2010). Our results therefore indicate the importance of considering spatial extension aspects of neural circuits in studying the dynamical effects of STDP.

The formation of wave patterns in two-dimensional neural circuits has been mainly studied in neural field models with consideration of mean firing rates only (Bressloff & Kilpatrick, 2011; Lu et al., 2011). These studies showed that firing rate adaptation or synaptic depression is essential in breaking the rotational symmetry of a two-dimensional pattern to form traveling wave patterns. For spiking neural circuits, it has been shown recently, that the refractory period of time plays an important role in facilitating the formation of spiking wave patterns (Gong & Robinson, 2012; Gong, Loi, Robinson, & Yang, 2013). In all these studies, some sort of tuning to system parameters was used to get symmetry breaking in pattern forms; for instance, in Gong and Robinson (2012) and Gong et al. (2013), the period of the refractory state was tuned, and in Bressloff and Kilpatrick (2011) and Lu et al. (2011), the coupling strength was tuned. In contrast to these studies, the symmetry breaking in pattern shapes revealed in our work happens in a self-organizing way in a broad parameter space. Our study, along with the previous studies, also indicates that a wide range of neurobiological factors, such as firing rate adaptation, synaptic depression, refractoriness and STDP, is capable of promoting the formation of localized wave patterns, suggesting that the resultant wave patterns are widespread in the brain and may play important functional roles in brain information processing.

A recent numerical study of a 2D neural field model showed that external inputs are able to control the propagating paths of traveling bumps (Lu et al., 2011). In our study, however, we have shown that an intrinsic neurophysiological factor, STDP, is able to learn the paths taken by traveling wave patterns. As far as a path has been learned, after applying external stimuli around its initial location, a wave pattern will return to the learned path and propagate along it. Importantly, this can happen in an error-correcting way. Such spiking wave patterns have been proposed to be basic building blocks for neural systems to carry out distributed dynamic computation (Gong & van Leeuwen, 2009); information is encoded in spiking wave patterns, information is communicated by propagation of these patterns, and information is processed when they collide with each other. In the future, it would be interesting to study whether some perceptual or cognitive functions can be understood by using interacting wave patterns regulated by STDP with error-correcting properties.