Refinement and Pattern Formation in Neural Circuits by the Interaction of Traveling Waves with Spike-Timing Dependent Plasticity

Traveling waves in the developing brain are a prominent source of highly correlated spiking activity that may instruct the refinement of neural circuits. A candidate mechanism for mediating such refinement is spike-timing dependent plasticity (STDP), which translates correlated activity patterns into changes in synaptic strength. To assess the potential of these phenomena to build useful structure in developing neural circuits, we examined the interaction of wave activity with STDP rules in simple, biologically plausible models of spiking neurons. We derive an expression for the synaptic strength dynamics showing that, by mapping the time dependence of STDP into spatial interactions, traveling waves can build periodic synaptic connectivity patterns into feedforward circuits with a broad class of experimentally observed STDP rules. The spatial scale of the connectivity patterns increases with wave speed and STDP time constants. We verify these results with simulations and demonstrate their robustness to likely sources of noise. We show how this pattern formation ability, which is analogous to solutions of reaction-diffusion systems that have been widely applied to biological pattern formation, can be harnessed to instruct the refinement of postsynaptic receptive fields. Our results hold for rich, complex wave patterns in two dimensions and over several orders of magnitude in wave speeds and STDP time constants, and they provide predictions that can be tested under existing experimental paradigms. Our model generalizes across brain areas and STDP rules, allowing broad application to the ubiquitous occurrence of traveling waves and to wave-like activity patterns induced by moving stimuli.


Introduction
After an initial stage of activity-independent construction [1], the developing nervous system undergoes a period of refinement that is strongly influenced by spontaneous and evoked patterns of neural activity [2,3]. Traveling wavefronts are a striking feature of these activity patterns [4][5][6][7][8][9]. Within short temporal windows, wavefronts induce strong interneuronal correlations that can act through Hebbian mechanisms of synaptic plasticity to build structure into the connectivity of neural circuits [10,11]. This has prompted the hypothesis that correlated activity plays an instructive role for circuit refinement in the developing brain (reviewed in [2] and [3]). One Hebbian mechanism that is well suited to this role and is widely reported in the brain is spike-timing dependent plasticity (STDP), for which synaptic connections are strengthened or weakened depending on the relative timing of pre-and postsynaptic spikes that arrive at the synapse, typically within tens of milliseconds of each other [12][13][14]. In this article, we undertake a mathematical analysis of the interaction between traveling wave activity patterns and STDP, and explore the types of connectivity patterns that emerge as a result of this interaction.
Past studies have demonstrated that STDP could translate correlated input patterns into structured neural circuits [15][16][17][18][19][20][21], and could mediate the construction of realistic receptive fields (RFs) with properties that resemble those found in the visual cortex [22,23]. These models primarily focussed on spatial and temporal correlations as separable features when considering their interaction with STDP. However, the spatiotemporal correlations induced by traveling waves are space-time inseparable, providing additional information that may be utilized during circuit building. Space-time inseparable activity patterns map the temporal profile of the STDP rule into a spatial profile of synaptic strength changes [24], which can be used to build circuits that mimic neuronal sensitivity to visual motion during repeated exposure to moving visual stimuli [25][26][27][28]. But despite the demonstrated applicability of STDP to specific cases of neural circuit development, a more general, formal analysis of the interaction of STDP and wave-like activity patterns is still lacking.
Here, we derive a mathematical expression that accounts for the interaction of a variety of traveling wave activity patterns and STDP rules, and we examine the analytical predictions in a simple yet biologically plausible model of spiking neurons. We show that, for a broad class of experimentally observed STDP rules, such interactions build highly structured and periodic connectivity patterns into feedforward circuits, analogous to a Turing instability in reaction diffusion systems [29], which has been applied to diverse cases of biological pattern formation [30,31]. We then demonstrate the robustness of this pattern formation process and how it can be utilized to construct and refine the size and shape of RFs. Our results offer theoretical insights that may advance the understanding of the role played by traveling wave phenomena in different areas of the brain. We highlight particular insights into visual system development and outline a number of predictions that may be tested experimentally.

Results
Our results are organized as follows. First, we describe how traveling waves can interact with STDP to influence synaptic strengths, derive a mathematical expression for this effect and verify analytical solutions to the derived equations using a simulation of spiking neurons. Second, we use further simulations to explore the robustness of our analytical results. Third, we demonstrate the properties of waves and STDP rules that allow for different types of refinement in downstream receptive fields (RFs).

Illustration of network dynamics
To understand how correlated activity caused by traveling waves could influence synaptic strengths via a STDP mechanism, we consider a reduced model (Fig 1A) consisting of a onedimensional (1D) layer of presynaptic input cells, all of which connect via excitatory synapses, w i,j , onto a single, postsynaptic output cell. In later sections, we extend the input layer to two dimensions. When a wave traverses the input population (Fig 1A), each input neuron is recruited by the wavefront (red colored unit with rightward arrow) and discharges a burst of spikes, which drives spiking in the output neuron. Temporal differences between the spike times of a given input neuron and the output neuron, Δt = t in −t out , determine how the respective synapse is modified in strength according to a STDP rule, K(Δt) (see Methods, Eqs 21 and 22). The set of spikes for all input neurons leads to a simple diagonal band structure in space and time (Fig 1B).
To illustrate how this input spike pattern leads to spatially structured changes in synaptic strength, we consider a STDP rule that is asymmetric in Δt (Fig 1C, top). In Fig 1D, we show a snapshot of the traveling wave and its influence on the surrounding synaptic strengths. Here, an input neuron (colored dark red) is recruited by the wavefront and fires a burst of spikes, which in turn elicits excitatory postsynaptic potentials (EPSPs) that drive the output neuron to spike. For waves traveling left to right, spikes generated by input neurons to the left of the wavefront always precede output spikes. Consequently, their respective synapses are strengthened, because the STDP rule specifies strengthening for negative Δt. Likewise, synapses connecting input neurons to the right of the wavefront will be weakened. In fact, due to the wave motion, the dependence of the STDP rule on relative spike times, Δt, is mapped onto the spatial axis of the input layer, Δx, and thus the relative spike locations, as shown in Fig 1D and 1G. Consequently, as depicted in Fig 1E, all input neurons (colored light red) that surround a synapse will influence the net change in strength at that synapse. One might predict this change to be proportional to the integral of the STDP rule, K 0 ¼ R 1 À1 dDt KðDtÞ. However, the extent to which an input influences surrounding synapses depends on its influence over the output firing rate, and thus its own synaptic strength. For example, with an asymmetric STDP rule and waves traveling left to right, a synapse will strengthen relative to K 0 if the surrounding synapses are stronger to the right than to the left, as in Fig 1F. On the other hand, a synapse will be relatively weakened if the surrounding synapses are stronger to the left than to the right. By the same argument, a synapse will be relatively weakened by a temporally symmetric STDP rule if the surrounding synapses are stronger on both sides (Fig 1G), and relatively strengthened if the surrounding synapses are weaker. In this way, strong synapses increasingly dictate changes to local synaptic strengths as more waves pass. Eventually, the connectivity pattern begins to form islands of strong synapses that are flanked by regions of weak synapses.
In the following analytical derivation and simulations, we show how this process of waveinduced STDP results in a distinct type of pattern formation in the network for a broad class of STDP rules.

Deriving synaptic strength dynamics
Here, we derive a description for the dynamics of synaptic strengths that are driven by pairs of input and output spikes acting through a STDP rule and resulting from traveling wave activity patterns traversing the input layer. Within our framework, input and output spike trains are The model consists of a layer of input neurons, along which wavefronts of spiking activity propagate. Input neurons undergo a burst of spikes, generated stochastically, when the wavefront passes (red unit). Input spikes elicit excitatory postsynaptic potentials (EPSPs) in the output neuron (blue unit) that are scaled by the synaptic strength, w ij . The output neuron generates spikes stochastically with a firing rate linearly proportional to its summed EPSP. B. Example spike raster from the simulation, showing one wave traveling past input neurons 1-10. C. Synaptic strengths are modified by either an asymmetric (top) or symmetric (bottom) STDP rule. D. Schematic for the influence of the wavefront on modifications at surrounding synapses. An input neuron (red) is recruited by a wave, which travels from left to right, and increases the firing rate of the output neuron (blue). When an asymmetric STDP rule is at play, output spikes at the current time point will cause synapses behind the wavefront to increase in strength ('+' symbol) because they were active at an earlier time. Similarly, synapses in front of the wave will decrease ('−' symbol), because they will become active at a later time. Thus, traveling waves map the STDP rule onto space. E. Schematic for the influence of surrounding inputs (colored light red), on synaptic modifications at the wavefront (dark red). F. Same as E, except that inputs to the right of the wavefront induce even greater synaptic strengthening as their respective synapses are stronger. G. Same as E, except for a symmetric STDP rule. Here, the greater strength of synapses either side of the wavefront induce more synaptic weakening. generated by a stochastic process with a time-dependent firing rate. The stochastic arrival of spikes results in stochastic changes to the synaptic strengths, thus posing a challenge when seeking a description for the spatial structure of synaptic strengths that evolves slowly over long periods of time, during which many traveling waves occur. It is therefore useful to separate the slow dynamics from the fast, stochastic dynamics under the assumption that, during a limited period of time, ΔT, individual changes in synaptic strengths are negligible, but accumulate slowly over multiple periods of ΔT as a result of the time-averaged input and output activity. By approximating synaptic strengths as being constant during the period ΔT, Kempter et al. [16] showed that changes in synaptic strength over this period could be described by the inner product of the STDP rule, K(Δt), and the cross-correlation function, C ij (Δt,t), between the spike trains of input neuron i and output neuron j: where η is a small, positive constant that sets the required slow rate of change in synaptic strengths. The cross-correlation is given by: where S i(j) (t) are ensemble averages, for example over multiple waves, of the input (output) spike trains and can thus be identified with the input (output) firing rates. It is important that the firing rates be sufficiently high for C ij (Δt,t) to accurately portray wave-induced correlations.
In addition, η must be sufficiently small for wave-induced correlations to be recovered over several waves. Moreover, without small η, calculating C ij (Δt,t) becomes difficult, as S j (t) would depend on stochastically changing synaptic strengths, w ij (t). As such, Eq 1 implements the approximation by averaging over the small stochastic fluctuations, thus providing only the mean drift in w ij (t). Given that we will deal with discrete waves that pass one-by-one across the input layer, it is convenient to relate the time scale, ΔT, to the passage time of just a single wave. Further assumptions are now required to uphold the validity of Eq 1. First, in order to ensure that multiple waves do not mutually influence changes in synaptic strength, ΔT must include an amount of time, K, both before and after the wave, where K is the temporal width of the STDP rule which contains most of its power. More formally, we require that R K ÀK dDt j KðDtÞ j ) [16]. Because we are effectively considering Δw ij (t) for a wave in isolation, we can extend the integral limits in Eqs 1 and 2 to ±1. Second, with w ij (t) effectively constant during a wave, and because Δw ij (t) is small, we will analyze changes in w ij on a slower time scale, T, which is discretized in ΔT increments. Thus, we approximate w ij (t) with w ij (T) and C ij (Δt,t) with C ij (Δt,T), and therefore approximate the left-hand side of Eq 1 with @w ij (T)/@T = @ T w ij (T).
By considering a simple 1D chain of input neurons with a single output neuron, we replace all subscripts i with the argument, x. That is, we replace w ij (T) with w(x,T), S i (t) with S in (x,t), and C ij (Δt,T) with C(x,Δt,T). Eq 1 then becomes: where Cðx; Dt; TÞ ¼ R 1 À1 dt S in ðx; t þ DtÞS out ðtÞ, with S in and S out the input and output firing rates, respectively. As a final matter of notation, we will hereafter refrain from explicitly writing the dependence of w and C on T, for brevity.
We make two further assumptions to simplify the analytical derivation, then relax these for a more general case: i) the input firing rate at the wavefront can be described as a short, traveling pulse using a Dirac delta function: S in (x,t) = δ(x−vt), where δ(y) = 1 if y = 0 and is zero otherwise, and v is the wave speed; ii) the output neuron's response to its input is instantaneous: S out ðtÞ ¼ R 1 À1 dx S in ðx; tÞwðxÞ ¼ wðvtÞ. Using these forms for S (in)out , C(x,Δt) = w(x−vΔt) and we can write Eq 3 as where K v (x) has been introduced as a rescaled copy of K, K v (x) = v −1 K(x/v), and Ã denotes convolution. There are two key features to Eq 4. Firstly, as illustrated above, the STDP rule can be reinterpreted as a spatial kernel as a result of the wavefront's constant velocity, which maintains a strict relationship between space and time. Secondly, the wave's effect on the synaptic dynamics is described by a convolution of the STDP rule with the synaptic strengths. By deriving a solution for w(x) in Eq 4, we will demonstrate how convolution plays an important role in the type of connectivity patterns that w(x) acquires, but first we relax the two simplifying assumptions used to reach Eq 4. Incorporating finite input bursts and the dependence of output firing rates on EPSPs adds a simple modification to Eq 4, which becomes where α(t) describes the time dependent firing rate during an input burst, and thus captures the shape of the wavefront, and (t) is the EPSP. In our model, both α(t) and (t) are positive valued for t > 0 (Methods), and act as low pass filters on the STDP rule. The full derivation for Eq 5 is provided in S1 Text. Note that the firing rate of the output neuron (R out in Methods) acts simply as a coefficient of in Eq 5 and hence plays a similar role to η by varying the rate at which synaptic strengths are modified. Amalgamating all terms in Eq 5, except for w(x), we have where κ(x) is the effective spike-location dependent plasticity rule, which incorporates the dynamics of input bursts and EPSPs: When a typical STDP rule (Fig 2A, top; see legend for parameters) is transformed into a spatial kernel (Fig 2A, bottom) by a wave traveling at 3 mm/s (the speed of spontaneous waves in the mouse cerebellum [8], or a 25°/s stimulus on the kitten retina using the visual angle to space conversion in Methods), the kernel extends over approximately 1 mm of input space. Note that κ(x) preserves the asymmetric shape of the STDP rule, but is low pass filtered by α(t) and (t).
A solution for w(x) as T ! 1 is more easily found in the frequency domain by taking the Fourier transform, such that convolution becomes multiplication: where '~' denotes the Fourier transform and k is the spatial frequency. Using w 0 to describe the The dominant spatial frequency, k*, lies wherekðkÞ has a global maximum. Inset:kðkÞ for a symmetric STDP rule (τ + = 20 ms, v = 3 mm/s). C. Evolution of the synaptic strengths, initialized with strength 0.5, from 500 input neurons to a single output neuron, under the influence of an asymmetric (left column) and symmetric (right column) STDP rule. Waves sculpt a periodic connectivity pattern into the synaptic strengths with a spatial frequency equal to k* (see B and Eq 5). Larger values for τ + (controls STDP rule width) or v (wave speed) produce lower spatial frequencies in the periodic pattern. Values for τ + and v are given between the columns. D-E. Log-log plots for the spatial frequency of the steady state connectivity pattern as a function of τ + and v, for both the symmetric (left panels) and asymmetric (right panels) STDP rules. Grey curves: predicted k*; blue: mean spatial frequency ± SEM of the final connectivity pattern in the spiking simulation; purple: mean spatial frequency ± SEM of w(x) in the numerical solution to Eq 6. D. Dependence of the dominant spatial frequency on τ + , keeping v = 3 mm/s fixed. E. Dependence of the dominant spatial frequency on v, keeping τ + = 20 ms fixed. F-G. Same as D and E, but for an integrate and fire output neuron model. initial synaptic strengths at time T = 0, the solution to Eq 8 forwðkÞ is We can reconcile Eq 9 with several previous studies of STDP by writingkðkÞ ¼ lðkÞ þ iðkÞ to explicitly express real and imaginary components: Here, the imaginary component gives rise to a spatial phase shift in w(x), a feature that has been exploited to drive a spatial redistribution of synaptic strengths by asymmetric STDP rules in several models of sequence learning [26,32], asymmetric shifts in hippocampal place fields [33,34], and the formation of direction selective cells [24,27,28].
The focus of our results, however, is the stability ofwðkÞ, which is determined by the real component, λ(k). We therefore consider waves that travel in both directions so that, for sufficiently small η, κ(x) is effectively symmetric and contains no imaginary component. This can be shown by integrating Eq 8 over two wave events, with the first wave traveling left to right, and the second wave traveling right to left (replacing v with −v). After the first wave, at T = T 1 , Usingwðk; T 1 Þ as the initial condition for the second wave, we have at wherekðkÞ is the complex conjugate ofkðkÞ, with conjugation resulting from the sign reversal in the wave speed. Expanding Eq 12, we havẽ wðk; T 2 Þ ¼w 0 ðkÞe ZlðkÞT 1 e iZðkÞT 1 e ZlðkÞðT 2 ÀT 1 Þ e ÀiZðkÞðT 2 ÀT 1 Þ where ΔT = T 1 = T 2 −T 1 is the time taken for one wave to cross the input layer. We can therefore write an approximation to Eq 9 for the special case in which waves travel in both directions in equal numbers:w If the wave direction were random, instead of alternating, Eq 14 would be valid only if a large number of waves traversed the input layer during the interval ΔT/η.
Thus, for any k such that Re½kðkÞ < 0,wðkÞ will be stable and decay to zero, whereas for any k such that Re½kðkÞ > 0,wðkÞ will become unstable and grow exponentially. For the STDP rules used in this study (Fig 1C),kðkÞ has a positive valued maximum at k Ã (Fig 2B). Therefore, if the synaptic strengths are initially random, such that the expected form ofw 0 ðkÞ is flat, thenwðk Ã Þ will be the fastest growing eigenmode of the synaptic strengths, andwðkÞ will asymptotically approximate a delta function, δ(k Ã ). If only a single spatial frequency dominates wðkÞ, w(x) will be sinusoidal. Thus, our derivation predicts that, for STDP rules like those reported experimentally, the synaptic strengths will develop a connectivity pattern that is periodic in space with a period of 1/k Ã , under the influence of traveling waves.
A similar model of pattern formation describes the development of ocular dominance columns in primary visual cortex [35], in which it is proposed that short range excitatory and long range inhibitory lateral interactions give rise to the spatially periodic dominance of eye specific afferents across primary visual cortex. In our model, it is the mapping of the STDP rule onto space that provides the lateral interactions necessary for pattern formation. More generally, the solution derived above is analogous to pattern forming solutions that result from Turing instabilities in reaction-diffusion systems [29], whereby the initially homogeneous distribution of synaptic strengths becomes unstable, allowing the fastest growing eigenmode to dominate the resulting pattern. In systems composed of a diffusible activator and an inhibitor, Turing instabilities frequently occur when the inhibitor diffuses over greater distances than does the activator [30]. Here, the decay constants of the positive and negative STDP lobes bear similarities to the length scales of diffusion in reaction-diffusion systems. Thus, STDP rules with narrow windows for strengthening and wide windows for weakening are good candidates for pattern formation in neural circuits that support traveling wave activity patterns.

Stability of the connectivity pattern
If the synapses are unbounded, the STDP rules used here will always yield a sinusoidal connectivity pattern when driven by traveling waves, given sufficient time, and the mean synaptic strength will always be zero (the mean of a sine wave). This is far from a physiologically plausible scenario: without bounds, the synaptic strengths would tend towards ±1. When bounds are imposed, however, care is needed to formulate the STDP rule, so that the synaptic strengths do not reach the upper or lower bound before the dominant eigenmode at k = k Ã takes hold. This can be achieved by ensuring that there is not too strong a bias for either synaptic weakening or strengthening in the STDP rule: where B L and B U are, respectively, the magnitudes of the lower and upper bounds to the STDP bias. Because synaptic strengths are constrained to be positive, it must be that B L > B U . In other words, the range of biases for synaptic weakening that will yield stable periodic patterns is greater than that for strengthening. Nevertheless, it is important to note that we need not have B U < 0. This contrasts with previous studies that utilized different input activity patterns [16,17,19], and that required a bias for synaptic weakening by setting R 1 À1 dDt KðDtÞ < 0 to stabilize the connectivity pattern. That is, in our model, it remains possible for periodic patterns to emerge even if there is a bias for synaptic strengthening in the STDP rule, so long as the synaptic strengths are not pushed to the upper bound before the dominant spatial frequency takes hold. The stability of the connectivity will also be sensitive to the learning rate, which is scaled by η in Eq 6, and the initial conditions of the synaptic strengths. For example, we later explore formulations ofkðkÞ for which there are multiple peaks that are similar in amplitude, or a single, broad peak. The stochastic dynamics introduced by spiking neurons may therefore move the slower dynamics of plasticity along a spectrum of eigenmodes, resulting in more aperiodic connectivity patterns. We later introduce a robustness measure to quantify the periodicity of a connectivity pattern and that takes these potential features ofkðkÞ into account.

Spatial pattern formation in a spiking neuron simulation
To test the predictions of the analytical derivation above, we examined whether a periodic connectivity pattern would develop as a result of traveling waves and STDP in a simulated network of linear Poisson, spiking neurons. The simulation captures the architecture and function of the network illustrated in Fig 1A, consisting of a 1D layer of input neurons that all synapse onto a single output neuron. During simulations, synaptic strengths were modified by one of two distinct forms of STDP rule (see Methods) that have been reported in the literature, one that is asymmetric (Eq 21) [13,17,36] and the other symmetric (Eq 22) [37,38] in Δt. Spiking activity is generated by a traveling wavefront that moves back and forth along the input layer, eliciting EPSPs and spikes in the output neuron. We varied two parameters in the simulation: i) the temporal window of K(Δt), which we control with the decay time constant for synaptic strengthening, τ + , and ii) the wave speed, v. These parameters modulate the shape of κ(x) and, therefore, the spatial frequency of the predicted periodic connectivity pattern.
As predicted by our derivation (Eq 9), the synaptic strengths in the simulation indeed developed a periodic connectivity pattern (Fig 2C) in the presence of traveling waves for both the asymmetric and symmetric rules. The periodic pattern was robust over a range of STDP time constants (τ + = 20, 30, 40, 50, 60 and 70 ms), keeping the wave speed constant (v = 3 mm/s), and for a range of wave speeds (v = 1, 2, 3, 4, 5, and 6 mm/s), keeping τ + constant (τ + = 20 ms). As predicted, increasing τ + or v caused the period of the connectivity pattern to increase ( Fig  2C). To test the accuracy with which Eq 9 predicts the resulting spatial frequency over this range of parameters, we computed the power spectrum of the steady state synaptic strengths (mean subtracted) to which a Gaussian curve was fit. The spatial frequency of the connectivity pattern was taken to be the centre of the fitted Gaussian. We repeated simulations sixteen times for each set of parameters, using a different seed for the random number generators. In Fig 2D and 2E, we compare the measured spatial frequencies (blue circles) with the predicted spatial frequencies, k Ã (grey curves), as a function of τ + ( Fig 2D) and v (Fig 2E), where the predicted values were found by numerically computing Re½kðkÞ and determining the spatial frequency at its peak. We measured the accuracy of our predictions by computing the coefficient of determination, R 2 , between the logarithms of the predicted and measured spatial frequencies. For each panel in Fig 2D and 2E, R 2 > 0.85. Thus, the assumptions made to derive Eq 9 appear to maintain a veritable description of the noisy dynamics in the simulation over the range of parameters tested. In addition to the simulations, we verified that solutions obtained by numerically integrating Eq 6 (incorporating nonlinearities such as hard bounds to w(x), see Methods) also produced periodic connectivity patterns with spatial frequencies that matched predictions (magenta circles, Fig 2D and 2E, R 2 > 0.92).
To investigate whether nonlinearities in the postsynaptic response influence the outcome of the connectivity patterns, we replaced the linear output neuron with a leaky integrate and fire (LIF) neuron (see Methods) that modeled absolute and relative refractory periods of 2 ms and 5 ms, respectively. These simulations yielded a similar relationship between the spatial frequency of the connectivity pattern, wave speed and STDP time scale (Fig 2F and 2G). This similarity might be expected, as the wave correlations extend over relatively long time scales and thus smooth out any ripples in the correlation function, C(x,Δt), that would be introduced on the short time scales of the refractory period or EPSP rise time. More noticeable differences may, however, be observed if wavefronts elicited very short bursts or only single spikes.
The development of periodic patterns in synaptic connectivity could have wide applications throughout the nervous system of many species, particularly because of the ubiquity of both traveling waves and STDP. However, it is first necessary to determine the extent to which pattern formation is influenced 1) over the range of space and time scales observed in biology, 2) in the presence of noise, and 3) by 2D waves. We explore these issues in the following sections.

Pattern formation over a wide range of spatial and temporal scales
In the previous section, we found that the wave speed and STDP time scale are key parameters that determine the spatial frequency of periodic connectivity patterns. Traveling waves in different areas of the brain are characterized by wave speeds that span at least two of orders of magnitude, from slow retinal waves with speeds on the order of 0.1 mm/s [39] to fast cortical waves with speeds reaching 17 mm/s [9]. On the other hand, time scales for STDP are typically 10-100 ms [40], but time scales on the order of seconds are predicted to be relevant to retinal waves [41]. We consider in the Discussion how the theoretical results above might apply to these different biological circuits. To do this, we first obtain a more complete picture of the spatial scales of pattern formation predicted by the theory over a wide range of wave speeds and STDP time scales.
We calculated the landscape of spatial frequencies as a function of τ + and v by numerically computingkðkÞ, and finding the spatial frequency, k Ã , that maximizes its real component. The k Ã landscapes for both the asymmetric and symmetric STDP rules (Fig 3A and 3B, respectively) reveal a remarkably similar dependence of the predicted spatial frequency on τ + and v, which is perhaps not surprising given the simple exponential functions underlying the STDP rules (Eqs In all plots, solid black contours denote spatial frequencies equal to 10 raised to integer exponents. A. Predicted spatial frequencies for an asymmetric STDP rule as a function of v and τ + . White triangle: when τ + is longer than the burst duration, scaling τ + or v by the same amount has the same effect on the spatial frequency. Green triangle: when τ + is shorter than the burst duration, the burst duration has a noticeable influence on the spatial frequency. Red rectangle: wave speeds recorded in the cerebellum, ventral cortex and hippocampus of neonatal mice occupy this region for STDP rules with τ + % 20 ms. B. Landscape of predicted spatial frequencies for a symmetric STDP rule. 21 and 22). In particular, for a large region of parameter space, scaling either v or τ + by a factor, f, simply scales k Ã by 1/f. We illustrate this scaling feature using triangles with two equal sides aligned to the axes, as multiplication in linear space is equivalent to addition in logarithmic space. Starting at one iso-frequency contour (1 cycle/mm, bottom left vertex of white triangle in Fig 3A), a constant step along the v-axis moves k Ã to the same iso-frequency contour (0.04 cycles/mm) as does an equal step along the τ + -axis. Because the input burst acts as a low pass filter on the STDP rule, this relationship does not continue (green triangle, Fig 3A) when τ + becomes shorter than the burst duration (here 0.1 s) of the input neurons. Thus, the burst duration has a strong influence on k Ã at time scales longer than τ + , which is particularly the case during retinal waves when burst durations are often as long as 3 s [42,43]. Addressing the impact of the burst duration is the focus of the next section.

Impact of burst duration on pattern formation
The warping of the spatial frequency landscapes above shows that the burst duration also plays a role in pattern formation. Furthermore, when correlated activity patterns comprise long burst durations, STDP rules with short time scales struggle to extract information from the correlations that might be relevant for circuit development [20,41,44]. Here, we further examine the influence of burst duration on pattern formation in our model by carrying out simulations over the range of burst durations that are observed for different types of traveling wave.
Waves in this set of simulations traveled with a fixed speed of 3 mm/s, and we used the asymmetric STDP rule with a fixed decay time of t + = 20 ms. In Fig 4A, we show examples of the evolving connectivity pattern for different burst durations. For bursts lasting 0.03 s, the connectivity pattern had only a very weak periodic structure ( Fig 4A, left panel). The power spectrum of connectivity patterns when the burst duration was 0.03 s, averaged over repeated trials, is shown in Fig 4B (blue). By normalizing the spectrum to the power at its peak, it is clear to see that power is spread over a broad range of spatial frequencies. Burst durations on the order of a few hundred milliseconds produced more distinct periodic connectivity patterns ( Fig 4A, middle panels), with power concentrated around the peak in the power spectrum ( Fig  4B, orange). However, in simulations with burst durations of 1 s or longer, the connectivity pattern became disordered ( Fig 4A, right panel), with some power concentrated at the lowest spatial frequencies and otherwise spread evenly across higher frequencies (Fig 4B, black).
To summarize the robustness with which a periodic connectivity pattern is produced under different conditions, we define robustness to be the ratio of the power spectrum amplitude at its peak to the total power in the discretized power spectrum (see Methods). All robustness measures are plotted in comparison to a reference value, which we take to be the mean robustness when parameters yielded a clear periodic connectivity pattern. Here, simulations with a 0.1 s burst duration were used as the reference. The robustness is plotted in Fig 4C for simulations using burst durations from 0.03 s to 5 s (black circles). In agreement with the shapes of the power spectra in Fig 4B, the robustness is relatively high for burst durations in the range 0.1-0.5 s, and relatively low outside of this range. Because non-linearities in the simulation, such as bounded synaptic strengths, allow for the building of multiple spatial frequencies into the connectivity pattern, we examined whether the concentration of power in Re½kðkÞ at the dominant spatial frequency, k Ã , could explain the relationship between the burst duration and robustness. That is, we numerically computed Re½kðkÞ over the full range of simulated burst durations, from which we calculated the ratio of power at k Ã to the total power (see Methods). We plot this theoretical measure of robustness in Fig 4C ( . For a 0.03 s burst, power is distributed across a wide range of high frequencies, whereas, for a 3.0 s burst, most of the power is concentrated in the negative dip at the lowest frequencies. However, for 0.3 s bursts, power is concentrated between these two extremes, around the dominant spatial frequency of k Ã . This trend is not specific to the choice of α(t). In S1A Fig, we show Re½kðkÞ for which α(t) is modeled using an alpha function.
Despite the poor robustness at short burst durations, the spatial frequencies that were measured from the connectivity patterns were well matched to the predicted k Ã for burst durations between 0.03-0.5 s ( Fig 4E). However, for burst durations exceeding 0.5 s, the connectivity patterns no longer yielded spatial frequencies that matched the theory. This is likely to be due to the similar amplitude of multiple peaks inkðkÞ for longer burst durations (for example, the black curve in Fig 4D), compounded by the overall loss in robustness ( Fig 4C). Butts & Rokhsar [41] have shown that waves with long burst durations provide more information that is relevant for refinement when the plasticity rule has a time scale much longer than is typically seen for STDP. It might therefore be possible to rescue pattern formation for waves with long burst durations by using a longer time scale for the STDP rule. To examine this possibility, we . The negative lobe in the black curve extends beyond the horizontal axis and has been cut for clarity. E. Spatial frequencies of connectivity patterns that developed in simulations as a function of the burst duration (black circles and error bars: mean ± SEM). Grey line: predicted spatial frequencies. The sawtooth fluctuations at longer burst durations correspond to the emergence of new dominant peaks inkðkÞ, which is compressed towards lower spatial frequencies as the burst duration is increased. F. Theoretical robustness as a function of burst duration for asymmetric STDP rules with τ + = 20 ms (pink) and τ + = 50 ms (green). computed the fraction of power atkðk Ã Þ over the range of simulated burst durations for a STDP rule with τ + = 50 ms (Fig 4F, green). The wider STDP rule exhibits a greater concentration of power at k Ã over a much wider range of burst durations, when compared with the STDP rule with τ + = 20 ms, and the same trend is seen for other α(t) kernels (S1B Fig). This includes bursts exceeding 1 s in duration. For burst durations as long as 3 s, which are common for retinal waves, STDP rules with even longer time scales would be expected, according to our model. Thus, our model supports the hypothesis of Butts & Rokhsar [41], and provides a new approach for estimating the required time scale of STDP with which retinal waves may refine developing neural circuits.
Additional contributions to wave-related correlations can arise from the noisy mechanisms of wave and spike generation during early development, which may deliver waves in quick succession as well as generate non-wave related input spikes. The extent to which our analytical results can be applied to biological systems therefore depends on their robustness to these additional contributions to C(x,Δt). In the following sections, we seek to understand how the presence of noise and multiple waves may impact upon the development of periodic connectivity patterns.

Sensitivity of the connectivity pattern to noise
Thus far, we have tested theoretical predictions using idealized wavefronts. In reality, spontaneous and sensory driven waves exist among continuous background activity and travel with variable wave speeds. Here, we use simulations to establish the sensitivity of pattern formation to these sources of noise. To examine the effect of variable wave speed, we conducted a set of simulations in which a new speed was assigned to each wave from a lognormal distribution, with a mean of 4 mm/s and standard deviation (SD) σ v . We imposed a minimum speed of 0.05 mm/s so that waves would not take too long to traverse the input layer. Examples of the evolving connectivity pattern in these simulations are provided in Fig 5A and, in Fig 5C, we plot the mean power spectra of the final connectivity patterns for σ v = 0.0, 2.0 and 4.0 mm/s. An unexpected feature of the resulting connectivity patterns was a decrease in the spatial frequency with increasing σ v (orange: σ v = 0.0 Hz, cyan: σ v = 2.0 Hz, black: σ v = 4.0 Hz). In these simulations, the developing connectivity pattern may experience a greater influence from the faster waves that recruit more input neurons per unit time and thus drive higher output firing rates. This would correspond to our earlier results in Fig 2, for which we showed that lower spatial frequencies are associated with higher wave speeds. In Fig 5E, we plot the robustness of the connectivity pattern to variation in wave speed, where we have referenced robustness to the case when σ v = 0. Despite the variation in wave speed, which corresponds to there being a spectrum of spatial frequencies impressed on the network, the connectivity retains a reasonably robust periodic structure.
We next tested the robustness of periodic connectivity patterning to the presence of background spiking noise, which diminishes the relative contribution of wave-activity to C(x,Δt) in Eq 3. To implement background noise in the simulation, input neurons that were not participating in a wave fired spontaneous spikes with a firing rate of R 0 , which was varied between simulations (see Methods). Examples of the evolving connectivity patterns are shown in Fig 5B  (top row) for different background firing rates. Mean power spectra of the final connectivity patterns when R 0 = 0.0 (cyan), 1.0 (orange) and 2.0 Hz (black) are plotted in Fig 5D, and show little variation in the dominant spatial frequency when R 0 < 2.0 Hz. However, the robustness degraded substantially when R 0 ! 2.0 Hz (black circles in Fig 5F, reference case R 0 = 0 Hz). In these simulations, a new wave traversed the input layer every 10 s, but the burst duration was only 0.1 s. Therefore, the ratio of wave-related spikes to background spikes, rather than the background rate, might better parameterize when the connectivity pattern will be robust. In this case, periodic patterns lacked robustness when the ratio went below approximately 1:4. Accordingly, the robustness could be significantly enhanced for almost all background rates by increasing the firing rate during a wave burst to 100 Hz. In this case, the robustness degraded substantially when R 0 ! 4 Hz (Fig 5B, bottom row; Fig 5F, grey circles). During early, spontaneous waves, low background firing rates are typical, and there is good reason to believe that background firing rates are low in early stages of sensory processing as well [45]. We review more evidence for this in the Discussion.
Spike-spike correlations induced by the wavefront, rather than time averaged firing rates, are the driving force behind periodic patterning in our model. Additional sources of spikespike correlations may therefore disrupt periodic patterning. We conducted a set of simulations in which inputs experienced instantaneous correlations with their neighbors (see Methods) while also varying the background firing rate. With non-zero background rates, increasing the correlation strength does indeed degrade the robustness of the periodic connectivity patterns (Fig 5F, orange circles), referenced to simulations with no local correlations and R 0 = 0.0 Hz. In the absence of background spiking, however, the additional local correlations act to enhance the robustness. Typically, correlations are not instantaneous but decay over time [46][47][48][49] and thus would have a reduced impact on plasticity due to the decaying amplitude of the STDP rules around Δt = 0 [18].
Sensitivity of the connectivity pattern to multiple, non-isolated waves The intervals between consecutive spontaneous waves can be very variable, and inter-wave intervals (IWIs) cover a broad range, from approximately 100 ms in the cerebellum [8] to tens of seconds in the retina [39,50]. IWIs during sensory driven waves may also be highly variable, and are likely to match the temporal pattern of naturally occurring stimulus features to which a population of input neurons are tuned. During vision, for example, multiple luminance contours can be cast across the retina in quick succession when tracking a moving object. In the following, we investigate how the presence of multiple waves, which simultaneously drive spiking in the output neuron, impact the development of the connectivity pattern. In so doing, we test our assumption in the above derivation that waves be sufficiently isolated in time.
We ran a set of simulations in which waves were generated with an approximately lognormal distribution of IWIs and varied the mean IWI, μ IWI , between simulations (see Methods). At a speed of 4 mm/s, waves took 2.5 s to traverse the input layer. Thus, to ensure that multiple waves were present for the majority of the simulation, we set μ IWI 2.0 s across the set of simulations.
Despite multiple waves driving the output cell, periodic connectivity patterns emerged in the simulations, examples of which are shown in Fig 6A (left and center panels). Power spectra of the connectivity patterns had distinct peaks near the predicted spatial frequency (0.91 cycles/mm) for μ IWI as short as 0.5 s (Fig 6C). The connectivity pattern began to degrade for μ IWI < 0.5 s, and had little structure when μ IWI 0.2 s (Fig 6A, right panel), with power spread broadly across the spectrum (Fig 6C). The loss of robustness with decreasing μ IWI is summarized in Fig 6E (black circles), where robustness is referenced to simulations in which a constant IWI of 5 s was used, thus ensuring that only one wave was present at a time.
The loss of robustness with decreasing μ IWI may be related to our requirement in the above derivation that waves be sufficiently isolated in time, with ΔT > 2K, where ΔT can be interpreted as the IWI and K as half the temporal width of the STPD rule. Although K is not precisely defined, values in the range 0.1-0.3 s would be reasonable for the STDP rule used here, for which the decay time of the longer, negative lobe was τ − = 0.04 s. An alternative explanation might be that successive waves with very short IWIs are not very different from single waves with long burst durations, as far as the timing precision of the STDP rule is concerned. Because the burst duration also influences robustness, these two possible factors would be difficult to disentangle, and efforts to do so go beyond the scope of this study. The important qualitative result is that, above a lower limit to the mean IWI of only a few tenths of a second, the succession of randomly occurring waves does not greatly degrade the robustness of periodic connectivity patterns.
We next asked whether a constant IWI between waves, which would contribute a strong frequency component tokðkÞ, would lead to a similar degradation in robustness with decreasing IWI. Because wavefronts of activity in the retina readily track moving luminance edges [7], regular waves could correspond to stimulation of the retina by luminance gratings, as used in experiments studying the development of orientation and direction selectivity [51]. We ran simulations in which waves traversed the inputs in a periodic fashion, with a constant IWI, but otherwise used the same parameters as for the simulations with irregular waves. Examples of the periodic patterns that developed for different IWIs are provided in Fig 6B, and their power spectra are shown in Fig 6D. It is clear that regular waves built robust periodic patterns, even with the shortest IWIs. Our robustness measure confirmed that this was the case across all of the IWIs tested ( Fig 6E, grey circles). In addition to the enhanced robustness, a distinct feature of the connectivity patterns with constant IWIs is the shift towards higher spatial frequencies for short IWIs (Fig 6D, power spectra for IWIs of 0.2 s and 0.15 s).
The emergence of robust periodic connectivity patterns, with increased spatial frequencies for fixed IWIs 0.2 s, appears at odds with our theoretical predictions. However, these features can be easily accounted for by considering how κ(x) is constructed in Eq 7. To represent the periodic input bursts elicited by multiple waves, we replace α(t) in Eq 7 with a p ðtÞ ¼ R 1 À1 dt 0 P n dðt 0 À ðn Â IWIÞÞaðt À t 0 Þ, the convolution of α(t) with a Dirac comb. Eq 7 expresses κ(x) as the convolution of K v (x), α p (x/v), α p (−x/v) and (x/v), and can be solved by taking its Fourier transform, whereby convolution in the real domain is equivalent to multiplication in the frequency domain. In Fig 6F, we plot the two functions that have the greatest influence in our simulations. In black is the Fourier transform of the STDP rule,K v ðkÞ. The remaining curves are the Fourier transforms of the input firing patterns,ã p ðvkÞ (with the zeroth frequency component removed), for three cases: a single wave in isolation (effectively, IWI = 1 grey curve) and periodic waves with IWI = 1.0 s (blue) and IWI = 0.15 s (pink). Note that the Fourier transform of the EPSP has little effect in this calculation because its time constant is much shorter than that of the other functions. The productK v ðkÞ jã p ðvkÞj 2 therefore yields a good approximation ofkðkÞ, and is drawn for the three cases ofã p ðvkÞ in Fig 6G. Because the Fourier transform of a Dirac comb is another Dirac comb,ã p ðvkÞ produces sharp peaks inkðkÞ (pink and blue) at spatial frequencies that may be higher (pink) than the peak in the case of isolated waves (grey). The large amplitude and narrow peaks inkðkÞ, when IWI = 0.15 s and 0.2 s, explain the dominance and robustness of a single spatial frequency in the respective connectivity patterns, and the positions of these peaks explain the shift of the power spectra toward higher spatial frequencies in Fig 6D. The strong dependence of the spatial frequency on the IWI for short IWIs only, as depicted in Fig 6F and 6G, suggests a critical value for the IWI. We can relate the critical value to the time taken for a wave to travel a distance equal to one cycle of the periodic pattern: IWI crit = 1/ vk Ã % 0.27 s. That is to say, when IWI < IWI crit , the spatial frequency with which waves are lined up along the input layer exceeds the dominant spatial frequency of κ(x) for a wave in isolation (S2 Fig). This interpretation fits well with our assumption in the derivation that waves be sufficiently isolated in time (and therefore isolated in space). Moreover, the existence of a critical IWI provides a means for quantifying K, and thus provides a possible explanation for the degradation in robustness when μ IWI 0.3 s in the simulations with irregular waves.
Traveling waves are instructive for shaping receptive fields via spiketiming dependent plasticity The filtering properties of a RF result from feedforward and recurrent synaptic inputs to the neuron, as well as its intrinsic cell properties. In this section, we used simulations of spiking neurons to examine whether the interaction of traveling waves and STDP rules, which can impose spatially periodic connectivity patterns on a uniform field of synaptic strengths (shown above), might be a plausible mechanism for shaping the feedforward component of RFs. As such, a RF in our model refers to the spatial pattern of synaptic strengths impinging on the output neuron from a small region of the input layer. During development, coarse RFs are thought to be set up by molecular cues before activity dependent refinement takes place [52]. Therefore, at the beginning of each simulation, we set synapses from the center of the input layer to maximum strength and the remaining synapses to zero, thus endowing the output neuron with an initial RF having a diameter, RF 0 . Despite this restricted initial arrangement of synaptic strengths, we found that a periodic pattern would nevertheless emerge across the entire input array (S3 Fig). This outcome is not typical of the brain, as topographic maps and the limited size of axonal arbors and dendritic trees restrict the spatial distribution of synaptic inputs to a neuron. We therefore applied an arbor function (see Methods) that was centered in the middle of the input layer and spanned a region greater than RF 0 . Synapses that were outside the arbor were disconnected from the output neuron and prevented from strengthening. A schematic of this new network architecture is provided in Fig 7A. In the majority of the following simulations, we set RF 0 = 0.8 mm, which approximates an area of retina that corresponds to RF sizes in developing areas of the visual system, including the cat primary visual cortex [53,54], the lateral geniculate nucleus (LGN) in ferrets [55,56], and the superior colliculus (SC) in adult mice that have had disrupted retinal waves. The conversion between RF size and retinal distance is described in the Methods. Furthermore, in the remaining sections, we increase the bias By interpreting the output neuron RF as a single cycle in the type of periodic connectivity patterns obtained above, we are able to examine how properties such as the wave speed and the STDP time scale may shape RF development by controlling the characteristic wavelength, 1/k Ã , of the connectivity pattern. Depending on the size of 1/k Ã with respect to RF 0 , one of three modes of RF modification were observed. In one set of simulations, the RF became larger ( Fig  7B, left), smaller (Fig 7B, middle) or split into subfields (Fig 7B, right) with progressively slower wave speeds and thus shorter characteristic wavelengths. The same modes of RF modification were achieved by decreasing the STDP time scales (Fig 7C) and holding the wave speed constant. If instead of changing k Ã we increased RF 0 , keeping the wave speed and STDP rule fixed, the RFs were similarly modified ( Fig 7D). Thus, characteristic wavelengths that are long relative to RF 0 caused the RF to grow, shorter wavelengths caused the RF to shrink, and even shorter wavelengths enabled multiple cycles of the periodic pattern to form within the arbor.
Our results suggest that 1/k Ã in fact determines the size of the final RF, or each subfield. To verify this, we varied RF 0 from 0.2 mm to 3.2 mm between simulations, with the wave speed and STDP rule fixed, and looked for the convergence of RF sizes, which we measured as the number of synapses with strengths > 0.5. To accommodate RF 0 , the width of the arbor function was set to the larger of 0.8 mm or RF 0 +0.4 mm. RF sizes converged to one of three sizes (Fig 7E), depending on whether the RF maintained a single field (purple), or split into two (yellow) or three (blue) subfields. The larger RF 0 , the more subfields that emerged. Regardless of RF 0 , the mean final size of each subfield was approximately the same. Integer multiples of the mean subfield size (measured from RFs with three subfields) are drawn on the right of Fig 7E to illustrate this fact. Small biases were evident, however, whereby larger values of RF 0 tended to give rise to larger RFs or subfields. This considered, the final number of strong synapses were easily distinguishable between cases in which one, two or three subfields developed. Thus, the size and shape of the final RF tightly corresponds with the wave, STDP and initial RF properties.
A quantitative relationship between periodic patterning and refinement of receptive fields We showed above, in Fig 7B and 7C, that different wave and STDP parameter combinations could yield the same mode of RF modification, including RF expansion, contraction, or splitting. Here, we focus on RFs that are made to contract, which is particularly relevant to the refinement of topographic maps such as the retinotopic map in superior colliculus, and examine the range over which the wave speed and STDP time scale achieves this mode of modification. To do this, we numerically integrated Eq 6 (see Methods) for the wide range of wave speeds and STDP time constants used to generate Fig 3, holding RF 0 = 0.8 mm constant. An example of w(x) during the numerical integration of Eq 6, starting with an initial RF, is shown in Fig 8A, for which parameters matched those in the central panels of Fig 7B and 7D. Using the solutions for w(x), we constructed a phase diagram for the size and shape of RFs that developed, which we call the refinement phase space.
In the refinement phase space, shown in Fig 8B, there is a region (shaded blue) corresponding to v and τ + parameters that caused the RF to contract, with darker shades of blue indicating greater contraction, measured as a percentage of RF 0 . In the region labelled α (unshaded), at which κ(x) has higher characteristic spatial frequencies, RFs split into subfields. In the region labelled β (unshaded), at which spatial frequencies are lower, RFs expanded. Overlaid are the iso-frequency contours of the predicted spatial frequencies, k Ã , as in Fig 3. The boundaries between each mode of RF modification closely follow the iso-frequency contours, with greater RF contraction (darker shades) occurring at higher spatial frequencies. Because the outcome of the final RF depends on the initial RF size, we recomputed the refinement phase space, setting RF 0 = 0.44 mm, which just exceeded the smallest size of RFs after contracting in solutions when RF 0 = 0.8 mm. For RF 0 = 0.44 mm, the region of phase space in which RFs contracted (colored purple) corresponded to slower wave speeds and shorter STDP time constants, i.e. higher spatial frequencies in κ(x). Taken together, these results confirm that smaller initial RF sizes require higher spatial frequencies in κ(x) to achieve a particular mode of RF modification. The area of the blue and purple regions of phase space also show that, in order for RFs to be refined to a smaller size by the asymmetric STDP rule, the maximum possible wave speed can be just over a factor of two times greater than the minimum possible wave speed. At least the same range in τ + can also enable RF contraction. Thus, the mechanism of RF contraction Relationship between periodic patterning and receptive field refinement. A. Example of RF refinement obtained by numerically solving Eq 6. Vertical teal lines denote the arbor function boundary. B. Refinement phase space for RFs as a function of wave speed and STDP time constants. Colored regions: areas of phase space in which RFs maintain a single field that is refined to smaller sizes; blue: RF 0 = 0.8 mm; purple: RF 0 = 0.44 mm. Shading indicates the total change in RF size as a percentage of RF 0 . Region α corresponds to higher spatial frequencies in which RFs split into subfields. Region β corresponds to lower spatial frequencies in which RFs expand. Overlaid are the iso-frequency contours of Fig 3. Solid contours correspond to spatial frequencies with integer exponents of 10 (lower left to upper right: 10 cyc./mm, 1 cyc./mm and 0.1 cyc./mm). C. Incremental changes in wave speed can continually refine RFs while maintaining a single subfield structure. Shown is the development of a RF, averaged over 16 trials, during which the wave speed is first set to 4 mm/s and then decreases to 2.5 mm/s, after which the RF decreases in size. Increasing the wave speed back to 4 mm/s returns the RF to its previous refined size. Purple trace: plot of the mean number of synapses with strength > 0.5 as an indicator of RF size. D. Refinement phase space using a symmetric STDP rule, for which RF 0 matches that for the blue region in B. Numbers next to the iso-frequency contours indicate their spatial frequency in cycles/mm. E. Refinement phase spaces, using the asymmetric STDP rule, for different biases towards synaptic weakening. Left: A − = 0.50, which corresponds to R 1 1 dDt KðDtÞ ¼ 0. Centre and right: A − = 0.55 and 0.60, respectively, which correspond to described here allows for some variation in STDP properties between synapses, variation in wave speed, or even fluctuations in both of these phenomena over time.
An interesting corollary of these results is that changes in wave speed and STDP properties can follow changes in RF size, thus allowing continued contraction beyond that of a fixed v and τ + combination. To demonstrate this point, we ran a simulation using the same parameters as those used for the central panels in Fig 7B and 7D. However, once the RF had refined to a smaller size, we decreased the wave speed from 4 mm/s to 2.5 mm/s. Note that, according to Fig  8B, v = 2.5 mm/s and τ + = 20 ms would cause a RF with RF 0 = 0.8 mm to split into subfields, as this point in the phase space lies to the lower left of the blue shaded region. However, as shown in Fig 8C, an initial period of contraction with v = 4 mm/s set up further contraction, after a decrease in wave speed to v = 2.5 mm/s, without splitting the RF into subfields. Increasing the wave speed back to 4 mm/s at a later time caused the RF to expand and return to its previous size. Analogously, changes in the time constants of STDP can achieve the same effect.
The phase diagram corresponds very well with the simulation results in Fig 7, confirming our hypothesis that κ(x) and RF 0 jointly determine the final shape and size of the RF. This was not unique to the asymmetric STDP rule. In Fig 8D, we show a similar correspondence between refinement and the iso-frequency contours for a symmetric STDP rule.
Tuning the characteristic wavelength by itself was not sufficient to achieve RF contraction. In addition, both asymmetric and symmetric STDP rules required a bias towards synaptic weakening. In Fig 8E, we plot the refinement phase space of a 0.8 mm RF for three asymmetric STDP rules that differed only in A − , which scaled the amplitude of the negative lobe in the rule. When A − = 0.5, we have R 1 À1 dDtKðDtÞ ¼ 0, such that synapses no longer compete for connections with the output. Consequently, increasing the dominant spatial frequency of κ(x) would not decrease the size of the RF before causing it to split into subfields (left panel). However, by increasing A − (A − = 0.55, middle panel; A − = 0.6, right panel), such that R 1 À1 dDtKðDtÞ < 0, RF contraction was possible. Furthermore, increasing A − moved the iso-frequency contours, making RF contraction possible at lower characteristic spatial frequencies and, therefore, larger v and τ + values. As such, for a given v and τ + pair (red circles in Fig 8E), greater contraction was achieved with increased bias for synaptic weakening.
Equipped with a basic understanding of the relationship between periodic patterning and RF modification, we examine how more realistic traveling waves and STDP impact the development of RFs in two dimensions in the following sections.

Refinement of receptive fields in two spatial dimensions
In two dimensions, the variables x and v in Eq 6 can now be considered as 2D vectors, x = (x,y) and v = (v x ,v y ), in the x-y plane. As such, we expected the connectivity pattern to adopt the characteristic spatial frequency, k Ã = (k x ,k y ), and thus for the type of RF modification along a particular axis to depend on the direction of wave propagation. To test this, we generated wave stimuli in a 2D input layer consisting of 64×64 units arranged on a square lattice. Plane waves were generated in a similar fashion to the 1D scenario, moving with a constant speed in alternating directions, from one side of the input layer to the opposite side with a wave speed of 4 mm/s. The single output neuron was provided with an initial, circular RF 0.8 mm in diameter (Fig 9A, left panel), and synapses were modified by an asymmetric STDP rule (τ + = 20 ms). When plane waves travelled along just one axis, the RF contracted along the same axis (Fig 9A; top panels: waves travel along the horizontal axis; bottom panels: waves travel along the vertical axis), corresponding with the 1D results above (Fig 8B). However, the RF grew along the orthogonal axis until it spanned the full extent of the arbor. This corresponds with an effectively infinite wave speed along the orthogonal axis, resulting in a periodic pattern with a spatial frequency of zero.
Exposing the network to a reduced set of wave directions in this way can be used to build RFs that exhibit periodic structure along the same directions. In S4 Fig, we provide examples of RFs that developed when waves traveled in both directions along one, two and three axes. Of particular interest are RFs that exhibit multiple parallel subfields (S4A Fig), which bare several similarities with oriented simple cell RFs in primary visual cortex. In the Discussion, we describe in more detail how periodic patterning might be applied to realistic simple cell RFs. When waves travelled in 16 possible directions, equally spaced around the compass and in a random sequence, the final RF maintained an approximately circular shape (Fig 9B, left). However, using the same wave and STDP properties as for Fig 9A, the RF area increased (Fig 9C,   Fig 9. Refinement of receptive fields with 2D plane waves. A. Refinement depends on the direction of wave propagation. If an initial RF is circular (left), then waves traveling along a horizontal direction will refine the RF along the horizontal axis only (top row). The RF necessarily expands to the arbor boundary along the orthogonal axis. Likewise, waves traveling along a vertical direction refine RFs along the vertical axis only (bottom row). Increasing time from left to right. Arrows: direction of wave propagation. Solid teal circle: arbor boundary. Dashed green circle: boundary of the initial RF. B. Final RFs that were exposed to plane waves traveling randomly in 16 possible directions. Left: using the same wave speed (4 mm/s) and STDP parameters (τ + = 20 ms, A − = 0.55) as in A. Centre: using a slower wave speed (3 mm/s) to lower the characteristic spatial frequency of κ(x). Right: using the faster wave speed (4 mm/s) but shifting refinement phase space to lower frequencies by increasing the bias for synaptic weakening (A − = 0.60). C. Summary of the development of RF sizes for each of the conditions in B, averaged over three trials for each condition. dark blue) because contraction along the axis parallel to the wave direction was outweighed by expansion along the orthogonal axis. RF expansion could be counteracted by changing parameters that enhance contraction, as explained in the results of Fig 8. Thus, increasing k Ã by decreasing the wave speed to 3 mm/s (Fig 9B, middle), or increasing the amplitude for synaptic weakening in the STDP rule from A − = 0.55 to A − = 0.6 ( Fig 9B, right), caused RFs to contract (Fig 9C, light blue and purple, respectively). This shows that the principle of using STDP and traveling waves to refine a RF extends from 1D to 2D networks for simple, idealized waves. In the following section, we will examine how this process holds up when waves follow more irregular and complex trajectories.

Refinement of receptive fields by complex wave patterns
To test the robustness of refinement under conditions in which waves are far from idealized plane waves moving with a constant velocity, we used a wave model developed by Feller et al. [39] to generate complex wave patterns in the input layer. In the Feller et al. model, waves are generated spontaneously in random locations, and propagate along winding trajectories on a 2D input layer (see Methods). Due to the complexity of the model, it was not possible to set a precise wave speed. We therefore controlled the mean wave speed by temporally rescaling precomputed wave patterns, and measuring the speeds of waves that were isolated by a center of mass (COM) tracking algorithm (Methods and S7 Fig). In this way, we generated slow, medium and fast waves with speeds of 1.29±0.01 mm/s, 2.58±0.02 mm/s and 3.87±0.03 mm/s (mean ± SEM), respectively. Examples of two isolated waves are shown in Fig 10A, and a 90 s movie of isolated waves is provided in S1 Movie.
Using these rescaled wave patterns as input to simulations of RF development, we found that, quite remarkably, complex waves could shape and refine RFs in much the same manner as simple plane waves. With an asymmetric STDP rule and holding t + = 20 ms fixed, RFs could be made to expand (Fig 10B, left), contract (Fig 10B, center) or split into subfields (Fig 10B,  right) by the fast, medium and slow complex waves, respectively. Although slow waves did not always split the RF, the final RF area was always smaller with slow waves than with medium wave speeds. The evolution of the RF area during the simulation, averaged over repeated trials, is plotted in Fig 10C for the different wave speeds (solid: fast waves; dashed: medium waves; dotted: slow waves). Similarly, STDP rules with shorter STDP time scales produced RFs with smaller areas (Fig 10D and 10E), corresponding to a shorter characteristic wavelength in a periodic connectivity pattern. These results recapitulate the dependence of RF development on wave speed and STDP time scales, even when the input wave patterns followed complex and noisy trajectories.

Discussion
Here we examined how spatiotemporal correlations in traveling waves, a prominent feature of activity in the developing nervous system, could in principle interact with Hebbian STDP to instruct neural circuit development. Building upon previous studies of the interaction of STDP with correlated activity [16][17][18][19][20][21], we have identified a novel process by which traveling waves are able to establish and refine highly structured connectivity patterns in neural circuits. Specifically, our analysis has led to the following major results: 1) common examples of experimentally observed STDP build periodic patterns into feedforward networks when driven by traveling waves; 2) the spatial frequency of the periodic pattern scales with the wave speed and characteristic time scale of the STDP rule; 3) robust periodic connectivity patterns are produced only when the duration of the wave-related burst falls within a range that depends on the time scale of the STDP rule; the longer the burst duration, the longer the time scale of STDP must be to yield robust periodic patterns; 4) periodic patterning is robust to variability in the wave speed and wave frequency, but is lost if the ratio of non-wave spikes to wave-spikes becomes too high; 5) the shape and size of a RF can be modified by the spatial frequency imposed upon it by the wave-STDP interaction. Below, we discuss in greater detail how our results build upon previous work, and reflect on their application to experimentally observed wave phenomena and circuit refinement. We point out predictions from our model that may be tested using existing experimental paradigms, and suggest how our results may be extended in future work.

Pattern formation in the brain and biology
Our theoretical results, expressed in Eqs 6 and 9, are analogous to those derived by Swindale [35], who modeled the development of ocular dominance columns in primary visual cortex. Swindale showed that short range excitatory and long range inhibitory lateral interactions could give rise to the spatially periodic dominance of eye specific afferents across primary visual cortex, a principle also used to model the development of orientation columns in visual cortex [57]. In our model, the lateral interactions between synapses that are necessary for pattern formation result from the timing dependence of STDP, which is mapped onto space by the spatiotemporal correlations of traveling waves. Pattern formation of this kind is analogous to a Turing instability in reaction-diffusion systems [29], which has been applied to diverse cases of biological pattern formation [31]. A useful recipe for Turing-like pattern formation includes the presence of a diffusible activator and inhibitor, whereby the inhibitor diffuses over greater distances than the activator [30]. STDP rules emulate this feature if τ − > τ + , i.e. the temporal window for synaptic weakening is longer than that for synaptic strengthening. Such STDP rules are widely reported throughout the brain [14,58]. In the alternative scenario, when τ + > τ − , we found that STDP tends to weaken local regions of synapses that had strengthened by chance, preventing islands of strong synapses from forming. It is not impossible, however, for periodic patterns to form when τ + > τ − , so long as we carefully choose A + and A − , which control the overall bias for strengthening and weakening. However, we found that the stability of the pattern is very sensitive to small changes in A + and A − in this case.

Influence of input correlations and synaptic competition on refinement
Models of circuit development by Hebbian plasticity require constraints on the synapses so as to keep their strengths within biologically realistic bounds. A suitable choice of constraint, such as subtractive normalization, provides competition between synapses so that when a subset of synapses strengthen, all other synapses are suppressed [11,59,60]. Under this constraint, correlations between nearby inputs over short temporal windows encourage a localized group of synapses to strengthen, yielding a RF-like connectivity pattern [11]. This is because synapses separated by larger distances are uncorrelated within short temporal windows, so are rarely strengthened but are often weakened together. STDP rules that are biased for synaptic weakening can achieve this type of competition between uncorrelated synapses, and are therefore capable of building RF-like connectivity patterns [18]. However, this competitive process is fundamentally different to the mechanism by which STDP builds structured connectivity patterns in this paper. We recognize that a single wave can induce strong correlations between any two input neurons that it passes, irrespective of their separation. They will be correlated with a specific time lag, Δt, which scales with the separation, Δx, between the two inputs as Δx = vΔt. If Δt corresponds to synaptic weakening in the STDP rule, the synapses of the two input neurons will compete as a result of their strong correlations. Note that waves cannot produce periodic connectivity patterns if synaptic weakening does not have a particular dependence on timing, for example when competition is provided by subtractive normalization, as this will not yield the bandpass filtering that we have demonstrated for experimentally observed STDP rules (Fig 2B). This underlies why the wave speed, and the range of Δt values for which the STDP rule specifies synaptic weakening or strengthening, determines the resulting spatial frequency of the connectivity pattern, and therefore influences the type of modification experienced by a RF.
In Fig 5, we showed how contributions to the correlation function that are not wave-related act to degrade the periodic connectivity pattern. The strength of wave-related correlations may be further suppressed by recurrent inputs arriving from other postsynaptic cells, which we have excluded from our model. Output spikes that are driven by recurrent activity are likely to be much less correlated with spikes in the input layer and will therefore contribute an approximately constant term to the correlation function, C(x,Δt), the effect of which will be to emphasize any bias in the STDP rule towards synaptic weakening or strengthening. Moreover, excitatory recurrent connections are well known for encouraging similarities between postsynaptic RF properties, whereas inhibitory connections are known to decorrelate RFs, without drastically changing the basic RF properties of individual postsynaptic neurons [18,57,[61][62][63].
Increasing the bias for synaptic weakening in our simulations enhanced the degree to which RFs could be refined to a smaller size (Fig 8E). A bias for weakening in STDP is frequently observed in experiments [14,64], but is not the only possible source of synaptic competition, which can also be achieved by homeostatic regulation and intrinsic plasticity [65,66]. As mentioned above, these additional sources of competition would not aid pattern formation if they do not have the appropriate spike timing dependence that yields bandpass filtering as do the STDP rules used in this study.

Previous STDP models with space-time inseparable inputs
Several studies have previously explored the interaction between space-time inseparable input patterns and STDP. They demonstrated how motion stimuli can act through asymmetric STDP rules to impart an asymmetry in the spatial profile of excitatory connection strengths. This can be utilized to endow a network with direction selectivity similar to that found in the visual cortex [25][26][27][28]. Through our analytical results, we have shown that this phenomenon is one component of a richer set of dynamics with which the connectivity pattern can evolve under space-time inseparable inputs. Specifically, synaptic strengths are modified according to the convolution of the spatially mapped STDP rule, κ(x), with the synaptic strengths, w(x). Because of this convolution, the imaginary component ofkðkÞ results in spatial shifts, whereas the real component results in the emergence of a periodic connectivity pattern. It is possible that periodic patterning was occluded in previous modeling studies due to the setup of those models for the specific application to the development of direction selectivity.

Patterning and refinement by retinal waves
Retinal waves have been reported in the developing retina of several species [4,[67][68][69][70][71], suggesting that they play an important developmental role that has been conserved over the course of evolution. In mice, for example, disrupting normal retinal wave propagation [42,72] has a profound effect on the refinement of retinal ganglion cell (RGC) afferents to, and RFs in, superior colliculus (SC) [73,74]. One way in which retinal waves may drive this refinement is by imposing a wavelength onto the spatial structure of retinocollicular connections. To illustrate the relationship between RF sizes in SC, retinal wave speeds and STDP time scales, we computed the landscape of dominant spatial frequencies (S5 Fig) as a function of v and τ + , as in Fig 3, but using a typical retinal wave burst duration of 2 s for α(t) [42,43]. Retinal waves travel at relatively slow speeds of 0.1-0.2 mm/s, such that if the STDP rule at retinocollicular synapses is asymmetric with τ + = 20 ms, the characteristic wavelength would be 1/k Ã % 0.09 mm (solid red rectangle in S5 Fig). By interpreting a RF as half a cycle in the periodic pattern, a wavelength of 0.09 mm corresponds to a RF size in mice of 1.8°, using the distance to visual angle conversion in Methods. This is much smaller than the reported RF sizes of~10° [74]. In order to obtain RFs~10°in diameter (equivalent to a wavelength of 0.51 mm) with slow retinal waves, the STDP time scale would have to be τ + = 0.15 s (dashed red rectangle in S5 Fig), almost an order of magnitude longer than time scales commonly reported in STDP studies. The long duration of retinal wave bursts imposes an additional constraint on the minimum STDP time scale. We showed in Fig 4 that, when τ + = 20 ms, periodic patterning degraded when the burst duration exceeded 0.5 s. However, as the STDP time scale increases, so too should the limiting burst duration (Fig 4F). Assuming that the limit imposed by the burst duration scales linearly with τ + suggests that a STDP time constant exceeding 0.1 s should operate at retinocollicular synapses during retinal waves.
Our prediction of long STDP time scales at developing retinofugal synapses supports previous work by Butts and Rokhsar [41], who showed that retinal waves convey the most information about the relative retinotopic positions of RGCs over time scales ranging 0.1-2 s. This was further supported by the observation of such a rule at developing retinogeniculate synapses in the rat [44]. Evidence for long STDP time scales was recently provided in the visual cortex of mice that had a greater than normal expression of NMDA receptors (NMDARs) with NR2B subunits [75]. The STDP rule in these animals was temporally asymmetric, but was sensitive to remarkably long spike time differences greater than 0.1 s. NR2B expression is high during the stage of development when retinal waves are important for refinement, which in mice is the first postnatal week [52]. Thereafter, the number of NR2B containing NMDARs reduces, being replaced by faster acting NR2A containing NMDARs [76,77]. The change in NMDAR composition may explain why many STDP studies in mice, which are typically conducted in the second or third postnatal week of development, reveal STDP rules with time scales not much longer than 20 ms. We speculate that such developmental changes in NMDAR composition may provide a gradual reduction in STDP time scales that enables greater levels of refinement, as we demonstrated in Fig 8B-8C. The effect of this developmental progression may be complemented by the observed reduction in wave speeds in mice during the same period [43]. Our results describe a potentially important role for these phenomena in refinement, which can be examined experimentally.
The noisy mechanism of wave generation in the developing retina adds considerable variability to the size of retinal waves [39]. In our model, it is important that the total distance travelled by a wave exceeds the characteristic wavelength, 1/k Ã . This is the spatial analog of requiring the wave duration to exceed the time scales of STDP in our analytical derivation. Our calculation above therefore suggests that retinal waves should travel distances greater than 0.51 mm in order to produce RFs 10°in diameter in the mouse SC. By tracking the center of mass of experimentally recorded retinal waves, Maccione et al. showed that a substantial majority of retinal waves indeed exceeded this distance throughout the period of major retinocollicular refinement [43].
Taken together, experimentally recorded retinal waves exhibit many of the properties required for retinocollicular refinement by a periodic patterning process. However, whether the mode of synaptic plasticity at retinocollicular synapses is suitable for this process has yet to be determined.

Suitability of traveling waves for refinement in other brain areas
Traveling waves occur in many different parts of the brain, and each location will have particular constraints that limit the types of patterns that could form based on the interaction with STPD rules. As discussed above, one constraint is the characteristic wavelength, 1/k Ã , which sets a lower bound to the distance traveled by waves if periodic patterning is to be achieved.
We computed the characteristic wavelengths associated with waves in brain areas other than the retina-including the cortex, cerebellum and hippocampus-and list them in Table 1, under the assumption of an asymmetric STDP rule like that used in many of the results above (τ + = 20 ms), and which is reported widely throughout the brain. The speeds of waves in these areas are 1-2 orders of magnitude faster than retinal waves [5,8,9], but yield characteristic wavelengths that are similar in scale to that required for RFs in SC because of the short STDP time scale (first four entries in Table 1). Fast wave speeds and short STDP time scales place these waves in a region near the bottom right hand corner of Fig 3A (solid red rectangle), assuming a burst duration of 0.1 s. Also listed in Table 1 are the maximum distances over which the waves in different areas may propagate. Comparing these distances with the characteristic wavelengths reveals that periodic patterning in the efferents of the cerebellum, dorsal cortex and ventral cortex would not be expected, as the maximum possible propagation distance is close to or less than the lower bound set by 1/k Ã . However, wave properties in the hippocampus, as in the retina, do satisfy the minimum distance constraint for the type of STDP rule assumed. More detailed analyses of waves in the hippocampus and STDP at its efferent synapses is necessary to determine whether these phenomena could drive periodic patterning in a manner useful for development.

Development of cortical receptive fields
RFs of simple cells in the primary visual cortex (V1) exhibit weak selectivity for orientated features in the visual environment at the onset of vision [78][79][80] and become more selective with visual experience [79,81]. Yet when an animal is raised in an environment comprising a restricted range of oriented contours, RF selectivity matures for those same orientations but not others [79][80][81][82][83][84]. A key feature of orientation tuned simple cell RFs is their composition of ON and OFF subfields, which are respectively sensitive to increments and decrements in luminosity, and are thought to reflect the visuotopic organization of ON and OFF inputs from the lateral geniculate nucleus (LGN) [85,86]. The orientation of these subfields is thought to confer the orientation preference of the simple cell [85,86]. However, it is not clear whether oriented features in the environment instruct the development of orientation selectivity, or whether they permit the maturation of RFs that are already selective for the same orientations. The spatial frequency tuning of simple cells may help to resolve this question, as it is well predicted by the periodic structure of ON and OFF subfields [87]. One model of simple cell RF development [57] posits that spatial frequency tuning is determined by the structure of spatial correlations between LGN cells. However, the necessary correlations were not observed in the developing LGN in later experiments [88]. Other models inspired by sparse coding schemes [89] posit that simple cell RFs result from learning the independent components of natural visual scenes [90], a consequence of which is that RFs would exhibit tuning for velocity and spatial frequency with an inversely proportional relationship [91]. However, despite recent insights [23,63], it is not yet well established how sparse coding is implemented in biological circuits. Table 1. Summary of different wave phenomena. Wavelengths, 1/k*, were computed using the wave speeds given, a burst duration of 0.1 s, and the same asymmetric STDP rule for the first four rows, with τ + = 20 ms, τ − = 40 ms, A + = 1.0 and A − = 0.51. For retinal waves in the last row, we used a burst duration of 2 s, and an asymmetric STDP rule with τ + = 0.2 s, τ − = 0.4 s, A + = 1.0 and A − = 0.55.

Brain area
Wave speed (mm/s) 1/k* (mm) Max. propagation distance (mm) Our results provide an alternative mechanism for simple cell RF development. We have shown how oriented wavefronts can build orientated RFs (Fig 9A), and how the same process can develop multiple, periodic subfields (Fig 7 and S4A Fig), akin to the ON and OFF subfields of simple cells, thus providing a means for spatial frequency tuning (Fig 8). Moreover, our prediction of an inverse relationship between the wave speed and spatial frequency of the connectivity pattern provides a novel test for the role of visual experience in shaping simple cell development, and bares striking similarity to the sparse coding schemes mentioned above. Using stimuli that consist of high contrast luminance contours to elicit wave-like activity in the retina and LGN may help to test our predictions, and experimental protocols using such stimuli with young animals are already well established (for example, see [92]).
We can use our model to predict a range of plausible STDP parameters that would build RFs with the spatial frequency tuning of simple cells. In adult cats, spatial frequency tuning ranges from~0.2-2 cycles/°within eccentricities of ±15° [54,93]. We further constrain the model by considering wave speeds that match the velocity tuning of cat simple cells (~0.5-20°/s, [94]). In S6 Fig, we illustrate the spatial frequencies that are obtained in selected parts of this large STDP parameter space. Assuming an asymmetric STDP rule in the cortex [36,95], realistic spatial frequencies can be obtained with STDP time scales ranging~1-100ms, and biases for synaptic weakening in the approximate range 0.3 ≲ A -/A + ≲ 0.7. Realistically, STDP need not operate over such a wide range of parameters, as variability in the speeds of natural stimuli should explain most of the variability in spatial frequency tuning.
The relationship between temporal correlations and the development of spatially structured connectivity patterns may be extended to other visual RF properties. For example, temporal delays between inputs with spatially offset RFs are the major components needed to build direction selective cells [96,97]. This kind of organization, and hence direction selectivity, can be learned with rate-dependent Hebbian plasticity [62], utilizing diverse response latencies in the LGN [98][99][100]. It is feasible that STDP should also yield spatial offsets between inputs with different response latencies, given its Hebbian nature. Combined with the capacity for periodic patterning, we speculate that the interaction of traveling waves with STDP could yield connectivity patterns that are direction selective, as well as orientation and spatial frequency tuned.

Tolerance to background spiking noise
We found that pattern formation in our simulations was sensitive to noise: patterns began to degrade when the ratio of the background spike rate to the wave-induced spike rate became too high. It is therefore essential to consider what is known about background firing rates reported in the literature. We concentrate on background noise in the developing visual system, for which good data are available. During retinal wave activity in mice, retinal ganglion cells rarely spike outside of a burst. However, bursts that occur outside of a wave event have approximately the same firing rate as wave related bursts ( Table 1 in [42]). Nevertheless, wave related bursts comprise approximately 90% of all bursts in the developing retina [42], suggesting that background spiking noise would have little impact on wave-induced correlations at this stage of visual development. Retinal waves also drive bursts in the LGN at firing rates of 10 Hz with long intervening periods of quiescence [101]. In later stages of development, visually evoked responses in the LGN reach firing rates of 10 Hz, whereas 1 Hz firing rates are typical of spontaneous activity [102]. This highly skewed distribution of firing rates is recapitulated in the developing visual cortex [103]. Thus, it is likely that, during early developmental periods when RF properties undergo refinement, spontaneous spikes contribute little to the overall activity in the visual system, and that patterned retinal waves and visual stimulation are propagated throughout the early visual system.

Modulation of STDP by other factors
The precise dependence of STDP on spike timing can be sensitive to many factors that can modify the shape of the STDP rule (reviewed extensively in [14] and [58]). Whatever the mechanisms that underly these additional interactions, they must preserve a bandpass dependency on spike times in order for periodic patterns to form. Biophysical models of cellular activity that may underlie STDP have been proposed to account for some of these additional interactions [104,105]. These models rely on the summed pre-and postsynaptic contributions to the intracellular calcium concentration, the amplitude of which at any given time determines whether synapses are strengthened or weakened. As such, the shape of the STDP rule changes for different spike trains and thus would not adhere to the spike-timing dependence that is necessary for bandpass filtering and pattern formation. If pattern formation as outlined in this paper does influence circuit development, then we predict that STDP must maintain a bandpass profile at synapses that mediate wave-like activity patterns. To test this, a typical STDP experimental protocol could be performed, in which paired spike bursts are stimulated in connected pre-and postsynaptic cells over a range of temporal offsets, and the resulting change in synaptic strength measured. A key requirement of such an experiment would be to extend the temporal offsets well beyond the burst durations, in contrast to previous studies [44,106,107], so as to detect the full temporal profile of the plasticity rule and ensure that synaptic changes decay to zero with larger offsets. If the bandpass property of synaptic plasticity is present at these synapses, then it should be reflected in the resulting STDP curve. Retinocollicular or geniculocortical synapses would be ideal substrates for testing this.

Future directions
The robustness of our results to various types of noise suggests that pattern formation should also be achieved in more complex models that incorporate details specific to the circuit. The developing visual system provides a promising arena in which to test our theory of pattern formation for the reasons discussed above, and we outline here future work that will facilitate this investigation.
The dominant spatial frequency that characterizes a periodic pattern is strongly influenced by the stimulus response properties of the input neurons, as demonstrated by the effect of the burst duration in Fig 4. Incorporating accurate space-and time-dependent response properties of neurons in the visual pathway, as determined by their spatiotemporal RFs [55,[108][109][110][111], is therefore essential to making accurate predictions about pattern formation in downstream visual targets such as the SC and V1. For theoretical analysis, these features can be easily incorporated by reformulating α(t) as α(x,t) to take account presynaptic RF structures. Application of our model to simple cell RF development will require modeling both ON and OFF response types in the input layer, which must become spatially segregated as a result of learning. Previous proposals for the mechanism of this segregation typically rely on correlations within each ON or OFF population exceeding those between the two populations [20,57,88]. The extent to which STDP will be sensitive to ON and OFF correlations, in addition to the strong spatiotemporal correlations induced by traveling waves, requires further investigation.
The linearity of our current model, though useful for theoretical analysis, prevents RFs from developing an orientation preference when waves travel in multiple directions (Fig 9B), as is the case for retinal waves and natural visual environments. Any bias towards one orientation is eventually averaged away by the influence of waves with other orientations. The addition of a nonlinearity, either in the sensitivity of STDP to the output firing rate, or in the transfer function for the output firing rate itself, should help to break the symmetry in wave directions and bias the RF towards any asymmetric shape that is built into it by chance. In this way, a population of output neurons might be able to acquire different orientation preferences when exposed to the same input pattern. In a similar fashion, a spectrum of spatial frequency preferences could be acquired throughout the population, rather than the uniform preference for a single spatial frequency, as we observed in simulations in which the wave speed was varied.
It is worth noting that traveling waves are just one incarnation of activity patterns that exhibit the space-time inseparable correlations necessary for periodic patterning. Our results may be generalized to any reliable temporal sequence of activity, some examples of which include navigation codes in the hippocampus [112], pre-motor coding in songbirds [113], and spontaneous 'synfire chains' in the cortex [114].

Overview
We model a reduced feedforward network consisting of a presynaptic layer of input neurons and a single postsynaptic output neuron. Traveling waves of activity traverse the input layer and recruit input neurons, which discharge bursts of spikes and provide excitatory synaptic inputs to the output neuron. Each feedforward synapse is modified according to the time delay between spikes from its corresponding input neuron and spikes from the output neuron. A schematic of the network is provided in Fig 1A. Spiking neuron simulation Our model is simulated using a time step of 1 ms. During each step, a number of processes are simulated. First, input neurons generate spikes from one of two wave models described below. Each input spike elicits in the output neuron an excitatory postsynaptic potential (EPSP), which is weighted by the synaptic strength. EPSPs modulate the output membrane potential, which determines the generation of output spikes according to either a linear or nonlinear process, described below. Once all spikes have been generated for a given time step, synaptic strengths are modified by a STDP rule.
Traveling waves in the input layer.
1D and 2D plane waves The input layer comprised point neurons, arranged 20 μm apart in a 1D chain or 2D square lattice. Plane waves were always initiated at one edge of the input layer and travelled the entire length to the opposite edge. For most simulations using plane waves, a blank period lasting at least 5 s was inserted between the end of one wave and the beginning of the next so as to allow the postsynaptic activation to decay to zero before the next wave traversed the inputs. The only simulations in which waves were not interleaved with a blank period were those in which the presence of multiple waves was tested.
In most 1D simulations, the wave direction was alternated between each wave. In simulations with multiple waves, the wave direction was alternated after every 100 waves. In the 2D plane wave simulations, the wave direction was selected randomly out of N possible directions that were equally spaced around the compass, where N was either 2, 4, 6 or 16.
When an input neuron is recruited by the wavefront, it undergoes a burst of spikes. The input firing rate during a burst, α(t), is modeled as a boxcar function: where H(t) is the Heaviside step function, R in is the input firing rate, and d is the burst duration. During an input burst, spikes were generated in time steps of 1 ms using a Bernoulli process. Unless otherwise stated, R in = 50 Hz.
In some simulations, input neurons were given background firing rates, R bg in , which are stated in the relevant results sections. In these simulations, the background rate was not added to the firing rate of inputs undergoing a wave-related burst. Otherwise, when an input neuron was not being recruited by a wave, it was silent.
In all cases, we repeated trials of plane wave simulations with the same model parameters 16 times, using different seeds for the random number generators, such that each trial would produce different spike trains. In simulations where a parameter that was normally fixed was instead a random variable, such as the wave speed, or wave direction in the case of 2D plane waves, the random number generator for that random variable was reseeded for repeated trials. Otherwise, the random number generator for spike generation was reseeded. In all cases, original spike sequences were always generated for each trial.
Complex wave patterns To generate more complex wave patterns in 2D, we implemented a simulation based on a previously published model of spontaneous depolarizations that propagate across the developing mammalian retina [39]. In short, cells in the model that undergo spontaneous depolarizations provide excitatory input to nearby cells, triggering cascades of excitatory synaptic transmission that propagate throughout the network, while strong refractory currents quickly return active cells to the rest state. As such, wavefronts of activity propagate along winding trajectories between refractory domains of the network.
The model comprises a layer of starburst amacrine cells (SACs) and a layer of retinal ganglion cells (RGCs). Both layers are arranged on 64×64 square lattices with a cell spacing of 34 μm. SACs receive excitatory synaptic input of equal strength from other SACs that are within a radial distance of 120 μm. RGCs also receive excitatory input from SACs over the same distance but do not provide synaptic input to other RGCs or SACs. None of the synapses within the model retina undergo synaptic plasticity.
Every cell in the retina model is characterized by a dynamic excitation variable, X i , which is analogous to its membrane potential, for i = 1. . .N where N is the number of SACs or RGCs. The simulation evolves in time steps of 100 ms and X i is updated for cells in the SAC (S) and RGC (R) layers according to where n S i is the number of active SACs at time, t, that are connected to cell i (from a maximum of 38 connected SACs). Thus, active SACs increase the excitation of connected SACs and RGCs. In the absence of input from any SACs, the excitation variables for both SACs and RGCs decay exponentially with a decay constant of τ R,S = 0.1 s. A SAC becomes active when X S i exceeds a threshold, θ S = 6. Similarly, RGCs become active when X R i exceeds a threshold, θ R = 10. Synaptic transmission from a SAC persists for 1 s, after which the SAC becomes refractory. Every time a SAC becomes refractory, the refractory period is chosen from a normal distribution with mean τ ref = 40 s and standard deviation σ ref = 20 s. After the refractory period, X S i is set to zero. To model the spontaneous initiation of waves, SAC activation variables, X S i , were forced to exceed θ S at each time step with a probability of 0.0035.
There are several possible ways to manipulate the wave speed in the wave model by tuning different parameters [39,115]. However, this approach has the undesirable effect of changing the overall spatiotemporal structure of the waves, which complicates the interpretation of the effect wave speed has on receptive field development. In order to manipulate the wave speed without altering the wave patterns, we simulated long periods of wave activity and rescaled these activity patterns in time by a factor, F t . This involved recording the set of times, T i , at which the RGC activation, X R i , exceeded θ R . We then rescaled T i by F t to compress the wave patterns in time, and used the rescaled T i as burst onset times for the i th RGC. RGC spike bursts 0.1 s in duration were then generated using a Bernoulli process with a firing rate of 50 Hz and temporal precision of 1 ms. The RGCs in this model retina act as the input neurons in the final set of 2D simulations in this paper. To produce slow, medium and fast waves, we used F t = 0.21, 0.11 and 0.07, respectively. Thus, medium and fast waves were respectively two and three times faster than the slow waves.
We repeated trials of complex wave simulations with the same model parameters three times, using different seeds for the τ ref random number generator, such that each trial would produce a different sequence of wave patterns.
Excitatory postsynaptic potentials. Each input spike elicits an EPSP in the output neuron. If at time, t j,n , input j fires its n th spike, then at some later time t, the resulting EPSP is described by a difference of exponentials: where τ d = 5 ms and τ r = 1 ms are the EPSP decay and rise times, respectively, based on neurotransmission through AMPA receptors, and j,n (t) = 0 for t < t j,n . Linear output neuron. Spikes in the linear output neuron are generated by a Bernoulli process: at time, t, the output fires a single spike with a probability that is determined by its time varying firing rate, λ(t), which is proportional to the summed EPSPs: Here, R out is a constant of proportionality, chosen to constrain the output neuron to realistic firing rates in the range 10-100 Hz, and w j is the strength of the synapse from input neuron j at the time of its n th spike.
Nonlinear output neuron. We model a nonlinear output neuron with leaky integrate and fire (LIF) dynamics. The output neuron fires a spike at timet when its membrane potential, u(t), is pushed above a threshold, ϑ, by incoming EPSPs, (t). EPSPs follow Eq 17, which models the leaky response of u(t) to an inward synaptic current, I(t), that has an instantaneous rise and exponential decay: I(t) = H(t)exp(−t/τ r )/τ r , with H(t) the Heaviside step function. After an output spike, an instantaneous outward current, ÀWdðt ÀtÞ, returns u(t) to the resting potential, u rest = 0. The full response of u(t) to incoming EPSPs and the outward current follows the equation where ξ(t) is a refractory kernel that models the leaky response of u(t) to the outward current. The kernel ξ(t) is therefore a decaying exponential that imposes a relative refractory period ρ rel = τ d = 5 ms. In addition, we incorporate into ξ(t) an absolute refractory period, ρ abs = 2 ms, during which it is impossible for the output to spike. The full refractory kernel is therefore: Note that ξ(t) depends only on the time since the last output spike, whereas EPSPs summate over all previous input spikes.
Synaptic plasticity. After the spike generating process, synaptic strengths are updated using a spike-timing dependent plasticity (STDP) rule. Two rules are used in this paper, one asymmetric [13], K asym , and the other symmetric [38] in time, K sym ( Fig 1C): Here, Δt = t in −t out is the time difference between spikes belonging to input and output neurons, respectively. The size of the temporal windows for synaptic strengthening and weakening are scaled by τ + and τ − , respectively, and A + and A − scale the relative degree of strengthening and weakening, respectively. The size of synaptic modifications, or learning, is further scaled by η, which is set to be a small number. Every time an input (output) neuron spikes, all output (input) spikes in the previous 5τ − seconds are used to update a synapse. Unless otherwise stated, all asymmetric STDP rules are constructed such that τ − = 2τ + , and all symmetric STDP rules have τ − = 1.6τ + . In the results, we therefore almost always refer to the STDP time scale using just τ + . Furthermore, unless otherwise stated, A + = 1.0 and A − = 0.51 for all asymmetric rules, whereas A + = 3.2 and A − = 2.1 for all symmetric rules. In the simulations with complex wave patterns, an asymmetric rule is used with A + = 1.0 and A − = 0.55.
To prevent synaptic strengths from increasing ad infinitum, or from becoming negative and thus inhibitory, both of which are biologically implausible scenarios, we imposed hard bounds at synaptic strengths of 0 and 1.
Initial synaptic strengths. At the beginning of every simulation, the synaptic strengths were initialized with one of two configurations: 1) all synapses were given a strength of 0.5; 2) the central synapses, spanning a distance of RF 0 , were given the maximum strength of 1.0, and all other synapses given strengths of zero, to simulate a coarse RF structure. When the synaptic strengths were initialized with a RF, an arbor function, A(x), was applied to maintain the RF structure and prevent a periodic connectivity pattern from emerging across all synapses. The arbor function was implemented by multiplying A(x) with @ T w(x), where where a 0 is the radius of the arbor and x 0 is the center of the arbor. Local input correlations. In some simulations, we induced additional correlations between neighboring inputs by generating spikes in two stages for each time step. In the first stage, spikes were generated by a Bernoulli process as usual but with a modified input firing where c is the correlation strength and takes values between 0 and 1. In the second stage, correlated spikes were generated with a probability, c, in all cells neighboring inputs that spiked in the first stage. In simulations with background firing rates, the background rate was likewise modified so thatR bg in ¼ R bg in = 1 þ 2c ð Þ. Dividing by the factor 1+2c in the first stage kept the firing rates constant after the second stage, allowing independent manipulation of the correlation strength.
Wave speed and inter-wave interval distributions. In simulations with variable wave speeds or IWIs, a new speed or IWI was drawn from a lognormal distribution for each wave. For the i th wave, the lognormal random variable, y, was: are respectively the mean and standard deviation (SD) of ln(y), μ and σ are respectively the mean and SD of y, and z is a normally distributed random variable with zero mean and unit variance.
For simulations with variable speeds, μ v = 4 mm/s was fixed and σ v was varied between simulations. We imposed a minimum speed of 0.05 mm/s, such that the i th wave speed was: For simulations with variable IWIs, the IWI distribution was modified slightly to ensure that IWIs were no shorter than the burst duration, d. That is, successive waves never overlapped. Each new IWI was therefore: where y was computed using m(μ−d,σ) and s(μ−d,σ). For these simulations, σ = 1 s was fixed, and μ was varied between simulations. Because a downtime of at least 5 s is needed for waves to clear the input layer before switching direction, we switched the wave direction in these simulations after every 100 waves. That is,~99% of the waves sampled from the IWI distribution.

Analysis of complex wave patterns
To measure the speed and size of simulated complex wave patterns, we analyzed a 2000 s segment of RGC spiking activity from a simulation of the slow waves only, as these waves were less compressed in time and therefore lasted longer, enabling a more accurate measure of the wave properties. Spike times were first converted into a firing rate movie, M(x,y,t), with dimensions 64 × 64 × 20000, where the (i,j,k) th bin in M contained the number of spikes fired by a single RGC, at location (x i ,y j ), during a 100 ms period, t k . During the retinal wave simulation, multiple waves were occasionally present at the same time. Concurrent waves would occasionally collide into or split from each other, but they were mostly well separated in space. We sought to isolate concurrent waves to accurately measure their individual properties, such as speed and size, using a center of mass (COM) tracking algorithm, of which a schematic is drawn in S7 Fig. First, we computed a smoothed firing rate movie, M s (x,y,t) = (M Ã G)(x,y,t), where Gðx; y; tÞ ¼ e t is a Gaussian filter with a spatial standard deviation (SD), σ x = σ y = 34 μm, and temporal SD, σ t = 100 ms. After smoothing, pixel values below 10 were set to zero. Thus, a single time bin in M s might have several domains of non-zero pixels. To track the COM of different waves, a boundary was drawn around each domain of activity in time bin k and the COM within the boundary calculated, where mass refers to the firing rates within the same boundary in M(x,y,t = k). Two domains, colored blue and orange, and their COMs (purple and green, respectively) are shown in S7 Fig, and the firing rates given in greyscale. Each domain was assigned a wave identification number (ID) but, if the COMs at two domains were less than 680 μm apart, they were given the same ID. If the COM of one domain was within a 680 μm radius from the COM of another domain in the previous time bin (purple lines extending from the COM in S7 Fig), it was given the same ID. In this way, wave COMs were tracked according to their ID. The COM trajectory of each wave was zero-padded and smoothed in time by a Gaussian filter with a SD of 0.1 s, and the first two and last two COM locations discarded. Thus, wave speeds and sizes were only computed for waves that lasted for more than 0.5 s. We used the path length of a COM trajectory as the distance travelled by a wave, and computed the wave speed from this by dividing it by the time taken to travel that trajectory.
RGCs near the edges of the model retina received fewer lateral inputs from SACs and were less active than those in the centre. Within the central 44 × 44 RGCs, the total level of spiking activity was comparably uniform. Accordingly, waves with a time-averaged centre of mass (COM) that resided in the outer 10 neurons were discarded from any further analysis. Occasionally, small segments of waves were isolated from larger waves. These were discarded from the analysis by removing waves that covered fewer than 1000 space-time bins.
The COM tracking algorithm performed well in separating waves that eventually merged with, or split from, other waves, and allowed us to perform analysis on almost every wave that had a unique ID by removing concurrent waves from that analysis. A total of 797 waves were isolated in the 2000 s segment, of which 779 lasted long enough to compute the wave speed and distance travelled.
All analysis was performed with Matlab R2012a using built-in and custom built functions.

Computing numerical solutions
Concordance of the spiking model with analytical results is verified by numerically integrating Eq 6 using the forward Euler method. Synaptic strengths were initialized in the same way as in the spiking simulation. However, as there was no spiking noise when solving Eq 6 numerically, low amplitude Gaussian white noise was added to the initial synaptic strengths. To further align numerical integration with the simulations, nonlinearities that were present in the simulations were also incorporated in the numerical integration by: 1) restricting synaptic strengths to the range [0, 1]; 2) alternating the direction of the wave by replacing the wave speed, v, with (−1) n v for the n th iteration; 3) applying the arbor function in Eq 23 when the initial condition of w(x) supported a RF structure. To prevent numerical solutions from becoming chaotic, the learning rate, η, was varied depending on the values of τ + and v.

Robustness measures
To measure the robustness, C w , with which a periodic structure emerged in a connectivity pattern, we computed the discrete power spectrum, P(k i ), of the connectivity pattern, with the DC component removed, and determined the ratio of the power at the dominant spatial frequency to the total power: where k N is the Nyquist limit. As a predictor of robustness in the connectivity patterns, a theoretical robustness measure, C κ , was computed for the power spectrum of the real component ofkðkÞ: j Re½kðkÞj 2 , which was obtained by first numerically computing κ(x) with a spatial resolution of 0.6 μm. The theoretical robustness was therefore: Converting between visual angle and retinal distance To set up simulations of RF refinement, we use data from experimental studies regarding receptive fields, or traveling wave speeds, for example, which require conversion between degrees of visual angle and units of distance along the retinal surface. The precise conversion relationship between visual angle and retinal distance depends on the geometry of the eye, which differs for different animals. We therefore use a general conversion factor to provide a useful estimate. We assume that the retina covers two thirds of the circumference of a horizontal section through the center of the eye [116], and that this maps to 180°of visual angle. This means that distance in millimeters, R, corresponds to visual angle, A, according to R = πrA/ 135, where r is the radius of the eyeball. We used the following values for the radii of eyes in different animals at different stages of development: mouse, 1-7 days old: r % 1.1 mm ( Table 1 in [117]); mouse, adult: r % 1.7 mm ( Table 1 in [117]); ferret, eye opening: r % 2.8 mm (Fig. 8 in [118]); cat, 4 weeks old: r % 5 mm (Fig. 3 in [119]).
Supporting Information S1 Movie. A 90 s extract of simulated slow waves demonstrate the performance of the COM tracking algorithm. Black and colored pixels correspond to RGCs that fired at least one spike in a 100 ms time bin. Each isolated wave is assigned a new color. Activity that was not assigned to any wave is in black.
(AVI) S1 Text. Full derivation of Eq 5 in the main text. , where α(t) has been modeled using a smooth function: aðtÞ ¼ t d e Àt=d . All other parameters are the same as for Fig 4D. The power is distributed across a broad range of high frequencies for a 0.01 s burst, and within a small band of the lowest frequencies for a 1.0 s burst, yielding the strong negative lobe in the black curve. The negative lobe in the black curve extends beyond the horizontal axis and has been cut for clarity. However, for 0.1 s bursts, the power is concentrated in between these two extremes, around the dominant spatial frequency, k Ã . (TIFF) for burst durations of 2 s. Solid black contours denote spatial frequencies equal to 10 raised to integer exponents. Spatial frequencies were obtained by locating the maximum inkðkÞ as a function of v and τ + , using an asymmetric STDP rule and a burst duration of 2 s for α(t). Sharp transitions in spatial frequency along the τ + axis are due to κ(x) having several peaks of near equal amplitude (c.f. black curve in Fig 4D), such that small changes in τ + can change which peak is the global maximum. Solid red rectangle: given a typical STDP rule with τ + = 20 ms, the connectivity pattern associated with retinal wave speeds would have a dominant spatial frequency of~11 cycles/mm, or a wavelength of 0.9 mm. Dashed red rectangle: RFs in the SC require a characteristic wavelength of~0.51 mm, which corresponds to a spatial frequency of~2 cycles/mm. Given the speed of retinal waves, the required STDP time scale is predicted to be 0.1-0.2 s. (TIFF) S6 Fig. Spatial frequency and velocity tuning of cat simple cells dictate a plausible range of STDP parameters at geniculocortical synapses. Using a burst duration of 100 ms, which matches the impulse response duration of immature LGN cells in the kitten [111], we computed a spatial frequency map as a function of the STDP decay times, τ + and τ -, for four conditions: the STDP amplitudes were either A − /A + = 0.3 (top row) or A + /A − = 0.7 (bottom row), and the wave speed was either 3°/s (left column) or 10°/s (right column) [94]. We further restricted our analysis to STDP rules that had a reasonable bias for either weakening or strengthening by ignoring any rule for which the DC power exceeded that at the dominant spatial frequency, i.e. cases when jkð0Þj 2 > jkðk Ã Þj 2 . We also restrict the spatial frequency maps to frequencies that lie in the range observed in adult cats: 0.2-2 cycles/° [54,93].
(TIFF) S7 Fig. Schematic for COM tracking algorithm. The two images depict activity that was generated with the complex wave model at time bins k (left) and k + 1 (right). The solid blue areas mark the domains in M s (x, y, t) that were assigned to one isolated wave, and the solid orange areas mark domains that were assigned to another isolated wave. Greyscale pixels illustrate the firing rates of RGCs in M(x, y, t) within each domain. Solid purple and green dots denote the COM for the first and second wave, respectively, in the current time bin. To illustrate how the COMs moved between time bins, the open purple and green dots (right) denote the COM of each wave in the previous time bin (left). Domains with COMs that are separated by less than 680 μm (purple circle around the blue wave), in the same time bin or in adjacent time bins, are assigned to the same wave. (TIFF)